# 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 ( x0 , y0 ) ,  ( x1 , y1 )  , ........................ ,  ( xN , yN )
-We define orthogonal polynomials Pk,N(t)  by                                                                                  ( t = (x-x0)/h  with  h = xi+1 - xi )

Pk,N(t) = SUM i=0,1,...,k (-1)i  Cki Ck+ii  t(i)/N(i)

where   t(i) = t(t-1)(t-2)....(t-i+1)  and  N(i) = N(N-1)(N-2).....(N-i+1)      ( = 1 if  i = 0 )
and    Cki = k!/(i!(k-i)!)

-The polynomial  p(x) can then be written  p(x) = a0 P0,N(t) + a1 P1,N(t) + .................... + am Pm,N(t)

with   ak = [ SUM t=0,1,...,N  yt Pk,N(t) ]/[ SUM t=0,1,...,N  (Pk,N(t))2 ]

-These coefficients are computed and stored when the program reaches line 149
-Then, the program calculates the coefficients in the more standard form:   c0 + c1 t + c2 t2 + ...... + cm tm

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

P-1,N(t) = 0 ,  P0,N(t) = 1 ,  (k+1)(N-k) Pk+1,N(t) = (2k+1)(N-2t) Pk,N(t) - k(N+k+1) Pk-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 = y0  •  R11 = y1  •  R12 = y2   .........................   •  R10+N = yN

>>>>    at the end          R11+N = a0  ,  R12+N =  a1  ,  .............................. ,  R11+N+m = am
R12+N+m = c0  , R13+N+m = c1  , ........................... ,  R12+N+2m = cm

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(ai) X m bbb.eee(ci)

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  ci's            --- execution time = 1mn25s ---
RDN   17.021 = control number of the coefficients  ai's

and we get:

p(t) = 5.771429 - 3.567857 P1,N(t) - 1.005952 P2,N(t) - 0.016667 P3,N(t) + 0.037013 P4,N(t)
p(t) = 1.217965 + 2.088023 t + 0.093561 t2 - 0.083586 t3 + 0.007197 t4

Notes:

-With m = N , you have the Lagrange Polynomial.
-This program calculates the Stirling number of the 1st kind to get the ci's from the ai'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  Pk+1,N(t)  Pk,N(t) & Pk-1,N(t) to compute the ci'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 = y0  •  R11 = y1  •  R12 = y2   .........................   •  R10+N = yN

>>>>    at the end          R11+N = a0  ,  R12+N =  a1  ,  .............................. ,  R11+N+m = am
R12+N+m = c0  , R13+N+m = c1  , ........................... ,  R12+N+2m = cm

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(ai) X m bbb.eee(ci)

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  ci's                        --- execution time = 1mn10s ---
RDN   17.021 = control number of the coefficients  ai's

and we get:

p(t) = 5.771429 - 3.567857 P1,N(t) - 1.005952 P2,N(t) - 0.016667 P3,N(t) + 0.037013 P4,N(t)
p(t) = 1.217965 + 2.088023 t + 0.093561 t2 - 0.083586 t3 + 0.007197 t4

Notes:

-With m = N , you have 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 [ Sumk=0,...,N ( Tk - Ak )2 / 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 = y0  •  R11 = y1  •  R12 = y2   .........................   •  R10+N = yN

R12+N+m = c0  , R13+N+m = c1  , ........................... ,  R12+N+2m = cm

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   c0 + c1 t + c2 t2 + ...... + cm tm

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 [ Sumk=0,...,N ( Tk - Ak )2 / ( 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