The Museum of HP Calculators
The Museum of HP Calculators
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.
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.
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.
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
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)
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)
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