বিশেষ ম্যাট্রিক্স — ট্রাইডায়াগোনাল ও স্পার্স সিস্টেম

class OpCounter:
    def __init__(self):
        self.n = 0
    def mul(self, a, b):
        self.n += 1
        return a * b
    def div(self, a, b):
        self.n += 1
        return a / b
    def addsub(self, a, b, sign=1):
        self.n += 1
        return a + sign * b

def make_tridiagonal_system(n):
    # ডায়াগোনালি-ডমিন্যান্ট ট্রাইডায়াগোনাল সিস্টেম -- a=নিচের, d=প্রধান, c=উপরের ডায়াগোনাল
    a = [0.0] + [1.0] * (n - 1)
    d = [4.0] * n
    c = [1.0] * (n - 1) + [0.0]
    b = [float(i + 1) for i in range(n)]
    return a, d, c, b

def thomas_solve(a, d, c, b, counter):
    n = len(d)
    cp = [0.0] * n
    dp = [0.0] * n
    cp[0] = counter.div(c[0], d[0])
    dp[0] = counter.div(b[0], d[0])
    for i in range(1, n):
        m = counter.addsub(d[i], counter.mul(a[i], cp[i - 1]), sign=-1)
        if i < n - 1:
            cp[i] = counter.div(c[i], m)
        dp[i] = counter.div(counter.addsub(b[i], counter.mul(a[i], dp[i - 1]), sign=-1), m)

    x = [0.0] * n
    x[n - 1] = dp[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = counter.addsub(dp[i], counter.mul(cp[i], x[i + 1]), sign=-1)
    return x

def tri_to_dense(a, d, c, n):
    A = [[0.0] * n for _ in range(n)]
    for i in range(n):
        A[i][i] = d[i]
        if i > 0:
            A[i][i - 1] = a[i]
        if i < n - 1:
            A[i][i + 1] = c[i]
    return A

def gaussian_solve_counted(A, b, counter):
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for k in range(n - 1):
        for i in range(k + 1, n):
            factor = counter.div(M[i][k], M[k][k])
            for j in range(k, n + 1):
                M[i][j] = counter.addsub(M[i][j], counter.mul(factor, M[k][j]), sign=-1)
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = M[i][n]
        for j in range(i + 1, n):
            s = counter.addsub(s, counter.mul(M[i][j], x[j]), sign=-1)
        x[i] = counter.div(s, M[i][i])
    return x

print(f"{'n':>4} | {'Thomas অপারেশন':>16} | {'সাধারণ Gaussian অপারেশন':>24} | {'অনুপাত':>10}")
for n in [5, 10, 20, 50, 100]:
    a, d, c, b = make_tridiagonal_system(n)
    A = tri_to_dense(a, d, c, n)

    ct1 = OpCounter()
    x_thomas = thomas_solve(a, d, c, b, ct1)

    ct2 = OpCounter()
    x_gauss = gaussian_solve_counted(A, b, ct2)

    max_diff = max(abs(x_thomas[i] - x_gauss[i]) for i in range(n))
    ratio = ct2.n / ct1.n
    print(f"{n:>4} | {ct1.n:>16} | {ct2.n:>24} | {ratio:>9.2f}x   (সমাধান পার্থক্য সর্বোচ্চ {max_diff:.2e})")

# n=5 কেসে মূল সমীকরণে বসিয়ে যাচাই
n = 5
a, d, c, b = make_tridiagonal_system(n)
A = tri_to_dense(a, d, c, n)
ct = OpCounter()
x = thomas_solve(a, d, c, b, ct)
print()
print(f"n={n} সমাধান (Thomas):", [round(v, 6) for v in x])
print("যাচাই -- মূল Ax=b-তে বসিয়ে অবশিষ্টাংশ:")
for i in range(n):
    lhs = sum(A[i][j] * x[j] for j in range(n))
    print(f"সমীকরণ {i+1}: LHS = {lhs:.10f}   b = {b[i]:.6f}   residual = {lhs - b[i]:.2e}")

   n |   Thomas অপারেশন |  সাধারণ Gaussian অপারেশন |     অনুপাত
   5 |               33 |                      135 |      4.09x   (সমাধান পার্থক্য সর্বোচ্চ 5.55e-17)
  10 |               73 |                      895 |     12.26x   (সমাধান পার্থক্য সর্বোচ্চ 4.44e-16)
  20 |              153 |                     6290 |     41.11x   (সমাধান পার্থক্য সর্বোচ্চ 8.88e-16)
  50 |              393 |                    89475 |    227.67x   (সমাধান পার্থক্য সর্বোচ্চ 1.78e-15)
 100 |              793 |                   691450 |    871.94x   (সমাধান পার্থক্য সর্বোচ্চ 7.11e-15)

n=5 সমাধান (Thomas): [0.167949, 0.328205, 0.519231, 0.594872, 1.101282]
যাচাই -- মূল Ax=b-তে বসিয়ে অবশিষ্টাংশ:
সমীকরণ 1: LHS = 1.0000000000   b = 1.000000   residual = 0.00e+00
সমীকরণ 2: LHS = 2.0000000000   b = 2.000000   residual = 0.00e+00
সমীকরণ 3: LHS = 3.0000000000   b = 3.000000   residual = 0.00e+00
সমীকরণ 4: LHS = 4.0000000000   b = 4.000000   residual = 0.00e+00
সমীকরণ 5: LHS = 5.0000000000   b = 5.000000   residual = -8.88e-16