From: dfk@hp-lsd.cos.hp.com (David F. Kurth) Newsgroups: comp.sources.hp48 Subject: v01i011: polyfit10 - Polyfit v1.0: produce polynomial given N data points, Part01/01 Date: 12 Aug 91 01:49:19 GMT Followup-To: comp.sys.hp48 Checksum: 3154341908 (verify with brik -cv) Submitted-by: David F. Kurth Posting-number: Volume 1, Issue 11 Archive-name: polyfit10/part01 BEGIN_DOC polyfit.doc I've always thought there should be a good use for a program like this, but I still haven't found it. It was a good challenge to program the algorithm, and maybe someone else can find a good use for it. POLYFIT produces the coefficients of an N-1 degree polynomial given N data points that lie on the curve of that polynomial. For example, consider the following table of data points: X Y ----------- -2 1 -1 -1 0 1 1 -1 2 1 There are 5 data points, and the Y value a very symmetric. POLYFIT produces the following polynomial: f(x) = 2/3*x^4 + 0*x^3 - 8/3*x^2 + 0*x^1 + 1*x^0 If you plot this with an X range between -2 and 2, you will see the curve "oscillate" between +1 and -1 as it hits all five data points. POLYFIT is not a least square's type curve fit program. There is no minimization of errors. The polynomial coefficients found by POLYFIT will exactly hit the given data points (to within the accuracy of the calculator). Using POLYFIT Enter the following code into your calculator. After loading, enter the new directory if not already there. Push the CST key load the custom menu which should reveal the following keys: POLY - Executes the main program which takes the input array of X-Y points from the stack and produces the K array of polymonial coefficients. VEIW - Displays the coefficients of the resulting polynomial. FX - Calculates f(x) using level 1 of the stack for x. X-Y - Recalls the input data points to the stack for editing or viewing. To begin, you must first enter the input data points array. This is easiest to do by using the Matrix Writer (Right-shift ENTER). Enter your first X-Y pair of points, then push the down arrow which informs the calculator that the array is only 2 elements wide. Continue entering the rest of your points until done. Then just push ENTER to return the array to the stack. Now push POLY to calculate the coefficients. The display should produce "DOING L:" which can be ignored (the program calculates a triangular array of difference values used for the rest of the computations). After this, the Ki coefficient values are calculated and briefly displayed. When done, the K values can be displayed again by using the VIEW key. Output values can be calculated by entering an X value on the stack and pressing the FX key. Plotting the equation can be done by pressing Left-shift PLOT, PLOTR, and then AUTO. For a small number of data points, the calculations are fairly quick. For each additional data point entered, calculation time roughly doubles. I've calculated 15 data points, but it took the calculator about 7.5 hours to finish. However, with sufficient memory and batteries, have at it. Dave Kurth dfk@hp-lsd.cos.hp.com END_DOC BYTES: #B610h 2267.5 BEGIN_RPL polyfit %%HP: T(3)A(D)F(.); 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 \>> END_RPL