hp41programs

MinMaxParabola Min-Max Polynomial for the HP-41

Overview

1°) Min-Max Parabola
2°) Min-Max Polynomial

1°) Min-Max Parabola

-Given a set of n data points ( x₁ , y₁ ) ..... ( x_n , y_n ) ( n > 3 )
the min-max parabola: y = a.x² + b.x + c minimizes the maximum "error" h = sup_i | y_i - a x_i² - b x_i - c | over the entire set of points.

-The following program uses the "exchange method":

1°) The min-max parabola is found for the first 4 points.
2°) The greatest error is computed at all data points for this parabola.
3°) One of the first 4 points is exchanged with a data point at which the greatest error occurs ( so that the errors are still of alternating sign for this new quadruple )

-The process is repeated until the greatest error occurs at one of the first 4 points ( step 3 is skipped if it happens the first time )

Data Registers: • R00 = n = number of points ( Registers R00 thru R2n are to be initialized before executing "MMP2" )

• R01 = x₁ • R03 = x₂ .......... • R2n-1 = x_n ( it's necessary to have x₁ < x₂ < x₃ < x₄ )
• R02 = y₁ • R04 = y₂ .......... • R2n = y_n

Flag: F10 is cleared at the end
Subroutines: /

-Lines 122-160-175-186-196-206-214 are three-byte GTOs

01 LBL "MMP2"
02 LBL 01
03 RCL 06
04 RCL 02
05 -
06 RCL 07
07 RCL 03
08 -
09 *
10 RCL 08
11 RCL 04
12 -
13 RCL 05
14 RCL 01
15 -
16 *
17 -
18 RCL 01
19 RCL 03
20 -
21 RCL 07
22 -
23 RCL 05
24 ST+ Y
25 RCL 01
26 -
27 *
28 RCL 07
29 RCL 03
30 -
31 *
32 /
33 STO M
34 RCL 01
35 RCL 05
36 +
37 *

38 RCL 06
39 RCL 02
40 -
41 RCL 05
42 RCL 01
43 -
44 /
45 X<>Y
46 -
47 STO N
48 RCL 03
49 RCL 05
50 +
51 *
52 RCL 03
53 X^2
54 RCL 05
55 X^2
56 +
57 RCL M
58 *
59 +
60 RCL 04
61 RCL 06
62 +
63 X<>Y
64 -
65 2
66 /
67 STO O
68 RCL 01
69 RCL M
70 *
71 RCL N
72 +
73 RCL 01
74 *

75 +
76 RCL 02
77 -
78 CF 10
79 X<0?
80 SF 10
81 ENTER^
82 ABS
83 RCL 00
84 ST+ X
85 E3
86 /
87 8
88 +
89 STO P
90 STO Q
91 LBL 02
92 ISG P
93 X<0?
94 GTO 03
95 CLX
96 RCL IND P
97 RCL M
98 *
99 RCL N
100 +
101 RCL IND P
102 *
103 RCL O
104 +
105 ISG P
106 RCL IND P
107 -
108 ABS
109 X<=Y?
110 GTO 02
111 LASTX

112 X<>Y
113 RCL P
114 STO Q
115 GTO 02
116 LBL 03
117 X<> Q
118 INT
119 STO L
120 8
121 X=Y?
122 GTO 09
123 R^
124 X<0?
125 FS?C 10
126 FS? 30
127 SF 10
128 LASTX
129 DSE X
130 STO Y
131 RCL IND X
132 RCL 07
133 X<Y?
134 GTO 04
135 CLX
136 RCL 05
137 X<Y?
138 GTO 05
139 CLX
140 RCL 03
141 X<Y?
142 GTO 06
143 CLX
144 RCL 01
145 X<Y?
146 GTO 07
147 FC? 10
148 GTO 08

149 X<> 03
150 X<> 05
151 X<> 07
152 X<> IND T
153 STO 01
154 RCL 02
155 X<> 04
156 X<> 06
157 X<> 08
158 X<> IND L
159 STO 02
160 GTO 01
161 LBL 04
162 FS? 10
163 GTO 00
164 X<> 05
165 X<> 03
166 X<> 01
167 X<> IND T
168 STO 07
169 RCL 02
170 X<> IND L
171 X<> 08
172 X<> 06
173 X<> 04
174 STO 02
175 GTO 01
176 LBL 05
177 FC? 10
178 GTO 06
179 LBL 00
180 RCL 07
181 X<> IND T
182 STO 07
183 RCL 08
184 X<> IND L
185 STO 08

186 GTO 01
187 LBL 06
188 FS? 10
189 GTO 00
190 RCL 05
191 X<> IND T
192 STO 05
193 RCL 06
194 X<> IND L
195 STO 06
196 GTO 01
197 LBL 07
198 FS? 10
199 GTO 00
200 LBL 08
201 X<> IND T
202 STO 01
203 RCL 02
204 X<> IND L
205 STO 02
206 GTO 01
207 LBL 00
208 RCL 03
209 X<> IND T
210 STO 03
211 RCL 04
212 X<> IND L
213 STO 04
214 GTO 01
215 LBL 09
216 R^
217 RCL O
218 RCL N
219 RCL M
220 CF 10
221 CLA
222 END

( 321 bytes / SIZE 2n+1 )

STACK	INPUTS	OUTPUTS
T	/	+/-h
Z	/	c
Y	/	b
X	/	a

Example: With the following set of n = 11 points ( 11 STO 00 )

x	0	1	2	3	4	5	6	7	8	9	10
y	1.2	6	9.7	13.7	16.8	20	21.7	24	24.9	25.7	26

-Store these 22 numbers into

R01 R03 .............................................................................................. R21
R02 R04 .............................................................................................. R22 respectively

           XEQ "MMP2" yields          a = -0.25          ( in 30 seconds )
                                    RDN          b = 5
                                    RDN          c = 0.975
                                    RDN    +/- h = -0.275

>>> whence the equation of the min-max parabola is y = -0.25 x² + 5 x + 0.975 and the maximum "error" is 0.275

Notes:

-When the program stops, the order of the points will have probably changed.
-The abscissas can be unequally spaced.
-Due to roundoff errors, it's perhaps not impossible that the algorithm oscillates for ever between 2 very close | h |-values
-Add VIEW Z after line 120 to check
-If it ever happens, set flag F21 ( the program will stop each time VIEW Z is executed ) and press XEQ 09 to get the results.

2°) Min-Max Polynomial of degree d

-The same "exchange method" may also be used to find the min-max polynomial of degree d
-The data must contain at least (d+2) points

WARNING:

-The set of n data points ( x₁ , y₁ ) ..... ( x_n , y_n ) must be stored so that x₁ < x₂ < x₃ < ............... < x_d+2

Data Registers: • R00 = n = number of points < 45 ( Registers R00 thru R2n are to be initialized before executing "MMP" )

• R01 = x₁ • R03 = x₂ .......... • R2n-1 = x_n with x₁ < x₂ < x₃ < ............... < x_d+2
• R02 = y₁ • R04 = y₂ .......... • R2n = y_n

R89 to R99: temp R100 to R105+d^2+5.d = the coefficients of the successive linear systems

Flag: CF01 <-> pivoting
SF 01 <-> no pivoting

Subroutine: "LS" ( cf "Linear and non-linear systems for the HP-41" )

-Lines 121-169-214-235-245-273 are three-byte GTOs

01 LBL "MMP"
02 STO 97
03 STO 99
04 X^2
05 LASTX
06 5
07 *
08 +
09 105
10 +
11 STO 98
12 2
13 ST+ 97
14 RCL 97
15 1
16 %
17 RCL 98
18 +
19 .1
20 %
21 100
22 +
23 STO 95
24 RCL 97
25 ENTER
26 X^2
27 .1
28 %
29 +
30 DSE Y
31 +
32 100.1
33 +
34 STO 94
35 RCL 00
36 STO 96
37 GTO 12
38 LBL 11
39 RCL IND 89
40 X<> IND 92
41 STO IND 89
42 CLX
43 SIGN

44 ST+ 89
45 ST+ 92
46 RCL IND 89
47 X<> IND 92
48 STO IND 89
49 LBL 12
50 RCL 98
51 RCL 97
52 ST+ X
53 2 E-5
54 +
55 LBL 00
56 RCL IND X
57 STO IND Z
58 RDN
59 DSE Y
60 DSE X
61 GTO 00
62 FRC
63 RCL 97
64 ST+ X
65 DSE X
66 +
67 STO N
68 CLX
69 RCL 99
70 STO O
71 LBL 01
72 CLX
73 RCL IND N
74 RCL O
75 Y^X
76 STO IND Y
77 DSE Y
78 DSE N
79 GTO 01
80 CLX
81 RCL N
82 CHS
83 RCL 97
84 DSE X
85 ST+ X
86 +

87 STO N
88 DSE O
89 GTO 01
90 RCL 97
91 STO T
92 SIGN
93 LBL 02
94 STO IND Z
95 DSE Z
96 DSE L
97 GTO 02
98 LBL 03
99 CHS
100 STO IND Z
101 DSE Z
102 DSE T
103 GTO 03
104 RCL 95
105 0
106 XEQ "LS"
107 RCL 96
108 STO 00
109 CLX
110 STO 91
111 STO 92
112 RCL 00
113 ST+ X
114 .1
115 %
116 RCL 97
117 ST+ X
118 +
119 ISG X
120 X=0?
121 GTO 10
122 STO 93
123 STO M
124 SIGN
125 RCL 94
126 INT
127 +
128 RCL 97
129 -

130 RCL IND X
131 STO 90
132 VIEW X
133 LBL 04
134 RCL 94
135 ENTER
136 CLX
137 LBL 05
138 RCL IND M
139 *
140 RCL IND Y
141 +
142 DSE Y
143 GTO 05
144 ISG M
145 RCL IND M
146 -
147 STO O
148 ABS
149 RCL 91
150 ABS
151 X>Y?
152 GTO 04
153 RCL O
154 STO 91
155 RCL M
156 STO 92
157 LBL 04
158 ISG M
159 GTO 04
160 DSE 92
161 CLX
162 RCL 90
163 ABS
164 E-5
165 +
166 RCL 91
167 ABS
168 X<=Y?
169 GTO 10
170 RCL 93
171 INT
172 STO Y

173 2
174 -
175 RCL IND 92
176 ENTER
177 LBL 06
178 CLX
179 RCL IND Z
180 X<Y?
181 GTO 06
182 X<> Z
183 STO T
184 X<> Z
185 LBL 06
186 DSE Z
187 CLX
188 DSE Z
189 GTO 06
190 R^
191 2
192 -
193 STO N
194 ABS
195 STO 89
196 RCL 94
197 ENTER
198 CLX
199 LBL 07
200 RCL IND 89
201 *
202 RCL IND Y
203 +
204 DSE Y
205 GTO 07
206 ISG 89
207 CLX
208 RCL IND 89
209 -
210 RCL 91
211 *
212 DSE 89
213 X>0?
214 GTO 11
215 1

216 CHS
217 RCL N
218 X#Y?
219 GTO 14
220 RCL 92
221 INT
222 1
223 +
224 STO Y
225 2
226 -
227 LBL 08
228 RCL IND Y
229 X<> IND Y
230 STO IND Z
231 RDN
232 DSE Y
233 DSE X
234 GTO 08
235 GTO 12
236 LBL 14
237 RCL 93
238 INT
239 2
240 -
241 X=Y?
242 GTO 14
243 LASTX
244 ST+ 89
245 GTO 11
246 LBL 14
247 RCL 01
248 X<> IND 92
249 STO 01
250 RCL 92
251 1
252 +
253 RCL 02
254 X<> IND Y
255 STO 02
256 RCL 93
257 INT
258 DSE X

259 .1
260 %
261 3
262 +
263 1
264 LBL 09
265 RCL IND Y
266 X<> IND Y
267 STO IND Z
268 CLX
269 SIGN
270 +
271 ISG Y
272 GTO 09
273 GTO 12
274 LBL 10
275 RCL 94
276 RCL 00
277 ST+ X
278 1
279 +
280 STO O
281 LBL 13
282 RCL IND Y
283 STO IND Y
284 CLX
285 SIGN
286 +
287 DSE Y
288 GTO 13
289 1
290 -
291 E3
292 /
293 RCL O
294 +
295 RCL IND Y
296 ABS
297 X<>Y
298 CLA
299 CLD
300 END

( 486 bytes / SIZE 106+d^2+5.d )

STACK	INPUTS	OUTPUTS
Y	/	\| h \|
X	d	bbb.eee

   Where       d       = degree of the polynomial
               bbb.eee = control number of the min-max polynomial ( in fact bbb = 2 n + 1 )
     and        | h |     = sup_i | y_i - p( x_i ) |

Example: With the following set of 7 data points, find the min-max polynomial of degree d = 3

x	1	2	4	7	9	10	12
y	2	5	3	7	6	12	17

n = 7 STO 00 CF 01

1 STO 01 2 STO 03 4 STO 05 7 STO 07 9 STO 09 10 STO 11 12 STO 13
2 STO 02 5 STO 04 3 STO 06 7 STO 08 6 STO 10 12 STO 12 17 STO 14

SIZE 130 ( or greater )

3 XEQ "MMP" >>>> The HP-41 displays the successive h-values
( namely 1.32 -1.506666667 -1.705 1.75 ) and eventually

bbb.eee = 15.018 ---Execution time = 4mn43s---
X<>Y | h | = 1.75

-The coefficients of the min-max polynomial of degree 3 are in R15-R16-R17-R18 and we get:

p(x) = 0.033333333 x³ - 0.45 x² + 2.016666666 x + 1.75 = x³ / 30 - 0.45 x² + ( 121 / 60 ) x + 1.75

and | h | = 1.75

Notes:

-Since R89 is used for temporary data storage, the number of points must be < 45

-When the program stops, the order of the points will have probably changed.

-The abscissas can be unequally spaced.

-Due to roundoff errors and to avoid a possible infinite loop, I've added lines 164-165
-So, the iterations will stop if 2 consecutive h-values differ by less than 0.00001
-Change these lines ( or delete them ) if you prefer another value...

-If you set flag F01 ( SF 01 ), the linear systems will be solved without pivoting.
-With the example above, you will get ( almost ) the same results in 4m15s, thus saving 28 seconds.
-However, there is a risk of more roundoff-errors.

-If d = 2, "MMP2" is much faster and for d = 1, cf "Min-Max Line for the HP-41" ( though "MMP" will work too... )

Reference:

[1] F. Scheid - "Numerical Analysis" - Mc Graw Hill ISBN 07-055197-9