The Museum of HP Calculators

# for HP-41CX with W&W CCD Module

This program is © 2000 by Hans Brueggemann and is used here by permission.

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

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

1. HP-41CX (or HP-41CV with XF/M)
2. W&W CCD A (or B)
3. HP 82143A or,
4. (alternative 1) HP 82162A with HP 82160A
5. (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

1. 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.
2. 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.
3. 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").

1. *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.
2. After starting "POLFITD" with degree > 4, go get a coffee.
3. All input data will be kept after program has finished. See register table of "POLFIT" for storage location of particular data.
4. For another complete run with same data, GTO. 041( "in POLFITD"); R/S
5. For another printout, XEQ "COUT"
6. 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
7. The names of printed results are consistent with those found in the 'Family Regression' library.
8. 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.
9. 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.
10. Improve the code and inform the author.

## Free Software

POLFITD, POLFIT and POLYN, Polynomial Curve Fitting for the HP-41C series Calculator

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
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>:"
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>:"
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
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
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"
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)
```