how do you solve a system of linear equations in ThinBasic?

[jack]

hi Bob and Dan
Bob, Octave does much beter than thinBasic, Dan, I did change the script to make it work.
anyway, you can also drop the first column in the matrix and just make the first coefficient = 1,
(also decrease the size of the dimensions by one) the rest of the coefficients will differ
slightly because of this change, but it approximates the gamma function almost as well, (it needs to be checked).

example thinBasic script, look carefully for the changes.

Uses "console"
Uses "math"
Const n = 6 
Dim As Ext g = 5 
Dim As Ext a(n,n), b(n,n), c(n,1) 
Dim As Integer i, j 
'for n = 6, the following sets up matrix a as: 
' [[ 1,   1/2, 1/3, 1/4,  1/5,  1/6 ] 
' [  1/2, 1/3, 1/4, 1/5,  1/6,  1/7 ] 
' [  1/3, 1/4, 1/5, 1/6,  1/7,  1/8 ] 
' [  1/4, 1/5, 1/6, 1/7,  1/8,  1/9 ] 
' [  1/5, 1/6, 1/7, 1/8,  1/9,  1/10 ] 
' [  1/6, 1/7, 1/8, 1/9,  1/10, 1/11 ]] 
 
For i=1 To n
                       '<< note the change 
  For j=1 To n  '<< note the change
    a(i,j)=1/(i+j-1)  '<< note the change
  Next
Next
For i=1 To n 
  c(i,1) = Factorial(i-1)/((i+g-0.5)^(i-0.5)/Exp(i+g-0.5)*Sqr(2*M_PI))-1 '<< note the change
Next
MAT b() = INV(a()) 
MAT a() = b() * c()
Console_WriteLine "C( 1) = 1" 
For i=1 To n 
  Console_WriteLine "C("&Str$(i+1)&") = "&a(i,1) 
Next 
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey

here are the coefficients as calculated with Maple

c[0] := 1
c[1] := 76.180091728331374539146745
c[2] := -86.505320289513654614503042
c[3] := 24.014097921606006186204527
c[4] := -1.2317386147754424755938420
c[5] := .120745388047503859955038e-2
c[6] := -.4868518292851569517614e-5

[danbaron]

It seems that the new values disagree more with the values from Numerical Recipes.

' code -----------------------------------------------------------------------------

Uses "console"
Uses "math"
Const n = 6
Dim As Ext a(n), b(n), c(n)
Dim As Integer i

' From Johan.

a(1) = 76.1800917283313588
a(2) = -86.5053202895132074
a(3) = 24.0140979216034793
a(4) = -1.23173861476925595
a(5) = 1.2074538708475302 * 10^-3
a(6) = -4.8685153615224408 * 10^-6

' From Numerical Recipes.

b(1) =  76.18009172947146
b(2) = -86.50532032941677
b(3) =  24.01409824083091
b(4) =  -1.231739572450155
b(5) =   0.1208650973866179 * 10^-2
b(6) =  -0.5395239384953 * 10^-5

' abs(difference).

For i = 1 To n
c(i) = Abs(a(i)-b(i))
Next

For i = 1 To n
PrintL c(i)
Next

Console_WaitKey

' output ---------------------------------------------------------------------------

1.14010120094709677E-9
3.99035626000299182E-8
3.19227430699633996E-7
9.57680899050031287E-7
1.19710301864880011E-6
5.267240234305592E-7