'program lin_pol5.bas
'pol5 least square
' y=a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5
 dim C(50,50): dim x(40): dim y(40)

a$(1)="a0 = "
a$(2)="a1 = "
a$(3)="a2 = "
a$(4)="a3 = "
a$(5)="a4 = "
a$(6)="a5 = "

x(1)=0 :x(2)=1 :x(3)=1.5:x(4)=2.3:x(5)=2.5:x(6)=4
y(1)=.2:y(2)=.8:y(3)=2.5:y(4)=2.5:y(5)=3.5:y(6)=4.3

x(7)=5.1 :x(8)=6 :x(9)=6.5:x(10)=7  :x(11)=8.1 :x(12)=9
y(7)=3.0 :y(8)=5 :y(9)=3.5:y(10)=2.4:y(11)=1.3 :y(12)=2

x(13)=9.3 :x(14)=11   :x(15)=11.3:x(16)=12.1  :x(17)=13.1 :x(18)=14
y(13)=-.3 :y(14)=-1.3 :y(15)=-3.0:y(16)=-4.0  :y(17)=-4.9 :y(18)=-4

x(19)=15.5:x(20)=16.0 :x(21)=17.5:x(22)=17.8:x(23)=19.0 :x(24)=20.0
y(19)=-5.2:y(20)=-3.0 :y(21)=-3.5:y(22)=-1.6:y(23)=-1.4 :y(24)=-0.1

     a=0:a1=0:a2=0:a3=0:a4=0:a11=0:a12=0:a54=0:a50=0:a51=0
    a5=0:a6=0:a7=0:a8=0:a9=0:a10=0:a13=0:a55=0:a52=0:a53=0
    n=24
    for i=1 to n
    x=x(i):y=y(i)
    a50=a50+y
    a51=a51+x*y
    a52=a52+x^(2)*y
    a53=a53+x^(3)*y
    a54=a54+x^(4)*y
    a55=a55+x^(5)*y

    a4=a4+x^(1)
    a5=a5+x^(2)
    a6=a6+x^(3)
    a7=a7+x^(4)
    a8=a8+x^(5)
    a9=a9+x^(6)
    a10=a10+x^(7)
    a11=a11+x^(8)
    a12=a12+x^(9)
    a13=a13+x^(10)
    next i

    M=6
 C(1,1)=n :C(1,2)=a4:C(1,3)=a5 :C(1,4)=a6 :C(1,5)=a7 :C(1,6)=a8:  C(1,7)=a50
 C(2,1)=a4:C(2,2)=a5:C(2,3)=a6 :C(2,4)=a7 :C(2,5)=a8 :C(2,6)=a9:  C(2,7)=a51
 C(3,1)=a5:C(3,2)=a6:C(3,3)=a7 :C(3,4)=a8 :C(3,5)=a9 :C(3,6)=a10: C(3,7)=a52
 C(4,1)=a6:C(4,2)=a7:C(4,3)=a8 :C(4,4)=a9 :C(4,5)=a10:C(4,6)=a11 :C(4,7)=a53
 C(5,1)=a7:C(5,2)=a8:C(5,3)=a9 :C(5,4)=a10:C(5,5)=a11:C(5,6)=a12 :C(5,7)=a54
 C(6,1)=a8:C(6,2)=a9:C(6,3)=a10:C(6,4)=a11:C(6,5)=a12:C(6,6)=a13 :C(6,7)=a55

     N=M+1
     FOR J=1 TO M
         MYI=J
  [LOOPMYI]
         I=MYI
         IF C(I,J)<>0 GOTO [CALC]
         MYI=MYI+1
         IF MYI<=M THEN [LOOPMYI]

        PRINT "SORRY, NO UNIQUE SOLUTION"
        END

  [CALC]
      FOR K=1 TO N
         X=C(J,K)
         C(J,K)=C(I,K)
         C(I,K)=X
      NEXT K
     Y=1/C(J,J)
     FOR K=1 TO N
       C(J,K)=Y*C(J,K)
     NEXT K

     FOR I=1 TO M
       IF I=J GOTO [LEAP]
       Y=C(I,J)*(-1)
         FOR K=1 TO N
           C(I,K)=C(I,K)+Y*C(J,K)
         NEXT K
      [LEAP]
     NEXT I
  NEXT J

     FOR J=1 TO M
   PRINT a$(J);C(J,N)'solution vector
   REM
 NEXT J
 END