Lagrangian Interpolation

[Thomas Klemm]

This program uses the method mentioned in my previous post.

Memory Map
\(
\begin{matrix}
I & : & x \\
R_0 & : & x_1 \\
R_1 & : & x_2 \\
R_2 & : & x_3 \\
R_3 & : & y_1 \\
R_4 & : & y_2 \\
R_5 & : & y_3 \\
R_6 & : & w_1 \\
R_7 & : & w_2 \\
R_8 & : & w_3 \\
R_{S4} & : & \Sigma x \\
R_{S8} & : & \Sigma xy \\
\end{matrix}
\)

Code:

001  31 25 11  : LBL A
002        01  : 1
003        84  : R/S
004     33 00  : STO 0
005     35 52  : x<>y
006     33 03  : STO 3
007        02  : 2
008        84  : R/S
009     33 01  : STO 1
010     35 52  : x<>y
011     33 04  : STO 4
012        03  : 3
013        84  : R/S
014     33 02  : STO 2
015     35 52  : x<>Y
016     33 05  : STO 5
017        01  : 1
018        42  : CHS
019     33 06  : STO 6
020     33 07  : STO 7
021     33 08  : STO 8
022     34 00  : RCL 0
023     34 01  : RCL 1
024        51  : -
025  33 81 06  : STO/ 6
026  33 81 07  : STO/ 7
027     34 01  : RCL 1
028     34 02  : RCL 2
029        51  : -
030  33 81 07  : STO/ 7
031  33 81 08  : STO/ 8
032     34 02  : RCL 2
033     35 00  : RCL 0
034        51  : -
035  33 81 08  : STO/ 8
036  33 81 06  : STO/ 6
037     31 42  : P<>S
038        84  : R/S
039  31 25 12  : LBL B
040     31 43  : CL REG
041     35 33  : ST I
042     31 42  : P<>S
043     34 03  : RCL 3
044     34 06  : RCL 6
045     35 34  : RC I
046     34 00  : RCL 0
047        51  : -
048        81  : /
049        21  : Σ+
050     34 04  : RCL 4
051     34 07  : RCL 7
052     35 34  : RC I
053     34 01  : RCL 1
054        51  : -
055        81  : /
056        21  : Σ+
057     34 05  : RCL 5
058     34 08  : RCL 8
059     35 34  : RC I
060     34 02  : RCL 2
061        51  : -
062        81  : /
063        21  : Σ+
064     31 42  : P<>S
065     34 08  : RCL 8
066     34 04  : RCL 4
067        81  : /
068     35 22  : RTN

Explanation

Lines 001-016 are the same as in Namir's orginal program.

Lines 017-036 calculate the weights

\(w_j = \frac{1}{\prod_{i=0,i \neq j}^k(x_j-x_i)}\).

However this formula is slightly changed so that we only need the differences \(x_1-x_2\), \(x_2-x_3\) and \(x_3-x_1\):

\(w_1=\frac{-1}{(x_1-x_2)(x_3-x_1)}\)

\(w_2=\frac{-1}{(x_2-x_3)(x_1-x_2)}\)

\(w_3=\frac{-1}{(x_3-x_1)(x_2-x_3)}\)

As you can see the same factor can be used in two of the weights. Nothing fancy it's just to get the sign correct. That's why we start with \(-1\) in lines 017-021.

Lines 043-063 recall \(y_j\) and calculate \(\frac{w_j}{x-x_j}\) and then use \(\Sigma +\) to calculate both \(\sum_{j=0}^k \frac{w_j}{x-x_j}\) and \(\sum_{j=0}^k \frac{w_j}{x-x_j}y_j\) in one single step.

Lines 065-067 finally recall these values and calculate \(\frac{\Sigma xy}{\Sigma x}\).

Cheers
Thomas


For those with a HP-15C here's the corresponding program:

Code:

001 - 42,21,11  LBL A
002 -        1  1
003 -       31  R/S
004 -    44  8  STO 8
005 -       34  x<>y
006 -    44 .1  STO .1
007 -        2  2
008 -       31  R/S
009 -    44  9  STO 9
010 -       34  x<>y
011 -    44 .2  STO .2
012 -        3  3
013 -       31  R/S
014 -    44 .0  STO .0
015 -       34  x<>y
016 -    44 .3  STO .3
017 -        1  1
018 -       16  CHS
019 -    44 .4  STO .4
020 -    44 .5  STO .5
021 -    44 .6  STO .6
022 -    45  8  RCL 8
023 - 45,30, 9  RCL- 9
024 - 44,10,.4  STO/ .4
025 - 44,10,.5  STO/ .5
026 -    45  9  RCL 9
027 - 45,30,.0  RCL- .0
028 - 44,10,.5  STO/ .5
029 - 44,10,.6  STO/ .6
030 -    45 .0  RCL .0
031 - 45,30, 8  RCL- 8
032 - 44,10,.6  STO/ .6
033 - 44,10,.4  STO/ .4
034 -       31  R/S
035 - 42,21,12  LBL B
036 -    44  0  STO 0
037 -    42 32  CLEAR Σ
038 -    45 .1  RCL .1
039 -    45 .4  RCL .4
040 -    45  0  RCL 0
041 - 45,30, 8  RCL- 8
042 -       10  /
043 -       49  Σ+
044 -    45 .2  RCL .2
045 -    45 .5  RCL .5
046 -    45  0  RCL 0
047 - 45,30, 9  RCL- 9
048 -       10  /
049 -       49  Σ+
050 -    45 .3  RCL .3
051 -    45 .6  RCL .6
052 -    45  0  RCL 0
053 - 45,30,.0  RCL- .0
054 -       10  /
055 -       49  Σ+
056 -    45  7  RCL 7
057 - 45,10, 3  RCL/ 3
058 -    43 32  RTN

Added a card of the program that can be used with Jacques Laporte's HP-67 Microcode Simulator.