hp41programs

LeastSquaresPolynomials Least-Squares Polynomials for the HP-41

Overview

1°) Program#1
2°) Program#2
3°) Root-Mean-Square Error

-The 2 following routines only deal with equally-spaced abscissas.

-A similar program is listed in "Least-squares approximation for the HP-41"
-The variants hereunder are faster.

1°) Program#1

-Given a set of (N+1) data points ( x₀ , y₀ ) , ( x₁ , y₁ ) , ........................ , ( x_N , y_N )
-We define orthogonal polynomials P_k,N(t) by ( t = (x-x₀)/h with h = x_i+1 - x_i )

P_k,N(t) = SUM_i=0,1,...,k (-1)ⁱ C_kⁱ C_k+iⁱ t⁽ⁱ⁾/N⁽ⁱ⁾

where t⁽ⁱ⁾ = t(t-1)(t-2)....(t-i+1) and N⁽ⁱ⁾ = N(N-1)(N-2).....(N-i+1) ( = 1 if i = 0 )
and C_kⁱ = k!/(i!(k-i)!)

-The polynomial p(x) can then be written p(x) = a₀ P_0,N(t) + a₁ P_1,N(t) + .................... + a_m P_m,N(t)

with a_k = [ SUM_t=0,1,...,N y_t P_k,N(t) ]/[ SUM_t=0,1,...,N (P_k,N(t))² ]

-These coefficients are computed and stored when the program reaches line 149
-Then, the program calculates the coefficients in the more standard form: c₀ + c₁ t + c₂ t² + ...... + c_m t^m

-In fact, the P_k,N(t) are computed by the recurrence relations:

P_-1,N(t) = 0 , P_0,N(t) = 1 , (k+1)(N-k) P_k+1,N(t) = (2k+1)(N-2t) P_k,N(t) - k(N+k+1) P_k-1,N(t)

Data Registers: • R00 = N ( N+1 data points ) ( Registers R00 and R10 thru R10+N are to be initialized before executing "LSPN" )

R01 thru R09: temp

• R10 = y₀ • R11 = y₁ • R12 = y₂ ......................... • R10+N = y_N

>>>> at the end R11+N = a₀ , R12+N = a₁ , .............................. , R11+N+m = a_m
R12+N+m = c₀ , R13+N+m = c₁ , ........................... , R12+N+2m = c_m

R13+N+2m to R14+N+3m: temp

Flags: /
Subroutines: /

-Line 132 is a three-byte GTO

01 LBL "LSPN"
02 CLA
03 STO 01
04 RCL 00
05 11
06 +
07 STO M
08 +
09 .1
10 %
11 +
12 RCL 01
13 -
14 ENTER
15 INT
16 RCL 01
17 ST+ X
18 +
19 3
20 +
21 RCL 10
22 1
23 LBL 00
24 STO IND Z
25 X<>Y
26 STO IND T
27 X<>Y
28 ST+ Z
29 ISG T
30 GTO 00
31 RCL 00
32 E3
33 /
34 11.01
35 +
36 STO 03
37 SIGN
38 RCL M
39 +
40 RCL 01
41 +
42 0
43 STO IND Y
44 SIGN
45 ST+ Y
46 STO IND Y

47 STO 02
48 LBL 01
49 RCL M
50 RCL 01
51 +
52 ENTER
53 SIGN
54 +
55 STO 06
56 LASTX
57 +
58 STO 05
59 RCL 01
60 +
61 .1
62 %
63 +
64 RCL 01
65 -
66 STO 04
67 CLX
68 STO 07
69 ISG 04
70 LBL 02
71 RCL 07
72 ST+ X
73 RCL 00
74 RCL 02
75 ST+ X
76 -
77 ST* Y
78 +
79 RCL IND 05
80 *
81 RCL 00
82 RCL 07
83 ST+ Y
84 ST* Y
85 +
86 RCL IND 06
87 *
88 -
89 RCL 07
90 RCL 00
91 RCL 07
92 -

93 ST* Y
94 +
95 /
96 STO IND 04
97 CLX
98 SIGN
99 ST+ 05
100 ST+ 06
101 ST+ 07
102 ISG 04
103 GTO 02
104 RCL 01
105 +
106 ST- 04
107 RCL M
108 STO 05
109 RCL 01
110 ST+ X
111 +
112 PI
113 INT
114 +
115 STO 06
116 LBL 03
117 RCL IND 04
118 X^2
119 ST+ IND 06
120 LASTX
121 RCL IND 03
122 *
123 ST+ IND 05
124 CLX
125 SIGN
126 ST+ 05
127 ST+ 06
128 ISG 04
129 GTO 03
130 ST+ 02
131 ISG 03
132 GTO 01
133 RCL 03
134 INT
135 RCL 05
136 DSE X
137 E3
138 /

139 +
140 RCL 04
141 INT
142 LBL 04
143 RCL IND X
144 ST/ IND Z
145 CLX
146 SIGN
147 +
148 ISG Y
149 GTO 04
150 RCL 05
151 RCL 04
152 INT
153 E3
154 /
155 +
156 CLRGX
157 RCL 05
158 DSE X
159 E3
160 /
161 RCL 03
162 INT
163 +
164 ENTER
165 STO 03
166 CLX
167 LBL 05
168 RCL IND Y
169 +
170 ISG Y
171 GTO 05
172 STO IND 05
173 RCL 04
174 INT
175 RCL 06
176 DSE X
177 E3
178 /
179 +
180 STO 02
181 STO 08
182 ENTER
183 SIGN
184 ST+ Y

185 ST+ 03
186 ENTER
187 STO IND Z
188 ST+ Z
189 LBL 06
190 RCL Y
191 *
192 CHS
193 STO IND Z
194 ISG Y
195 CLX
196 ISG Z
197 GTO 06
198 RCL 01
199 1
200 STO N
201 +
202 .1
203 %
204 +
205 RCL 03
206 STO 09
207 +
208 STO 04
209 GTO 11
210 LBL 12
211 RCL 09
212 +
213 STO 09
214 INT
215 RCL M
216 -
217 STO 03
218 RCL 02
219 ENTER
220 SIGN
221 ST+ Y
222 LBL 10
223 ENTER
224 X<> IND Z
225 X<>Y
226 RCL 03
227 *
228 -
229 ISG 03
230 CLX

231 ISG Y
232 GTO 10
233 RCL 09
234 STO 03
235 LBL 07
236 RCL O
237 RCL 00
238 -
239 ST* N
240 SIGN
241 ST- O
242 LBL 11
243 RCL 03
244 INT
245 RCL M
246 -
247 STO 05
248 LBL 08
249 RCL 04
250 INT
251 ENTER
252 SIGN
253 -
254 RCL 01
255 -
256 RCL M
257 -
258 RCL 05
259 E3
260 /
261 +
262 STO 06
263 RCL N
264 STO 07
265 CLX
266 LBL 09
267 RCL 05
268 RCL 05
269 RCL 06
270 INT
271 ST+ Z
272 ST- Y
273 FACT
274 X^2
275 X<>Y
276 FACT

277 *
278 X<>Y
279 FACT
280 X<>Y
281 /
282 RCL 06
283 INT
284 ENTER
285 SIGN
286 ST+ 02
287 -
288 RCL 00
289 -
290 RCL 07
291 *
292 STO 07
293 /
294 RCL IND 02
295 *
296 +
297 ISG 06
298 GTO 09
299 RCL IND 03
300 *
301 ST+ IND 04
302 RCL 08
303 STO 02
304 SIGN
305 ST+ 05
306 ISG 03
307 GTO 08
308 ST+ 02
309 ST+ 08
310 ISG 04
311 GTO 12
312 RCL 03
313 FRC
314 RCL M
315 +
316 RCL 03
317 INT
318 RCL 04
319 FRC
320 +
321 CLA
322 END

( 439 bytes / SIZE N+3m+15 )

STACK	INPUTS	OUTPUTS
Y	/	bbb.eee(a_i)
X	m	bbb.eee(c_i)

Example:

t	0	1	2	3	4	5	6
y	1.2	3.4	5.1	6.7	7.6	8	8.4

N = 6 STO 00

1.2 STO 10         6.7 STO 13        8.4 STO 16
3.4 STO 11         7.6 STO 14
5.1 STO 12           8   STO 15

-If we seek a polynomial of degree 4

4 XEQ "LSPN" >>>> 22.026 = control number of the coefficients c_i's --- execution time = 1mn25s ---
RDN 17.021 = control number of the coefficients a_i's

and we get:

p(t) = 5.771429 - 3.567857 P_1,N(t) - 1.005952 P_2,N(t) - 0.016667 P_3,N(t) + 0.037013 P_4,N(t)
p(t) = 1.217965 + 2.088023 t + 0.093561 t² - 0.083586 t³ + 0.007197 t⁴

Notes:

-With m = N , youhave the Lagrange Polynomial.
-This program calculates the Stirling number of the 1st kind to get the c_i's from the a_i's

-With N = 87 and m = 13 execution time = 44m25s ( about 4s with V41 in turbo mode )

2°) Program#2

-It is actually faster to use the recurrence relation between P_k+1,N(t) P_k,N(t) & P_k-1,N(t) to compute the c_i's
-On the other hand, it uses more registers. This method is employed by the variant below.

Data Registers: • R00 = N ( N+1 data points ) ( Registers R00 and R10 thru R10+N are to be initialized before executing "LSPN" )

R01 thru R09: temp

• R10 = y₀ • R11 = y₁ • R12 = y₂ ......................... • R10+N = y_N

>>>> at the end R11+N = a₀ , R12+N = a₁ , .............................. , R11+N+m = a_m
R12+N+m = c₀ , R13+N+m = c₁ , ........................... , R12+N+2m = c_m

R13+N+2m to R14+N+4m: temp

Flags: /
Subroutines: /

-Line 132 is a three-byte GTO

01 LBL "LSPN"
02 STO 01
03 RCL 00
04 11
05 +
06 STO 08
07 +
08 .1
09 %
10 +
11 RCL 01
12 -
13 ENTER
14 INT
15 RCL 01
16 ST+ X
17 +
18 3
19 +
20 RCL 10
21 1
22 LBL 00
23 STO IND Z
24 X<>Y
25 STO IND T
26 X<>Y
27 ST+ Z
28 ISG T
29 GTO 00
30 RCL 00
31 E-3
32 STO 09
33 *
34 11.01
35 +
36 STO 03
37 SIGN

38 RCL 08
39 +
40 RCL 01
41 +
42 0
43 STO IND Y
44 SIGN
45 ST+ Y
46 STO IND Y
47 STO 02
48 LBL 01
49 RCL 08
50 RCL 01
51 +
52 ENTER
53 SIGN
54 +
55 STO 06
56 LASTX
57 +
58 STO 05
59 RCL 01
60 +
61 .1
62 %
63 +
64 RCL 01
65 -
66 STO 04
67 CLX
68 STO 07
69 ISG 04
70 LBL 02
71 RCL 07
72 ST+ X
73 RCL 00
74 RCL 02

75 ST+ X
76 -
77 ST* Y
78 +
79 RCL IND 05
80 *
81 RCL 00
82 RCL 07
83 ST+ Y
84 ST* Y
85 +
86 RCL IND 06
87 *
88 -
89 RCL 07
90 RCL 00
91 RCL 07
92 -
93 ST* Y
94 +
95 /
96 STO IND 04
97 CLX
98 SIGN
99 ST+ 05
100 ST+ 06
101 ST+ 07
102 ISG 04
103 GTO 02
104 RCL 01
105 +
106 ST- 04
107 RCL 08
108 STO 05
109 RCL 01
110 ST+ X
111 +

112 PI
113 INT
114 +
115 STO 06
116 LBL 03
117 RCL IND 04
118 X^2
119 ST+ IND 06
120 LASTX
121 RCL IND 03
122 *
123 ST+ IND 05
124 CLX
125 SIGN
126 ST+ 05
127 ST+ 06
128 ISG 04
129 GTO 03
130 ST+ 02
131 ISG 03
132 GTO 01
133 RCL 03
134 INT
135 RCL 05
136 DSE X
137 RCL 09
138 *
139 +
140 RCL 04
141 INT
142 LBL 04
143 RCL IND X
144 ST/ IND Z
145 CLX
146 SIGN
147 +
148 ISG Y

149 GTO 04
150 .1
151 %
152 +
153 STO 08
154 X<> 04
155 INT
156 STO 07
157 X<> 03
158 INT
159 X<>Y
160 FRC
161 +
162 STO 02
163 FRC
164 RCL 01
165 1
166 +
167 3
168 *
169 RCL 09
170 *
171 RCL 05
172 +
173 +
174 CLRGX
175 RCL IND 02
176 STO IND 05
177 CLX
178 STO 06
179 RCL 04
180 ENTER
181 SIGN
182 ST- Y
183 STO IND Y
184 ST+ 02
185 LBL 05

186 ST+ X
187 RCL IND 04
188 RCL 00
189 *
190 -
191 RCL 06
192 ST+ X
193 ENTER
194 SIGN
195 +
196 *
197 RCL IND 04
198 X<> IND 03
199 RCL 06
200 ST* Y
201 RCL 00
202 +
203 ENTER
204 SIGN
205 ST+ 03
206 +
207 *
208 +
209 RCL 06
210 RCL 06
211 RCL 00
212 -
213 ST* Y
214 +
215 /
216 X<> IND 04
217 ISG 04
218 GTO 05
219 RCL 07
220 STO 03
221 RCL 08

222 RCL 05
223 RCL IND 02
224 SIGN
225 DSE Z
226 LBL 06
227 CLX
228 RCL IND Z
229 LASTX
230 *
231 ST+ IND Y
232 ISG Y
233 CLX
234 ISG Z
235 GTO 06
236 RCL 08
237 RCL 09
238 +
239 STO 04
240 STO 08
241 SIGN
242 ST+ 06
243 ISG 02
244 GTO 05
245 RCL 02
246 RCL 01
247 -
248 1
249 -
250 RCL 03
251 2
252 -
253 RCL 09
254 *
255 RCL 05
256 +
257 END

( 348 bytes / SIZE N+4m+15 )

STACK	INPUTS	OUTPUTS
Y	/	bbb.eee(a_i)
X	m	bbb.eee(c_i)

Example:

t	0	1	2	3	4	5	6
y	1.2	3.4	5.1	6.7	7.6	8	8.4

N = 6 STO 00

1.2 STO 10         6.7 STO 13        8.4 STO 16
3.4 STO 11         7.6 STO 14
5.1 STO 12           8   STO 15

-If we seek a polynomial of degree 4

4 XEQ "LSPN" >>>> 22.026 = control number of the coefficients c_i's --- execution time = 1mn10s ---
RDN 17.021 = control number of the coefficients a_i's

and we get:

p(t) = 5.771429 - 3.567857 P_1,N(t) - 1.005952 P_2,N(t) - 0.016667 P_3,N(t) + 0.037013 P_4,N(t)
p(t) = 1.217965 + 2.088023 t + 0.093561 t² - 0.083586 t³ + 0.007197 t⁴

Notes:

-With m = N , youhave the Lagrange Polynomial.
-With N = 87 and m = 13 execution time = 36m43s ( about 3s with V41 )

3°) Root-Mean-Square Error

"RMS" evaluates the root-mean-square error = SQRT [ Sum_k=0,...,N ( T_k - A_k )² / N ]

where T = true value and A = approximate value

Data Registers: • R00 = N ( N+1 points ) ( Registers R00 and R10 thru R10+N are to be initialized before executing "LSPN" & "RMS" )

R01 = RMSE

• R10 = y₀ • R11 = y₁ • R12 = y₂ ......................... • R10+N = y_N

R12+N+m = c₀ , R13+N+m = c₁ , ........................... , R12+N+2m = c_m

Flags: /
Subroutines: /

01 LBL "RMS"
02 INT
03 LASTX
04 FRC
05 E3
06 ST/ Z
07 *
08 +
09 STO 03

10 CLX
11 STO 01
12 STO 02
13 RCL 00
14 E3
15 /
16 10.01
17 +
18 STO 04

19 LBL 01
20 RCL 03
21 ENTER
22 CLX
23 LBL 02
24 RCL IND Y
25 +
26 RCL 02
27 *

28 DSE Y
29 GTO 02
30 RCL IND Y
31 +
32 RCL IND 04
33 -
34 X^2
35 ST+ 01
36 SIGN

37 ST+ 02
38 ISG 04
39 GTO 01
40 RCL 01
41 RCL 00
42 /
43 SQRT
44 STO 01
45 END

( 69 bytes )

STACK	INPUTS	OUTPUTS
X	bbb.eee	rms error

where bbb.eee is the control number of the polynomial c₀ + c₁ t + c₂ t² + ...... + c_m t^m

Example: With the example of the first 2 paragraphs and after executing "LSPN"

22.026 XEQ "RMS" >>>> RMSE = 0.0656 ---Execution time = 10s---

Notes:

-RMSE may also be defined by SQRT [ Sum_k=0,...,N ( T_k - A_k )² / ( N + 1 ) ]
-In this case, add 1 + after line 41 ( the result becomes 0.0608 in the example above )

-With N = 87 & m = 13, execution time = 5mn09s

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] F. Scheid - "Numerical Analysis" - Mc Graw Hill ISBN 07-055197-9