The Museum of HP Calculators

# Rational Function Interpolation for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

## Overview

-The following program uses Thiele's formula:

-Given n data points  ( x1 , y1 ) , ( x2 , y2 ) , .................. , ( xn , yn ) ,   the rational function f(x) is evaluated by the continued fraction:

x - x1     x - x2      x - x3                                       x - xn-1
f(x) =   y1 +  ------    -------    --------   ...............................  ---------
r12 +      r123 +     r1234 +                                      r12......n

where  r1k = ( xk - x1 )/( yk - y1 )  ,    r12k = ( xk - x2 )/( r1k - r12 )  ,   r123k = ( xk - x3 )/( r12k - r123 )  , .............. ,  r12....n = ( xn - xn-1 )/( r12...n - r12...n-1 )

-Thus,  f(x) = p(x)/q(x)     with   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 = x1    •  Rn+1 = y1
•  R02 = x2    •  Rn+2 = y2       R2n+1 = r12
---------        ----------         R2n+2 = r13     R3n = r123
---------        ----------         ------------       ---------

•  Rnn = xn    •  R2n   = yn        R3n-1 = r1n     R4n-3 = r12n  ---------  R(n2+3n)/2 = r12.....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: /

01  LBL "THIELE"
02  STO M
03  RCL 00
04  RCL 00
05  3
06  +
07  *
08  2
09  /
10  FS? 01
11  GTO 03
12  CF 10
13  FS? 02
14  SF 10
15  RCL 00
16   E3
17  /
18  STO O
19  .999
20  +
21  STO N
22  SIGN
23  RCL 00
24  +
25  STO P                ( synthetic )
26  STO Q
27  RCL 00
28  +
29  ISG Q
30  LBL 01
31  RCL N
32  INT
33  ISG X
34  CLX
35  RCL O
36  FRC
37  +
38  STO O
39  X<>Y
40  LBL 02
41  RCL IND N
42  RCL IND O
43  -
44  RCL IND P
45  FS? 10
46  1/X
47  RCL IND Q
48  FS? 10
49  1/X
50  -
51  /
52  STO IND Y
53  CLX
54  SIGN
55  ST+ Q
56  +
57  ISG O
58  GTO 02
59  CF 10
60  RCL Q
61  STO P                ( synthetic )
62  SIGN
63  ST+ Q
64  X<>Y
65  ISG N
66  GTO 01
67  LASTX
68  LBL 03
69  STO O
70  RCL 00
71  STO N
72  SIGN
73  ST- N
74  0
75  LBL 04
76  RCL IND O
77  +
78  RCL M
79  RCL IND N
80  -
81  X<>Y
82  /
83  X<>Y
84  1
85  +
86  ST- O
87  X<>Y
88  DSE N
89  GTO 04
90  RCL IND O
91  FS? 02
92  1/X
93  +
94  FS? 02
95  1/X
96  RCL M
97  SIGN
98  X<>Y
99  CLA
100  END

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

 STACK INPUTS OUTPUTS X x f(x) L / x

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      into registers     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 r12...-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  y1 equals another y-value, there will be a DATA ERROR , so always choose the first point such that y1 is different from all other yk
-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  and  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 ( xi , yi ) 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 ( xi , 1/yi ) 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:        F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
Abramowitz and Stegun  -  "Handbook of Mathematical Functions"  -  Dover Publications  -  ISBN  0-486-61272-4