//
// Gaussian Elimination Estimate
//
// FutureBasic 7.0.34 August 2025 R.W
//
//--------------------------------------------------
//
// Test data array. I wanted it to be base 1 and not zero
// hence the 99 which is just filler for the zeroeth element
//
double a(6,6) =  ¬
{99,  99,   99,   99,    99,    99,   99,¬
99, 1.00, 0.00, 0.00,  0.00,  0.00, 0.00,¬
99, 1.00, 0.63, 0.39,  0.25,  0.16, 0.10,¬
99, 1.00, 1.26, 1.58,  1.98,  2.49, 3.13,¬
99, 1.00, 1.88, 3.55,  6.70, 12.62, 23.80,¬
99, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28,¬
99, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02 }

double b(6) =  ¬
{99,-0.01, 0.61, 0.91, 0.99, 0.60, 0.02 }

double x(6)
double ab(6)
double ans(6),r(6),bb(6,6), f
int c, d, g, p1, p2, z
//

int row, col
int j, k, n = 6
//
// Pivot n times
local fn pivot(num as int)
  For p1  = num to n - 1
    For p2 = p1 + 1 To n
      If abs(bb(p1,num)) < abs(bb(p2,num))
        Swap r(p1),r(p2)
        For g = 1 To n
          Swap bb(p1,g),bb(p2,g)
        Next g
      End If
    Next p2
  Next p1
end fn

// The 'Gaussian' function takes two matrices,
// 'a' and 'b', with 'a' square, and it returns
// the determinant of 'a' and a matrix 'x'
// such that a * x = b.
// If 'b' is the identity, then 'x' = inverse of 'a'

local fn Gaussian(matrix(6,6) as double, rhs(6) as double)

  // copy matrix
  for c = 1 to n
    r(c) = rhs(c)
    for d = 1 to n
      bb(c,d) = matrix(c,d)
    next d
  next c

  // pivot it (hey that's what the task requires
  //            don't ask me why. It doesn't
  //            improve accuracy with these numbers)
  for k = 1 to n-1
    fn pivot(k)
    for row = k to n-1
      if bb(row + 1, k) == 0 then exit for
      f = bb(k, k) / bb(row+1, k)
      r(row + 1) = r(row + 1) * f - r(k)
      For g = 1 To n
        bb((row + 1), g) = bb((row + 1), g) * f - bb(k,g)
      Next g
    Next row
  Next k

  //back substitute
  For z = n to 1 step -1
    ans(z) = r(z) / bb(z,z)
    For j = n to z + 1 step -1
      ans(z) = ans(z) - (bb(z,j) * ans(j) / bb(z,z))
    Next j
  Next z
end fn

window 1,@"Gaussian Elimination"

print "TEST DATA USED"
print "                                        ax = transpose b"
print "[1.00 |0.00 |0.00 | 0.0000 | 0.00 |  0.00]   [??]   [-0.01]"
print "[1.00 |0.63 |0.39 | 0.2500 | 0.16 |  0.10]   [??]   [ 0.61]"
print "[1.00 |1.26 |1.58 | 1.9800 | 2.49 |  3.13] x [??] = [ 0.91]"
print "[1.00 |1.88 |3.55 | 6.7000 |12.62 | 23.80]   [??]   [ 0.99]"
print "[1.00 |2.51 |6.32 |15.8800 |39.90 |100.28]   [??]   [ 0.60]"
print "[1.00 |3.14 |9.87 |31.0100 |97.41 |306.02]   [??]   [ 0.02]"
print

fn Gaussian(a(0,0),b(0))

print "Solution:"
print "---------"
for j = 1 to 6
  print using "###.################";ans(j)
next j

HandleEvents