hp41programs

Lagrange Lagrange Polynomials for the HP-41
 

Overview
 

 1°) One Dimension Problem

   a) Univariate Interpolation
   b) Derivative of the Lagrange Polynomial
   c) Integration of the Lagrange Polynomial
   d) 7 Points with Equidistant Abscissas: f(x) = 0 & f'(x) = 0

 2°) Bivariate Interpolation

   a) non-synthetic program
   b) synthetic program


Latest Update:   §1°)d)

 

1°) One Dimension Problem
 

    a) Univariate Interpolation
 

-Let be a set of n data points:   M1 ( x1 ; y1 )  ,   M2 ( x2 ; y2 )  ,  ..........  ,  Mn ( xn ; yn )
 There is a unique polynomial L of degree < n  such that:     L(xi) = yi       i = 1 ; 2 ; ..... ; n
-In other words, L is the collocation polynomial for these arguments.
 ( xi values can be unequally-spaced.)

-The following program calculates  L(x) = y1 L1(x) + y2 L2(x) + ..... + yn Ln(x)

 where    Li(x) = ( x - x1 ) ...... ( x - xi-1 ) ( x - xi+1 ) ...... ( x - xn ) / [ ( xi - x1 ) ...... ( xi - xi-1 ) ( xi - xi+1 ) ( xi - xn ) ]
 are the Lagrange's multiplier functions.

-This formula is also applied to an extrapolation to the limit in example2.
 

Data Registers:

-You have to store the 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  into contiguous registers,  starting at any register  Rbb:
 

       x1         x2    ...........        xn
       y1         y2    ...........        yn

                                  into
 

        Rbb      Rbb+2      .............      Ree-1
      Rbb+1      Rbb+3      .............        Ree

                      ( ee = bb + 2n -1 )

Flags:  /
Subroutines:  /

-Thanks to synthetic programming,  "LAGR" only uses the registers containing the  xi and yi numbers,
  but you can replace status registers M N O P Q by the standard R00 R01 R02 R03 R04 ( in this case, add  CLX  STO 03  after line 06 and of course, choose bb > 04 ).
 
 

01  LBL "LAGR"
02  CLA
03  STO M         
04  X<>Y
05  STO O
06  STO Q
07  LBL 01
08  RCL Q		 09  STO N
10  SIGN
11  LBL 02
12  RCL M         
13  RCL IND O
14  RCL IND N
15  ST- Z
16  -		 17  X#0?
18  GTO 03
19  SIGN
20  STO Y
21  LBL 03         
22  /
23  *
24  ISG N		 25  ISG N          
26  GTO 02
27  ISG O
28  RCL IND O
29  *
30  ST+ P
31  ISG O
32  GTO 01		 33  RCL Q 
34  RCL M         
35  SIGN
36  X<> P
37  CLA
38  END

 
  ( 69 bytes / SIZE eee+1 )
 
 

    STACK      INPUTS    OUTPUTS
         Y       bbb.eee      bbb.eee
         X            x        L(x)
         L            /          x

 
Example1:     Evaluate  y(3) and y(5) for the collocation polynomial that takes the values prescribed below:
 
 

      xi       0       1       2       4       7
      yi       3       2       4       6       5

 
-For instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

 1.010  ENTER^
     3     XEQ "LAGR"  produces:  y(3) = 5.847619048    ( in 14 seconds )

 RDN    5   R/S     >>>>     y(5) = 4.523809520
 

Example2:

-Lagrange's polynomial can also be used for an extrapolation to the limit:
 Suppose we are using the trapezoidal rule to evaluate the integral:       I  =  §02  ex dx
 The interval  [0;2]  is divided into n subintervals and

  n = 1  yields   I1 = 8.389056101
  n = 2  yields   I2 = 6.912809880
  n = 4  yields   I4 = 6.521610110
  n = 8  yields   I8 = 6.422297820

We would like to know the result with n = infinity !  i-e  1/n = 0.
The error of the trapezoidal rule is nearly proportional to  1/n2  , therefore, we can use Lagrange interpolation formula with

   x1 = 1   x2 = 1/4   x3 = 1/16   x4 = 1/64
   y1 = I1   y2 = I2      y3 = I4       y4 =  I8      and  evaluate   L (0)

If these 8 numbers are stored in R11 thru R18 for instance,

  11.018  ENTER^
      0       XEQ "LAGR"  produces   6.389056387

Exact result is of course e2 - 1 = 6.389056099 to ten places
and the error of our Lagrange approximation is only  3. 10-7  although I8 was only exact to 1 decimal!

N.B:

- If you are using Simpson's rule instead of the trapezoidal rule, store  xn = 1/n4  : it's a 4th order method ( at least far from any singularity )
-  n is usually doubled at every step, but it's not necessary here.
-  The same idea can be applied to any other problem of this kind.
 

    b) Derivative of the Lagrange Polynomial
 

-"DLAGR" takes x in X-register and returns  L'(x)  in X-register.
 

Data Registers:

-The same 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  are to be stored into contiguous registers,  starting at any register  Rbb:
 

       x1         x2    ...........        xn
       y1         y2    ...........        yn

                                  into
 

      Rbb    Rbb+2   .............    Ree-1
    Rbb+1    Rbb+3   .............      Ree

                      ( ee = bb + 2n -1 )

Flags:  /
Subroutines:  /
 
 

 01  LBL "DLAGR"
 02  CLA
 03  STO M           
 04  X<>Y
 05  STO Q
 06  STO a
 07  LBL 01
 08  RCL a
 09  STO O
 10  CLX
 11  LBL 02
 12  RCL O
 13  RCL Q
 14  X=Y?
 15  GTO 02		  16  X<> Z
 17  RCL a
 18  STO N           
 19  SIGN
 20  LBL 03
 21  RCL N
 22  RCL O
 23  X=Y?
 24  GTO 03
 25  CLX
 26  RCL Q
 27  X=Y?
 28  GTO 03
 29  RDN
 30  CLX		  31  RCL M           
 32  RCL IND N
 33  -
 34  *
 35  ENTER^
 36  ENTER^
 37  LBL 03
 38  R^
 39  R^
 40  ISG N
 41  ISG N
 42  GTO 03
 43  +
 44  STO Z
 45  LBL 02		  46  X<> Z
 47  ISG O
 48  ISG O
 49  GTO 02          
 50  RCL a
 51  X<>Y
 52  LBL 04
 53  RCL IND Q
 54  RCL IND Z
 55  -
 56  X=0?
 57  SIGN
 58  /
 59  ISG Y
 60  ISG Y		  61  GTO 04
 62  ISG Q
 63  RCL IND Q
 64  *
 65  ST+ P
 66  ISG Q
 67  GTO 01          
 68  RCL a
 69  RCL M
 70  SIGN
 71  X<> P
 72  CLA
 73  END

 
     ( 123 bytes / SIZE eee+1 )
 
 

    STACK      INPUTS    OUTPUTS
         Y       bbb.eee      bbb.eee
         X            x        L' (x)
         L            /          x

 
Example:

  Evaluate  y(PI)  for the collocation polynomial that takes the values prescribed below:
 
 

    xi     0     1     2     4     7
    yi     3     2     4     6     5

 
-If, for instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

 1.010  ENTER^
    PI    XEQ "DLAGR"  >>>>  L' (PI) = 0.857353871                                  ---Execution time = 67s---

Notes:

-Synthetic register a is used, so this program cannot be called as more than a first level subroutine.
-Registers M N O P Q a  may of course be replaced by standard registers like R00 thru R05
-In this case, choose  bbb > 005.
 

    c) Integration of the Lagrange Polynomial
 

-This program is also listed in "Numerical Integration for the HP-41"

 n data points are given and stored in registers  R01 thru R2n
 
 

       x1         x2    ...........        xn
       y1         y2    ...........        yn

 
-The program below calculates its coefficients c0 , c1 , ........... , cn-1  when L(x) is written:

    L(x) = c0 + c1 ( x-x1 ) + c2 ( x-x1 )2 + .................. + cn-1 ( x-x1 )n-1

-First, x1 is subtracted from all other xi 's
-We have obviously  c0 = y1  and  the other yi 's  are replaced by  ( yi - y1 ) / ( xi - x1 )
-Then, we find the value of the new Lagrange polynomial at zero, thus producing c1 which is stored in register R04
-The process is repeated until  cn-1 is computed and stored in register R2n
-Finally, lines 42 to 57 evaluate  §x1xn  L(x).dx
 

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

                                  •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
                                  •  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

  >>>>    when the program stops,  R02 = c0 = y1 ,  R04 = c1 , ............... ,  R2n = cn-1

Flags:  /
Subroutine:   "LAGR"  ( cf paragrph 1-a) )
 
 

01  LBL "LGI"
02  RCL 00
03  ST+ X
04  STO Y
05  RCL 01         
06  LBL 00
07  DSE Z
08  ST- IND Z 
09  DSE Z
10  GTO 00
11  CLX
12   E3		 13  /
14  ISG X
15  STO 01
16  RCL 02
17  GTO 02         
18  LBL 01
19  2
20  ST+ 01
21  RCL 01
22  0
23  XEQ "LAGR"
24  LBL 02		 25  ISG Y
26  STO IND Y
27  ISG Y
28  FS? 30
29  GTO 04         
30  LBL 03
31  RCL IND Y 
32  ISG Z
33  X<>Y
34  ST- IND Z
35  X<>Y
36  ST/ IND Z		 37  RDN
38  ISG Y
39  GTO 03         
40  GTO 01
41  LBL 04
42  RCL 00
43  ST+ X
44  0
45  LBL 05
46  RCL IND Y 
47  RCL 00
48  /		 49  +
50  RCL IND 01 
51  *
52  2
53  ST- Z
54  RDN
55  DSE 00
56  GTO 05         
57  END

 
   ( 98 bytes / SIZE 2n+1 )
 
 

      STACK        INPUT      OUTPUT
           X              /   §x1xn  L(x).dx

 
Example:     Calculate  §18  L(x).dx   from the following values
 
 

      x      1     2.4      4     5.2      7      8
      y      1      4      6      5      4      2

 
-Store these 12 numbers into registers

     R01  R03  R05  R07  R09  R11   ( x )
     R02  R04  R06  R08  R10  R12    ( y )    respectively

-There are 6 points:     6   STO 00      XEQ "LGI"     >>>>    §18  L(x).dx  =   29.61789480  ( in 43 seconds )

-We have also the coefficients of the Lagrange polynomial ( as a sum of powers of ( x-x1 ) )  in registers  R02  R04  .............  R12
-Here, x1 = 1 and we find:

  L(x) =  1 - 0.362103178 ( x-1 ) + 3.623795356 ( x-1 )2 - 1.661873944 ( x-1 )3 + 0.272598127 ( x-1 )4 - 0.015381483 ( x-1 )5
 

Remarks:

-For large n-values, the coefficients of the Lagrange polynomial are often great in magnitude
  and of alternate signs, thus leading to inaccurate results.
-Moreover, the coefficients themselves are not obtained very accurately if n is too great.

-For instance, with the following data:
 
 

  x   1   3   5   7   9  11  13  15  17  19  21  23  25  27  29  31  33  35  37  39
  y   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36  38  40

 
 "LGI"  yields  797.9971774  wheras the exact result is 798

Execution time t:
 
 

       n        t
       4      15s
       6      43s
      10   2mn51s
      20  21mn13s

 

    d) 7 Points with Equidistant Abscissas:  f(x) = 0 & f'(x) = 0


-The following program uses the Lagrange 7-point interpolation formula.
-We assume that the xi's are equally spaced. ( xi+1 - xi = h  ,   i  = -3 , -2 , -1 , 0 , 1 , 2 )

-So,   f(x0+p.h) - f(x0) = A.p + B.p2 + C.p3 + D.p4 + E.p5 + F.p6   with:

    A = (1/720) [ 540 ( y1 - y-1 ) - 108 ( y2 - y-2 )  + 12 ( y3 - y-3 ) ]
    B = (1/720) [ -980 y0 + 540 ( y1 + y-1 ) - 54 ( y2 + y-2 )  + 4 ( y3 + y-3 ) ]
    C = (1/720) [ -195 ( y1 - y-1 ) + 120 ( y2 - y-2 )  - 15 ( y3 - y-3 ) ]
    D = (1/720) [ 280 y0 - 195 ( y1 + y-1 ) + 60 ( y2 + y-2 )  - 5 ( y3 + y-3 ) ]                         yk = f(xk)
    E = (1/720) [ 15 ( y1 - y-1 ) - 12 ( y2 - y-2 )  + 3 ( y3 - y-3 ) ]
    F = (1/720) [ -20 y0 + 15 ( y1 + y-1 ) - 6 ( y2 + y-2 )  + ( y3 + y-3 ) ]

Data Registers:              R00 = x0            R16 = h             ( Registers R01 thru R07 are to be initialized before executing "SLVD7" )

                                      •  R01 = f(x-3)    •  R02 = f(x-2)    •  R03 = f(x-1)    •  R04 = f(x0)    •  R05 = f(x1)    •  R04 = f(x2)    •  R05 = f(x3)                R08 to R15:  temp

Flag:   F01                       CF 01  to solve f(x) = 0
                                          SF 01  to solve f'(x) = 0   ( extrema )
Subroutines:  / 


 01 LBL "SLVD7"
 02 STO 00
 03 X<>Y
 04 STO 16
 05 RCL 01
 06 RCL 07
 07 +
 08 STO 14
 09 RCL 02
 10 RCL 06
 11 +
 12 6
 13 *
 14 STO 15
 15 -
 16 RCL 03
 17 RCL 05
 18 +
 19 STO 08
 20 15
 21 *
 22 +
 23 RCL 04          
 24 20
 25 *
 26 -
 27 STO 13		  28 RCL 15
 29 ST+ X
 30 RCL 14
 31 -
 32 RCL 08
 33 39
 34 *
 35 -
 36 RCL 04
 37 56
 38 *
 39 +
 40 5
 41 *
 42 STO 11
 43 RCL 14
 44 4
 45 *
 46 RCL 15
 47 9
 48 *
 49 -
 50 RCL 08          
 51 540
 52 *
 53 +
 54 RCL 04		  55 980
 56 *
 57 -
 58 STO 09
 59 RCL 07
 60 RCL 01
 61 -
 62 STO 14
 63 RCL 06
 64 RCL 02
 65 -
 66 STO 15
 67 4
 68 *
 69 -
 70 RCL 05
 71 RCL 03          
 72 -
 73 STO 08
 74 5
 75 *
 76 +
 77 3
 78 *
 79 STO 12
 80 RCL 15
 81 8		  82 *
 83 RCL 14
 84 -
 85 RCL 08
 86 13
 87 *
 88 -
 89 15
 90 *
 91 STO 10
 92 RCL 14
 93 RCL 15
 94 9
 95 *
 96 -
 97 RCL 08          
 98 45
 99 *
100 +
101 12
102 *
103 STO 08
104 RCL 04
105 720
106 *
107 FS? 01
108 X<>Y		 109 STO 14
110 RCL 09
111 ST+ X
112 RCL 08
113 FS? 01
114 X<>Y
115 STO 15
116 CLX
117 LBL 01
118 STO Y
119 RCL 13
120 6
121 FC? 01
122 SIGN
123 *
124 *
125 RCL 12        
126 5
127 FC? 01
128 SIGN
129 *
130 +
131 *
132 RCL 11
133 4
134 FC? 01
135 SIGN		 136 *
137 +
138 *
139 RCL 10
140 3
141 FC? 01
142 SIGN
143 *
144 +
145 *
146 RCL 09
147 FS? 01
148 CLX
149 +
150 *
151 FC? 01
152 *
153 RCL 14        
154 +
155 RCL 15
156 CHS
157 /
158 VIEW X
159 X#Y?
160 GTO 01
161 ENTER
162 STO Z		 163 RCL 13
164 *
165 RCL 12
166 +
167 *
168 RCL 11
169 +
170 *
171 RCL 10
172 +
173 *
174 RCL 09        
175 +
176 *
177 RCL 08
178 +
179 *
180 720
181 /
182 RCL 04
183 +
184 X<>Y
185 RCL 16
186 *
187 RCL 00
188 +
189 CLD
190 END


         ( 232 bytes / SIZE 017 )

 
           STACK           INPUTS         OUTPUTS
              Y                h              f(x)
              X               x0               x
 
    If CF 01   x is a solution of f(x) = 0
    If SF 01   x is a solution of f'(x) = 0   ( extremum )

Example1:     f(x) = 0   from these 7 values
  

    x
     1
     3
     5
    7
    9
   11
   13
    y
    16
    10
     6
    1
   -3
   -8
  -20

   16  STO 01   10  STO 02   6  STO 03   1  STO 04  -3  STO 05  -8  STO 06  -20  STO 07

    CF 01           h = 2  &  x0 = 7   so:

    2  ENTER^
    7  XEQ "SLVD7"   >>>>  x = 7.446084649                     ---Execution time = 20s---
                               

Example2:     Extremum from these 7 values
  

    x
    1
    3
     5
    7
    9
   11
   13
    y
   16
   10
     6
    4
    5
    8
   20

   16  STO 01   10  STO 02   6  STO 03   4  STO 04   5  STO 05   8  STO 06   20  STO 07

    SF 01           h = 2  &  x0 = 7   so:

    2  ENTER^
    7  XEQ "SLVD7"   >>>>      x  =  7.235889441                     ---Execution time = 16s---
                                   X<>Y  f(x) =  3.977546065

-So  f(x)  has a minimum  3.977546065  for x = 7.235889441

Notes:

-The method may fail , for example if  A = 0  when CF 01

-It's also possible to have an infinite loop with  X#Y?  ( line 159 )
-It would be better to replace line 159 by the M-Code routine  X#Y??  ( cf "A few M-Code routines for the HP41" )

-Similar programs with 5-point Lagrange polynomials are listed in "Linear & Non-Linear systems for the HP41" ( paragraph 2°)a4) )
  and "Extrema for the HP41" ( paragraph 1°)e) )


2°) Bivariate Interpolation
 

    a) Non-synthetic program
 

-Now we have a grid of  nxm  values:   f(xi,yj)   with   i = 1,2,......,n  &  j = 1,2,.....,m   and we want to estimate f(x,y)

-"LAGR2" uses the Lagrange Polynomial in the directions of x-axis and y-axis.
-You have to store the data in registers Rbb to Ree as shown below ( bb > 10 ):
 

Data Registers:           •  R00 = bbb.eeenn                             ( Registers R00 and Rbb thru Ree are to be initialized before executing "LAGR2"  )

                                  R01 = x    R03 = f(x,y)
                                  R02 = y    R04 thru R10: temp                                             We must choose bb > 10

                                                 x  \  y |        •  Rbb+n = y1                •  Rbb+2n+1 = y2       ..............................    •  Ree-n = ym
                                  ---------------------------------------------------------------------------------------------------------------------
                                •  Rbb         = x1 |   •  Rbb+n+1 = f(x1,y1)    •  Rbb+2n+2 = f(x1,y2)    ..............................     •  Ree-n+1 = f(x1,ym)
                                •  Rbb+1     = x2 |   •  Rbb+n+2 = f(x2,y1)    •  Rbb+2n+3 = f(x2,y2)    ..............................     •  Ree-n+2 = f(x2,ym)
                                   ....................    |     .............................................................................................................................................
                                                           |
                                •  Rbb+n-1 = xn  |   •  Rbb+2n  = f(xn,y1)     •  Rbb+3n+1 = f(xn,y2)   ...............................     •  Ree  = f(xn,ym)

Flags: /
Subroutines: /
 
 

01  LBL "LAGR2"
02  STO 01
03  X<>Y
04  STO 02
05  CLX
06  STO 03          
07  RCL 00
08  ISG X
09  STO 07
10  INT
11  DSE X
12  .1
13  %
14  RCL 00
15  INT
16  +
17  STO 05		 18  STO 06
19   E-5
20  ST+ 07
21  LBL 00
22  RCL 07          
23  STO 08
24  STO 09
25  FRC
26  RCL 06
27  INT
28  +
29  ISG X
30  STO 10
31  CLX
32  STO 04
33  LBL 01
34  RCL IND 10		 35  LBL 02
36  RCL 02
37  RCL IND 08
38  RCL IND 09
39  ST- Z
40  -
41  X#0?
42  GTO 02          
43  SIGN
44  STO Y
45  LBL 02
46  /
47  *
48  ISG 09
49  GTO 02
50  ST+ 04
51  RCL 07		 52  STO 09
53  ISG 10
54  CLX
55  ISG 08
56  GTO 01
57  RCL 05
58  STO 10
59  RCL 04          
60  LBL 03
61  RCL 01
62  RCL IND 06
63  RCL IND 10
64  ST- Z
65  -
66  X#0?
67  GTO 03
68  SIGN		 69  STO Y
70  LBL 03
71  /
72  *
73  ISG 10
74  GTO 03
75  ST+ 03
76  ISG 06
77  GTO 00 
78  RCL 01          
79  SIGN 
80  RCL 02
81  RCL 03
82  END

 
  ( 121 bytes / SIZE? )
 
 

    STACK      INPUTS    OUTPUTS
         Y            y          y
         X            x       f(x,y)
         L            /          x

 
Example:      A function f(x,y) is only known by the following values ( meaning for example:  f(1,4) = 3   f(4,6) = 9  ...etc... )

            >>>   if, for instance, you choose bb = 11 ,  store these  19 numbers  into

     x \ y |   2    3    4    6                                                          |   R14   R18   R22   R26
--------------------------                                              ---------------------------------
       1   |   4    3    3    5                                                  R11  |  R15   R19   R23   R27                             respectively
       2   |   3    1    2    6                                                  R12  |  R16   R20   R24   R28
       4   |   1    0    4    9                                                  R13  |  R17   R21   R25   R29

-The control number of this table is   11.02903  so   11.02903  STO 00
  ( note that nn = 03 is the number of x-values which may be different from m = the number of y-values ( here m = 4 ) )

-Then, to compute, say  f(3,5)

   5  ENTER^
   3  XEQ "LAGR2"  >>>>   f(3,5) ~  5.8333   ( in 29 seconds )

-Likewise, to obtain  f(1.6,2.7)

   2.7  ENTER^
   1.6     R/S        >>>>   f(1.6,2.7) ~  1.9038
 

Note:   Execution time is approximately proportional to nxm
 

    b) Synthetic program
 

-Synthetic programming allows to store the data into registers R01 thru Rn.m+n+m
-However, synthetic register a is used, so this program cannot be called as more than a first level subroutine.
 
 

Data Registers:           •  R00 = n.mmm = n + m/103         ( Registers R00 thru Rnm+m+n are to be initialized before executing "BVI"  )
 

                                          x  \  y |       •  Rnn+1 = y1             •  R2n+2 = y2       ..............................    •  Rnm+m = ym
                                  ---------------------------------------------------------------------------------------------------------------------
                                •  R01  = x1 |   •  Rnn+2 = f(x1,y1)    •  R2n+3 = f(x1,y2)    ..............................     •  Rnm+m+1 = f(x1,ym)
                                •  R02  = x2 |   •  Rnn+3 = f(x2,y1)    •  R2n+4 = f(x2,y2)    ..............................     •  Rnm+m+2 = f(x2,ym)
                                    ...........    |      .............................................................................................................................................
                                                   |
                                •  Rnn = xn  |   •  R2n+1  = f(xn,y1)   •  R3n+2 = f(xn,y2)   ..............................      •  Rnm+m+n  = f(xn,ym)

Flags: /
Subroutines: /
 
 

01  LBL "BVI" 
02  CLA
03  STO M          
04  X<>Y
05  STO N
06  RCL 00
07  INT
08  STO P
09  LBL 01
10  RCL 00
11  INT
12  LASTX
13  FRC
14  ENTER^
15  SIGN
16  10^X
17  CHS		 18  10^X
19  SQRT
20  +
21  ST* Y
22  +
23  ISG X
24  STO Q          
25  CLX
26  LBL 02
27  RCL Q
28  FRC
29  ISG X
30  STO a
31  SIGN
32  LBL 03
33  RCL N
34  RCL IND a		 35  -
36  RCL IND Q
37  ST- L
38  X<> L
39  CHS
40  X#0?
41  GTO 03
42  SIGN
43  STO Y
44  LBL 03
45  /
46  *
47  ISG a
48  GTO 03
49  RCL P          
50  RCL Q
51  +		 52  RDN
53  RCL IND T
54  *
55  +
56  ISG Q          
57  GTO 02
58  RCL 00
59  INT
60  STO Q
61  X<>Y
62  LBL 04
63  RCL M
64  RCL IND P
65  RCL IND Q
66  ST- Z
67  -
68  X#0?		 69  GTO 04
70  SIGN
71  STO Y
72  LBL 04
73  /
74  *
75  DSE Q          
76  GTO 04
77  ST+ O
78  DSE P
79  GTO 01
80  RCL M
81  SIGN
82  RCL N
83  RCL O
84  CLA
85  END

 
    ( 131 bytes / SIZE (n+1)(m+1) )
 
 

    STACK      INPUTS    OUTPUTS
         Y            y          y
         X            x       f(x,y)
         L            /          x

 
Example:      With the same data as above:

     x \ y |   2    3    4    6                                          |   R04   R08   R12   R16
--------------------------                              --------------------------------
       1   |   4    3    3    5                =                R01  |  R05   R09   R13   R17            respectively
       2   |   3    1    2    6                                  R02  |  R06   R10   R14   R18
       4   |   1    0    4    9                                  R03  |  R07   R11   R15   R19

-Here,  n = 3 and  m = 4  so   3.004  STO 00

-Then, to compute  f(3,5)

   5  ENTER^
   3  XEQ "BVI"  >>>>   f(3,5) ~  5.8333   ( in 31 seconds )

-And to obtain  f(1.6,2.7)

   2.7  ENTER^
   1.6     R/S        >>>>   f(1.6,2.7) ~  1.9038
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   Press, Vetterling, Teukolsky, Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4