hp41programs

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
 
 

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

 
Example:

 
  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 , you have 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
 
 

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

 
Example:

 
  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 , 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 [ 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: /
 
 

 
    ( 69 bytes )
 
 

   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