The Museum of HP Calculators

Polynomial Curve Fitting

for HP-41CX with W&W CCD Module

This program is supplied without representation or warranty of any kind. Hans Brueggemann and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Introduction

After finally getting my hands on the famous W&W CCD Module in 1985, I was wondering what could probably done with it in order to remove the 'nice-to-have' status from my HP-41 equipment. During that time I did lots of programming on the wonderful HP9845A, equipped with *all* peripherals you could dream of. We used it to calculate statistical data from experiments with power inverters and machines. Unfortunately, the HP9845A was too bulky to be moved around in the office or in the laboratory. In order to have the possibility to quickly analyze data on location, I decided to port a part of the HP-9845A Gen. Util. Routines Library to my HP41-CX, and put the CCD Module to some good work.

Overview

Please note:
In order to use the programs described herein, an extra (and not supplied in this article) subroutine "INV", capable of calculating the inverse of a quadratic matrix, is necessary. A listing of such a subroutine can be found in the W&W CCD manual.

What You will find here is a set of three programs that enables You to make a polynomial curve fit of a given data set x(*) and y(*). After supplying the number N of data sets , the x and y data, and the desired degree D of the polynomial , the coefficients C(*) of the polynomial are calculated. The resulting polynomial is of the form
f(x) = a0 + a1x + a2x^2 + a3x^3 + .. + aDx^D, where ak = C(k+1).

Once the coefficients are calculated, the polynomial can be evaluated at any x to check the fit.

Required Hardware

HP-41CX (or HP-41CV with XF/M)
W&W CCD A (or B)
HP 82143A or,
(alternative 1) HP 82162A with HP 82160A
(alternative 2) HP 82242A with HP 82240 ( change pgm lines 53, 94 in "POLFITD" to ADV respectively, because the PRL routine of the W&W CCD is not compatible with this printer- ROM.

Usage

In order to run the polynomial fit, "POLFITD", "POLFIT", "POLYN", and "INV" (a subroutine for matrix inversion, provided with the manual of the CCD ROM) must be present in calculator memory. Set SIZE to 031.
XEQ "POLFITD" and key in the following data: the degree of the polynomial, the number of data sets (x,y) you can provide, all x data, and all y data. Please note: this program is a rather crude implementation of data input routines. It neither checks the input nor does it provide any means of input correction. If it asks you for more data (x or y) than you have specified, just press R/S. the purpose of "POLFITD" is to gather all the data and call "POLFIT". After "POLFIT" has finished the calculation, "POLFITD" will automatically print the results.
Now key in any desired x - value and XEQ "POLYN" to calculate y = f(x). "POLYN" expects the coefficients vector C in XF/M (as delivered by "POLFIT").

General Comments

*always* key in numbers with maximum obtainable precision. on 'lazy' rounding of input data ( i.e. 3.142 instead of 3.141592654), "POLFIT" will almost certainly come up with a totally unreliable table of coefficients.
After starting "POLFITD" with degree > 4, go get a coffee.
All input data will be kept after program has finished. See register table of "POLFIT" for storage location of particular data.
For another complete run with same data, GTO. 041( "in POLFITD"); R/S
For another printout, XEQ "COUT"
Testing the program by using the data set from the HP-9845A Gen. Util. Routines Library 'Family Regression', (can be found elsewhere on the MoHP CDs) gives the following results (calc time approx. 5 min): a0: 0,84870195340; a1: 2,04504820000; a2: -0,33243480460; a3: 0,01832369164
The names of printed results are consistent with those found in the 'Family Regression' library.
You will find speckles of 'Synthetic Programming' anywhere in the code. See the comments on the respective lines, they will give You hints on workarounds.
Operations that involve the CCD- Module are commented, so there should be a chance to rewriting the code for those who do not own the CCD- Module. (surely not for the faint of the heart). Note, that CCD matrix operations generally involve the alpha- registers in order to adress the proper datafiles in XF/M.
Improve the code and inform the author.

Free Software

POLFITD, POLFIT and POLYN, Polynomial Curve Fitting for the HP-41C series Calculator
Copyright (c) 2000, Hans Brueggemann

These programs are free software; you can redistribute them and/or modify them under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 or any later version. These programs are distributed in the hope that they will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. This license is available from: http://www.gnu.org/copyleft/gpl.html

You can contact the author at: porkpiehat@gmx.net

Program Listing of "POLFITD"

 01 LBL "POLFITD"	(english translations where necessary)
 02 .03	
 03 CLRGX	
 04 "GRAD?"		"degree?"	
 05 PROMPT	
 06 STO 26	
 07 "ANZ. DATEN?"	"# of data sets?"	
 08 PROMPT	
 09 STO 25	
 10 "XA"	
 11 MDIM		DIM XA(N)	
 12 "Y"	
 13 MDIM		DIM Y(N)	
 14 LBL "IN"	
 15 "XA"	
 16 CLX			make XA the active vector	
 17 IJ=A		and seekpt to first element	
 18 LBL 40	
 19 ?IJ			get pointer of XA	
 20 "X"	
 21 ARCLI		concat "X" + pointer	
 22 "&127:"		"<app>:"	
 23 PROMPT		ask for Xn	
 24 SF 25	
 25 >R+			put Xn into vector & inc pointer	
 26 FS? 25	
 27 GTO 40	
 28 "Y"	
 29 CLX			make Y the active vector	
 30 IJ=A		and seekpt to first element	
 31 LBL 50	
 32 "Y"	
 33 ?IJ	
 34 ARCLI	
 35 "&127:"		"<app>:"	
 36 PROMPT		ask for Yn	
 37 SF 25	
 38 >R+			put Yn into vector & inc pointer	
 39 FS? 25	
 40 GTO 50	
 41 TIME		get the actual time just for fun	
 42 STO 30	
 43 XEQ "POLFIT"	do the trick...	
 44 ENG 9	
 45 F/E			FIX/ENG autoformat of display	
 46 TIME	
 47 RCL 30	
 48 HMS-		how long did it take ?	
 49 PRX			amaze user	
 50 LBL "COUT"		(the output routine)	
 51 "KOEFFIZIENTEN:"	"coefficients:"	
 52 PRA	
 53 PRL			PRint Line of "-" (PRL won't work with 82242A)	
 54 "C"	
 55 CLX			make C the active vector	
 56 IJ=A		and seekpt to first element	
 57 RCL 26		
 58 1		
 59 +		
 60 1 E3		(you could save a byte here and there)	
 61 /		
 62 1		
 63 +		
 64 STO 09		set up coeff indexer	
 65 ENG 9		show all digits in results	
 66 LBL 88		
 67 "a"		
 68 RCL 09		
 69 1		
 70 -		
 71 ARCLI		ARCL integer (stackx)	
 72 "&127: "		"<app>: "	
 73 R>+			get coeff to stackx & inc pointer	
 74 ARCL X		
 75 PRA			print coefficient	
 76 ISG 09		next index	
 77 GTO 88		
 78 ADV		
 79 "F: "		
 80 ARCL 16		
 81 PRA		
 82 "REGRESSION:"		
 83 17		
 84 XEQ 01		
 85 "RESIDUUM:"		
 86 18		
 87 XEQ 01		
 88 "TOTAL:"		
 89 SF 01		
 90 19		
 91 LBL 01		general formatting routine	
 92 ADV		
 93 PRA			print caption	
 94 PRL			PRint Line (PRL won't work with 82242A!!)
 95 "Df: "		
 96 ARCL IND X		
 97 PRA		
 98 6		
 99 -		
100 "SS: "		
101 ARCL IND X		
102 PRA		
103 FS?C 01		
104 RTN		
105 3		
106 +		
107 "MS: "		
108 ARCL IND X		
109 PRA		
110 END			(302 Bytes)

Program Listing of "POLFIT"

General Settings:

SIZE:	min. 027
Flags: none
Registers:
stack, alpha
 0	FNRM (M)
 1
 2
 3	I
 4	J
 5	K
 6	L
 7	INT(K)
 8	D+1
 9	Ia
10	temporary Sum
11	Reg SS
12	Res SS
13	Total SS
14	Reg MS
15	Res MS
16	F
17	DF Reg
18	DF Res
19	DF Total
20	X1 = SUMx/N
21	X2 = SUMx^2
22	Y1 = SUMy/N
23	Y2 = SUMy^2
24	Z = SUMxy
25	N
26	D

 01 LBL "POLFIT"	
 02 RCL 25	
 03 2	
 04 -	
 05 RCL 26		degree D	
 06 X>Y?	
 07 RTN			oops, not enough data sets...	
 08 E1			10
 09 ENTER^	
 10 ENTER^	
 11 RCL 25		# of data	
 12 "C"	
 13 MDIM		DIM C(N)	
 14 "Y,C"	
 15 MOVE		move vector Y to C	
 16 E-3			0.001	
 17 RCL 26	
 18 E1	
 19 +	
 20 STO 08		D+1	
 21 "B"	
 22 MDIM		DIM B(D+1)	
 23 *	
 24 LASTX	
 25 +	
 26 "M"	
 27 MDIM		reDIM M	
 28 RCL 25	
 29 RCL 26	
 30 STO 17	
 31 E	1	
 32 "B"	
 33 IJ=A		make B the active vector	
 34 +	
 35 -	
 36 STO 18	
 37 RCL 17	
 38 +	
 39 STO 19	
 40 LASTX	
 41 E3			1000	
 42 /	
 43 STO 05		for K = 0 to D	
 44 RCL 25	
 45 E3	
 46 /	
 47 ISG X	
 48 STO 09		prepare loop counter Ia
 49 LBL 66	
 50 RCL 05		get K	
 51 INT	
 52 STO 07	
 53 RCL 26	
 54 E3	
 55 /	
 56 +	
 57 STO 04		for J = K to D	
 58 LBL 07	
 59 RCL 09	
 60 STO 03		for I = 1 to N	
 61 RCL 07	
 62 RCL 04	
 63 +			K + J	
 64 INT	
 65 .			0.0	
 66 "XA"	
 67 IJ=A		XA(1)	
 68 LBL 08	
 69 R>+			get element & inc pointer	
 70 RCL Z	
 71 Y^X	
 72 +			stackx = SUM(I=1,N; XA(I)^(K+J)) 
 73 ISG 03	
 74 GTO 08		next I	
 75 STO 10	
 76 "M"	
 77 E-3	
 78 RCL 07	
 79 E	
 80 +			stackx = K,J	
 81 ST* Y	
 82 RCL 04		stacky = J,K	
 83 LASTX	
 84 +	
 85 INT	
 86 ST+ Z	
 87 E3	
 88 /	
 89 +	
 90 IJ=A	
 91 RCL 10		M(K,J) =	
 92 >R+	
 93 RCL Z		stackx =	
 94 IJ=	
 95 X<>Y	
 96 >R+			M(J,K)	
 97 ISG 04	
 98 GTO 07		next J	
 99 "C,XA,C"	
100 RCL 07	
101 X#0?		x not equal 0 ?	
102 M*			C*XA --> C	
103 SUM			sum up all elements of C	
104 "B"	
105 ?IJA	
106 X<>Y	
107 >R+			B(K) = SUM(I=1,N; Y(I)*XA(I)^K)	
108 ISG 05	
109 GTO 66		next K	
110 RCL 08	
111 "C"	
112 MDIM		reDIM coefficient vector C to C(D+1)	
113 "M,B,C"	
114 XEQ "INV"		invert M to M^(-1)	
115 M*M			M^(-1)* B --> C
116 RCL 25	
117 MDIM	
118 "B"	
119 PURFL		discard B	
120 "XA"	
121 SUM			sum up all elements of XA	
122 X<>Y	
123 /	
124 STO 20		X1 = (SUM(all;x))/N	
125 LASTX	
126 FNRM		Frobenius Norm of XA	
127 X^2	
128 STO 21		X2 = SUM(all;x^2)	
129 "Y"	
130 SUM	
131 RCL Z	
132 /	
133 STO 22		Y1 = (SUM(all;y))/N	
134 X^2	
135 RCL 25	
136 *	
137 CHS	
138 FNRM	
139 X^2	
140 STO 23		Y2 = SUM(all;y^2)
141 +	
142 STO 13	
143 "XA,Y,M"	
144 M*			M = XA * Y	
145 "M"	
146 SUM			sum up all elements of M	
147 STO 24		Z = SUM(all; x*y)
148 PURFL		ciao M!	
149 LBL "R"	
150 "XA"	
151 CLX	
152 IJ=A		stackx = XA(1)	
153 STO 11	
154 LBL 14	
155 "XA"	
156 ?IJA	
157 R>+			stackx = XA(Ia)	
158 RCL 08	
159 "C"	
160 IJ=A		select C(D)	
161 2	
162 -			loop counter	
163 STO 06		for Horner scheme	
164 X<>Y	
165 ENTER^	
166 ENTER^	
167 ENTER^		stack = XA(Ia)	
168 R>-			get element & dec pointer	
169 *	
170 LBL 15	
171 R>-			get element & dec pointer	
172 +			do Horner scheme	
173 *			for 
174 DSE 06		J = SUM(L=0,D; XA(Ia)^L*C(L))
175 GTO 15	
176 R>-	
177 +	
178 RCL 22	
179 -	
180 X^2	
181 ST+ 11		SUM(all;((J-(SUM(all;y))/N))^2)
182 ISG 09	
183 GTO 14		next Ia	
184 RCL 13	
185 RCL 11	
186 -	
187 STO 12	
188 LASTX	
189 RCL 17	
190 /	
191 STO 14	
192 X<>Y	
193 RCL 18	
194 /	
195 STO 15	
196 /	
197 STO 16	
198 END			(341 Bytes)

Program Listing of "POLYN"

01 LBL "POLYN"		expects x in stackx, returns y = f(x) in stackx
02 "C"	
03 DIM			get polynomial degree D	
04 INT	
05 IJ=A			make C the active matrix	
06 2	
07 -	
08 STO [		STO M	
09 X<>Y	
10 ENTER^	
11 ENTER^	
12 ENTER^		prepare for Horner scheme	
13 R>-	
14 *	
15 LBL 01	
16 R>-	
17 +			stackx = SUM(K=2,D;C(K)*x^K)	
18 *	
19 DSE [		DSE M	
20 GTO 01	
21 R>-	
22 +			add C1 (i.e., a0)	
23 END			(42 Bytes)

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall