%%HP: T(3)A(D)F(.);
@ POLYFIT by David F. Kurth
DIR
POLY
@ This program takes an array of N X-Y values from the stack
@ and calculates the coefficients of an N-1 polynomial that
@ intersect with all N input points.  The N X-Y values are most
@ easily entered using the Matrix Writer Application.

@ Directory
@ Checksum # B610h
@ Bytes    2267.5

\<< IF DEPTH                            @ Check for an object on stack
  THEN
    IF DUP TYPE 3 ==                    @ There is object, check if array
    THEN
      DUP 'XY' STO SIZE 1 GET 'N' STO   @ Save Array and Number of points
      CLLCD                             @ Clear screen for progess updates
      LijCAL                            @ Calculate L coefficients
      KCAL                              @ Calculate final K coefficients
    ELSE EMESS                          @ Display error message
    END
    {CST POLY VIEWK Fx}                 @ Restore main variable order
    ORDER                               @ after new variables are created
  ELSE
    EMESS                               @ Display error message
  END
\>>
CST
{                                       @ Custom Menu for Polyfit program
  {"POLY" POLY}                         @ Main program to find coefficients
  {"VIEW" VIEWK}                        @ View polynomial coefficients
  {"Fx" Fx}                             @ Given X (on stack) return f(X)
  {}                                    @ Nothing
  {}                                    @ Nothing
  { "X-Y" XY}                           @ Recall input array to stack for editing
}
Cab
\<< @ Generate terms to sum for C(a)(b) @
  IF DUP2 ==
  THEN DROP2 1                          @ If a=b, the Cab=1
  ELSE
    0 1 \-> a b SUM Index
    \<<
      { 1 }                             @ Initial list

      DO  @ For all possible terms @
        WHILE @ Build list up to a-b terms @
          Index a b - <                 @ Until we hit last term
        REPEAT
          DUP Index GET 1 + 1 \->LIST +  @ Add next term to list
          Index 1 + 'Index' STO
        END

        DO  @ Increment last terms until limit @
          FETCH                         @ Get X values and multiply
          SUM + 'SUM' STO               @ Add to sum and save
          DUP Index GET 1 + Index SWAP
          DUP 4 ROLLD PUT               @ Increment last term and save
        UNTIL SWAP b Index + >          @ Until we go over b+Index limit
        END

        IF Index 1 >  @ Increment previous terms @
        THEN                            @ If possible
          DO
            Index 1 - 'Index' STO       @ Decrement index
            1 Index SUB                 @ Shorten list
            DUP Index GET 1 + Index SWAP
            DUP 4 ROLLD PUT             @ Increment this term and save
          UNTIL
            SWAP b Index + \<=          @ This term at limit
            Index 1 == OR               @ or this is last term
          END
        END
      UNTIL
        DUP Index GET b Index + >
        Index 1 == AND
      END
      DROP                              @ Remove last list from stack
      SUM
      -1 a b - ^ *                      @ Fix sign of product
    \>>
  END                                   @ If a=b
\>>
EMESS
\<< @ Display error message if no stack value, or value is not an array @
  CLLCD
  "Input array not found." 1 DISP
  "Use Matrix Writer to " 3 DISP
  "enter your X Y values" 4 DISP
  "into an array on the" 5 DISP
  "stack.  Then start" 6 DISP
  "POLY again." 7 DISP
  7 FREEZE
\>>
EQ
Y.EQ                                    @ For PLOTR to plot f(x)
FETCH
\<< @ Use Index list on stack to fetch X values and multiply @
  DUP SIZE 1 \-> lst n product          @ Save term list, size, and product
  \<<
    XY
    1 n
    FOR i
      DUP lst i GET 1 XYindex GET       @ Get Xi term
      product * 'product' STO           @ Multiply and save result
    NEXT
    DROP                                @ Drop XY array
    lst product                         @ Restore term list and product
  \>>
\>>
Fx
\<<  @ Calculate F(x) using K coefficients and x value from stack @
  K N GET                               @ First coefficient
  N 1 - 1
  FOR i                                 @ For i = N-1 to 1
    OVER *                              @ x times
    K i GET +                           @ Add in next coefficient
    -1
  STEP
  SWAP DROP                             @ Remove x value
\>>
KCAL
\<< @ Calculate K coefficients (f(x) = Kn*x^n-1 + ... K2*x + K1) @
  1 N
  START                                 @ Create K array with N zeroes
    0
  NEXT
  N \->ARRY 'K' STO                     @ Save K array

  1 N
  FOR i                                 @ i=1 to N
    "DOING K" STD i \->STR + 1 DISP     @ Display Ki value
    0                                   @ Zero sum
    i 1 - N 1 -
    FOR j                               @ j=i-1 to N-1
      @ Ki = Ki + Lj1 * C(j)(i-1) @
      L j 1 Lindex GET                  @ Get L term from array
      j i 1 - Cab                       @ Get Cab coefficient
      * +                               @ Sum to Ki
    NEXT
    K i 3 ROLL PUT 'K' STO              @ Add K value to array
    "K" i \->STR + "=" + K i GET \->STR + 3 DISP
  NEXT
\>>
LijCAL
@ Calculate triangular table of L values @
\<<
  "DOING L:" 1 DISP
  1 TMAX                                @ Create L values array
  START                                 @ With TMAX zeroes
    0
  NEXT
  TMAX \->ARRY                          @ L array initialized
  1 N                                   @ Set Lj,1 = Yj
  FOR j                                 @ For i=0, j = 1 to N
                                        @ L array already on stack
    j                                   @ L index value to place Y value
    XY j 2 * GET                        @ Y value for L array
    PUT                                 @ Y value now in L array
  NEXT
                                        @ Initialized L array left on stack
  1 N 1 -
  FOR i                                 @ For i = 1 to N-1
    1 N i -
    FOR j                               @ For j = 1 to N-i
      DUP DUP i 1 - j 1 + Lindex GET    @ Get L i-1,j+1 value
      SWAP i 1 - j Lindex GET -         @ Get L i-1,j value and subtract
      XY DUP i j + 1 XYindex GET        @ Get X i+j value
      SWAP j 1 XYindex GET -            @ Get X j value and subtract
      / i j Lindex SWAP PUT             @ Save new L i,j into array
    NEXT
  NEXT
  'L' STO                               @ Save final array
\>>
Lindex
@ Return index into L array of i,j element @
\<<
  \-> i j                               @ Save index values
  \<<
    i 1 - N i 2 / - * N + j +           @ index=(i-1)(N-i/2)+j+N
  \>>
\>>
PPAR
{
  (-6.5,-3.1)                           @ X range
  (6.5,3.2)                             @ Y range
  X                                     @ Independent Variable
  0                                     @ Resolution (pixel in each column)
  (0,0)                                 @ Axes intersection
  FUNCTION                              @ Function type
  Y                                     @ Dependent Variable
}
TMAX
@ Maximum Number of elements in the L triangular array @
\<<
  N DUP 1 + * 2 /                       @ TMAX = N(N+1)/2
\>>
VIEWK
\<< @ Display the K coefficients of the resulting polynomial expression @
  STD CLLCD "Hit VIEW to cont..." 1 DISP
  N 1
  FOR i                                 @ For each coefficient
    "K" i \->STR + "=" +
    K i GET \->STR + "*x^" + i 1 -
    \->STR + 3 DISP 0 WAIT DROP
    -1
  STEP
  CLLCD                                 @ Clear display to let user know
                                        @ the program is still working
  KILL                                  @ Get rid of HALT annunciator
\>>
XYindex
@ Return index into XY array of i,j element @
\<<
  \-> i j                               @ Save index values
  \<<
    i 1 - 2 * j +                       @ index=2(i-1)+j
  \>>
\>>
Y.EQ
\<< @ Fx function requires X value from stack, not compatible with PLOTR. @
  @ This program takes 'X' from PLOTR function, calls Fx, and @
  @ returns with f(x) value to satisfy PLOTR requirements. @
  X Fx
\>>
