Rational Function Interpolation for the HP-41

Overview
 

-The following program uses Thiele's formula:

-Given n data points  ( x₁ , y₁ ) , ( x₂ , y₂ ) , .................. , ( x_n , y_n ) ,   the rational function f(x) is evaluated by the continued fraction:

                       x - x₁     x - x₂      x - x₃                                       x - x_n-1
  f(x) =   y₁ +  ------    -------    --------   ...............................  ---------
                       r₁₂ +      r₁₂₃ +     r₁₂₃₄ +                                      r_12......n
 

  where  r_1k = ( x_k - x₁ )/( y_k - y₁ )  ,    r_12k = ( x_k - x₂ )/( r_1k - r₁₂ )  ,   r_123k = ( x_k - x₃ )/( r_12k - r₁₂₃ )  , .............. ,  r_12....n = ( x_n - x_n-1 )/( r_12...n - r_12...n-1 )
 

-Thus,  f(x) = p(x)/q(x)   where  deg(p) = deg(q)         if  n is odd
                                        and   deg(p) = deg(q) + 1   if n is even      ( the degrees may of course be smaller if some leading coefficients equal zero )

-However, if n is even, you may prefer:   deg(p) = deg(q) - 1
-In this case, simply set flag F02 and the above formula will be applied to 1/y
-But this will work only if all y-values are different from zero.

-"THIELE" uses synthetic registers M N O P Q  which can be replaced by any unused data registers.
 

Program Listing
 

Data Registers:           •  R00 = n                                                     ( Registers R00 thru R2n  are to be initialized before executing "THIELE"  )

                                      •  R01 = x₁    •  Rn+1 = y₁
                                      •  R02 = x₂    •  Rn+2 = y₂       R2n+1 = r₁₂
                                         ---------        ----------         R2n+2 = r₁₃     R3n = r₁₂₃
                                         ---------        ----------         ------------       ---------

                                      •  Rnn = x_n    •  R2n   = y_n        R3n-1 = r_1n     R4n-3 = r_12n  ---------  R(n²+3n)/2 = r_12.....n

Flags:  F01-F02-F10   Set flag F01 after executing the routine once: it will avoid to re-calculate the reciprocal differences.
                                     Set flag F02 if you prefer that the degree of the numerator is smaller than the degree of the denominator ( if n is even )
Subroutines: /
 
 

 
   ( 159 bytes SIZE (n²+3n+2)/2 )
 
 

 
Example1:     Given  f(0) = 226 ,  f(2) = 58 ,  f(5) = 18 ,  f(10) = 6 ,  f(20) = 1     Calculate  f(1) ,  f(3)  ,  f(7)

-There are 5 points so   5  STO 00  and we store:

    0   226             R01   R06
    2    58              R02   R07
    5    18      in     R03   R08    respectively
   10    6               R04   R09
   20    1               R05   R10

-Then,  CF 01    1  XEQ "THIELE"   >>>>   f(1) =  109.4020   ( in 12 seconds )

-After executing "THIELE" once, the r_12...-values are in the required registers ( R11 R15 R18 R20 in this example ),
 and we can set flag F01 to speed up execution:

  SF 01   3   R/S  >>>>   f(3) =  35.8827   ( in 3 seconds )
              7   R/S  >>>>   f(7) =  10.8845

-Here, the number of points is odd ( n = 5 ), so setting or clearing flag F02 yields the same results ( provided no y-value equals zero )
 

Example2:     Given  f(0) = 226 ,  f(2) = 58 ,  f(10) = 6 ,  f(20) = 1      Evaluate again f(1) , f(3) , f(7)

-There are 4 points whence   4  STO 00  and

    0    STO 01    2   STO 02   10   STO 03   20  STO 04
  226  STO 05   58  STO 06    6    STO 07    1   STO 08

 1°)  With flag F02 clear CF 02, we find:

   CF 01     1  XEQ "THIELE"  >>>>   f(1) = 97.5178 ,  SF 01   3  R/S  >>>>   f(3) = 38.9928 ,  7  R/S  >>>>  f(7)  = 12.2710

 2°)  With flag F02 set SF 02, we find:

   CF 01     1  R/S  >>>>    f(1) = 100.1803 ,   SF 01   3  R/S  >>>>   f(3) = 37.8940 ,  7  R/S  >>>>  f(7)  = 11.5319
 

Notes:

-If  y₁ equals another y-value, there will be a DATA ERROR , so always choose the first point such that y₁ is different from all other y_k
-There is always a small risk that some denominators equal 0. If it ever happens, change the order of the data points and start again...

-This program only deals with the cases:

  deg(numerator) = deg(denominator)          if n is odd
  deg(numerator) = deg(denominator) +/-1  if n is even,  but there are many other possibilities.

-For instance, if you want to fit a set of data points ( x_i , y_i ) to a rational function of the type  1/p(x)  where p is a polynomial,
 you can use Lagrange's interpolation formula to the set ( x_i , 1/y_i ) and take the reciprocal of the results.

-In the first example above ( 5 points ) it gives  f(1) = 108.4804  ,   f(3) = 35.8694  ,  f(7) = 10.9406
 

References:

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