'program lin_pol6.bas
'pol6 least square
' y=a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5  +a6*x^6
 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 = "
a$(7)="a6 = "
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
    a56=a56+x^(6)*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)
    a14=a14+x^(11)
    a15=a15+x^(12)
    next i

    M=7
 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)=a9:   C(1,8)=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)=a10:  C(2,8)=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)=a11:  C(3,8)=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)=a12:  C(4,8)=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)=a13:  C(5,8)=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)=a14:  C(6,8)=a55
 C(7,1)=a9:C(7,2)=a10:C(7,3)=a11:C(7,4)=a12:C(7,5)=a13:C(7,6)=a14:C(7,7)=a15:  C(7,8)=a56
    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