# gauss.py - een stelsel oplossen met de methode van Gauss (met rijverwisseling), en vergelijken met numpy
import numpy as np


def gauss(A, b):
    """Lost A x = b op. Geeft x en het aantal vermenigvuldigingen en delingen terug."""
    A, b = np.array(A, float), np.array(b, float)
    n = len(b)
    bewerkingen = 0
    for k in range(n - 1):                          # heen: nullen onder de diagonaal
        p = k + np.argmax(abs(A[k:, k]))            # grootste spil in kolom k
        if A[p, k] == 0:
            raise ValueError("de matrix is singulier")
        A[[k, p]], b[[k, p]] = A[[p, k]], b[[p, k]]  # rijen verwisselen
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]
            A[i, k:] -= m * A[k, k:]
            b[i] -= m * b[k]
            bewerkingen += (n - k) + 2
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):                  # terug: van onder naar boven invullen
        x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
        bewerkingen += n - i
    return x, bewerkingen


x, ops = gauss([[3, -1, -2], [-1, 4, -1], [-2, -1, 4]], [10, 0, 0])
print("brug met ronde waarden (mA):", np.round(x, 4), " maal 7:", np.round(7 * x, 6))

R1, R2, R3, R4, R5 = 1.0, 2.2, 1.5, 2.2, 1.0          # de brug van hoofdstuk 21 van Intermediate (kΩ)
M = [[R1 + R2, -R1, -R2], [-R1, R1 + R3 + R5, -R5], [-R2, -R5, R2 + R4 + R5]]
I, _ = gauss(M, [10, 0, 0])
print(f"brug van hoofdstuk 21: stroom door R5 = {I[2] - I[1]:.4f} mA")

rng = np.random.default_rng(5)                      # hoe groeit het werk met het aantal onbekenden?
for n in (10, 20, 40, 80, 160):
    A, b = rng.normal(size=(n, n)), rng.normal(size=n)
    x, ops = gauss(A, b)
    verschil = np.max(np.abs(x - np.linalg.solve(A, b)))
    print(f"n = {n:3d}: {ops:8d} bewerkingen = {ops / n**3:.3f} n³, verschil met numpy {verschil:.0e}")