The Museum of HP Calculators

# Least-Squares Approximation for the HP-41

This program is 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

-Given a set of N data points ( xi ; yi ) , we seek a best fitting curve y = a1 f1(x) + a2 f2(x) + ...... + an fn(x)           ( n <= N )
where  f1 ,  f2 , ...... ,  fn      are n  known functions and  a1 , a2 , ...... , an   are unknown.
-The least-squares principle minimizes the sum of vertical distances to a curve:        SUMi   [ yi - a1 f1(xi) + a2 f2(xi) + ...... + an fn(xi) ]2
leading to a n*n linear system called the "normal equations".

-5 programs are listed below:

"LR" is a classic linear regression program
"LSL" is another ( and less known ) least-squares line
"LSF3" solves the problem for three functions
"LSF"   solves the general case.
and  "LSFXY"  performs the same calculation for a surface.

Linear Regression  ( SigmaREG 00 )

-Here, the curve is a straight line y = m.x + p  given by:

m = ( Sum xi yi - n µx µy ) / ( Sum xi2 - ( Sum xi )2 / n )                r = the coefficient of determination = ( Sum xi yi - n µx µy ) / (( n-1) sx sy)
p =   µy - m.µx                                                                and       µx , µy ,  sx , sy    are the means and standard deviations of x and y.

Data Registers:

•   R00 = Summation of x values
•   R01 = Summation of x2 values
•   R02 = Summation of y values
•   R03 = Summation of y2 values
•   R04 = Summation of x.y values
•   R05 = Number of data points

 STACK INPUTS OUTPUTS T / r Z / p Y / m X x y(x) = m.x + p L / x

Notes:

-If your favorite statistics registers are, for instance, R11 thru R16 , replace R00 thru R05 by these registers in the following listing.
-Synthetic register M can be replace by any other unused data register.
-This program calculates m , p , r each time: thus, no other data register is used.

01  LBL "LR"
02  STO M
03  MEAN
04  *
05  RCL 05
06  *
07  RCL 04
08  -
09  STO Z
10  SDEV
11  RCL 05
12  DSE X
13  *
14  *
15  /
16  CHS
17  X<>Y
18  RCL 00
19  X^2
20  RCL 05
21  /
22  RCL 01
23  -
24  /
25  MEAN
26  LASTX
27  *
28  ST- Y
29  X<> L
30  RCL M
31  SIGN
32  RDN
33  ST* M
34  X<>Y
35  ST+ M
36  X<>Y
37  0
38  X<> M
39  END

( 55 bytes / SIZE 006 )

Example:  Find the least-squares linear function for the following data:

 y 2 3 3 5 5 x 1 2 3 4 5

and evaluate the linear estimates   y(2.5)  and  y(7)

a-  [SREG 00]   [CLS]        ( I mean SigmaREG 00 and CLSigma )

b-  2  ENTER^     1      S+            ( S+ = Sigma+ )               ( enter y-values first )
3  ENTER^      2     S+
3  ENTER^             S+
5  ENTER^      4     S+
5  ENTER^             S+

c- 2.5 XEQ "LR"  yields   3.2                thus,        y(2.5) = 3.2
RDN              0.8
RDN              1.2               whence     y = 0.8 x + 1.2
RDN              0.9428         -------      r = 0.9428

7     R/S   ---->    y(7) = 6.8

Another Least-Squares Line   ( SigmaREG 00 )

-The following program minimizes the sum of squares of perpendicular distances instead of vertical distances to a line y = m.x + p
this leads to a quadratic equation and m and p are determined by:

m = -A/2 + sign (r) ( 1 + A2/4 )1/2             where  A = ( sx/sy - sy/sx ) / r    ;    r = the coefficient of determination = ( Sum xi yi - n µx µy ) / (( n-1) sx sy)
p =   µy - m.µx                                                                                    and       µx , µy ,  sx , sy    are the means and standard deviations of x and y.

Data Registers:

•   R00 = Summation of x values
•   R01 = Summation of x2 values
•   R02 = Summation of y values
•   R03 = Summation of y2 values
•   R04 = Summation of x.y values
•   R05 = Number of data points

 STACK INPUTS OUTPUTS T / r Z / p Y / m X x y(x) = m.x + p L / x

Notes:

-If your favorite statistics registers are, for instance, R11 thru R16 , replace R00 thru R05 by these registers in the following listing.
-Synthetic register M can be replace by any other unused data register.
-This program calculates m , p , r each time: thus, no other data register is used.

01  LBL "LSL"
02  STO M
03  SDEV
04  *
05  RCL 05
06  DSE X
07  *
08  RCL 04
09  RCL 00
10  RCL 02
11  *
12  RCL 05
13  /
14  -
15  X<>Y
16  /
17  RDN
18  SDEV
19  /
20  ENTER^
21  1/X
22  -
23  R^
24  ST+ X
25  /
26  STO Y
27  X^2
28  1
29  +
30  SQRT
31  R^
32  SIGN
33  *
34  +
35  MEAN
36  LASTX
37  *
38  ST- Y
39  X<> L
40  RCL M
41  SIGN
42  RDN
43  ST* M
44  X<>Y
45  ST+ M
46  X<>Y
47  0
48  X<> M
49  END

( 67 bytes / SIZE 006 )

Example:  Find this second least-squares line for the same data:

 y 2 3 3 5 5 x 1 2 3 4 5

and evaluate the linear estimates   y(2.5)  and  y(7)    ( skip to step c if registers R00 thru R05 are unchanged )

a-  [SREG 00]   [CLS]        ( I mean SigmaREG 00 and CLSigma )

b-  2  ENTER^     1      S+            ( S+ = Sigma+ )               ( enter y-values first )
3  ENTER^      2     S+
3  ENTER^             S+
5  ENTER^      4     S+
5  ENTER^             S+

c- 2.5 XEQ "LSL"  yields   3.1799                thus          y(2.5) = 3.1799
RDN              0.8402
RDN              1.0794               whence     y = 0.8402 x + 1.0794
RDN              0.9428               -------      r = 0.9428

7     R/S   -->     y(7) = 6.9608

Linear Combination of 3 functions

-The principle of the least-squares line can be extended to the least-squares parabola, but we are going to solve a slightly more general problem:
-Given a set of data points ( xi ; yi ), we seek a best fitting curve y = a f(x) + b g(x) + c h(x)
where f , g , h are 3 known functions and a , b , c are unknown.

Data Registers:     •  Registers R01 R02 R03 are to be initialized and R07 thru R15 are to be cleared

•  R01 = f name     R04 = a    R07 = Sum f.f     R10 = Sum g.g    R13 = Sum y.f
•  R02 = g name    R05 = b    R08 = Sum f.g    R11 = Sum g.h    R14 = Sum y.g    ( R00 and R16: temporary data storage )
•  R03 = h name    R06 = c    R09 = Sum f.h    R12 = Sum h.h    R15 = Sum y.h

Flags: none

Subroutines:   the 3 functions f , g , h  ( assuming x is in X-register upon entry )

Note:  No coefficient of determination is calculated by this program.

001  LBL A
002  X<>Y
003  STO 00
004  X<>Y
005  STO 16
006  XEQ IND 01
007  STO 04
008  X^2
009  ST+ 07
010  LASTX
011  RCL 00
012  *
013  ST+ 13
014  RCL 16
015  XEQ IND 02
016  STO 05
017  X^2
018  ST+ 10
019  LASTX
020  RCL 00
021  *
022  ST+ 14
023  RCL 16
024  XEQ IND 03
025  STO 06
026  X^2
027  ST+ 12
028  LASTX
029  RCL 00
030  *
031  ST+ 15
032  RCL 04
033  RCL 05
034  *
035  ST+ 08
036  RCL 06
037  LASTX
038  *
039  ST+ 11
040  RCL 04
041  RCL 06
042  *
043  ST+ 09
044  RCL 16
045  RTN
046  LBL "LSF3"
047  RCL 07                        lines 47 to 137 solve the normal equations: a 3*3 linear system.
048  RCL 10
049  RCL 12
050  *
051  RCL 11
052  X^2
053  -
054  STO 04
055  *
056  RCL 09
057  RCL 11
058  *
059  RCL 08
060  RCL 12
061  *
062  -
063  STO 05
064  RCL 08
065  *
066  +
067  RCL 08
068  RCL 11
069  *
070  RCL 09
071  RCL 10
072  *
073  -
074  STO 06
075  RCL 09
076  *
077  +
078  STO 00
079  RCL 04
080  RCL 13
081  *
082  RCL 05
083  RCL 14
084  *
085  +
086  RCL 06
087  RCL 15
088  *
089  +
090  RCL 00
091  /
092  STO 04
093  RCL 05
094  RCL 13
095  *
096  RCL 07
097  RCL 12
098  *
099  RCL 09
100  X^2
101  -
102  RCL 14
103  *
104  +
105  RCL 08
106  RCL 09
107  *
108  RCL 07
109  RCL 11
110  *
111  -
112  STO 16
113  RCL 15
114  *
115  +
116  RCL 00
117  /
118  STO 05
119  RCL 06
120  RCL 13
121  *
122  RCL 14
123  RCL 16
124  *
125  +
126  RCL 07
127  RCL 10
128  *
129  RCL 08
130  X^2
131  -
132  RCL 15
133  *
134  +
135  RCL 00
136  /
137  STO 06
138  RCL 05
139  RCL 04
140  RTN
141  LBL F
142  STO 00
143  XEQ IND 01
144  RCL 04
145  *
146  STO 16
147  RCL 00
148  XEQ IND 02
149  RCL 05
150  *
151  ST+ 16
152  RCL 00
153  XEQ IND 03
154  RCL 06
155  *
156  ST+ 16
157  RCL 16
158  RTN
159  GTO F
160  END

( 198 bytes / SIZE 017 )

Example:  Find 3 real numbers   a ; b ; c  such that  y = a + b x2 + c / ln x   best fits the following data:

 y 12.9141 0.4989 -18.2283 -40.5506 -68.2507 x 1.5 2 3 4 5

Then, estimate y(6) and y(7)

-Let's clear registers R07 thru R15:                        7.015  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-We load these 3 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  1
03  RTN
04  LBL "U"        ( any global label, maximum of 6 characters )
05  X^2
06  RTN
07  LBL "V"        ( any global label, maximum of 6 characters )
08  LN
09  1/X
10  RTN

-We store their names into registers R01 R02 R03     "T" ASTO 01  "U"  ASTO 02  "V"  ASTO 03
-and we store the data ( GTO "LSF3" )

( enter y-values first )

12.9141   ENTER^    1.5   [A]      ( the Sigma+ key in user mode )          x = 1.5  ends up in X-register
0.4989   ENTER^      2    [A]                                                                x =   2   ends up in X-register
-18.2283   ENTER^     3    [A]
-40.5506   ENTER^     4    [A]                                                                              ....  etc  ....
-68.2507   ENTER^     5    [A]

-Then  R/S  produces   2.4001               a =   2.4001 = R04
RDN               -3.0000               b = -3.0000 = R05
RDN                6.9999                c =  6.9999 = R06

Thus,    y = 2.4001 - 3.0000 x2 + 6.9999 / ln x

-To obtain y(6)       6   [F]  or   R/S   yields   -101.6933             6  is saved in R00
--------- y(7)        7   [F]  or   R/S   -----    -141.0029             7 ---------------

Linear Combination of n functions

Data Registers:     •  Registers    R00 R01 R02  ......  Rnn   are to be initialized and  R2n+1 thru Rn2+3n  are to be cleared

•  R00 = n
•  R01 = f1 name    Rnn+1 = a1    R2n+1 = Sum f1.f1  .....................  Rn2+n+1 = Sum fn.f1    Rn2+2n+1  = Sum y.f1
•  R02 = f2 name    Rnn+2 = a2    R2n+2 = Sum f1.f2  .....................  Rn2+n+2  = Sum fn.f2   Rn2+2n+2  = Sum y.f2
...................................................................................................................................

•  Rnn = fn name      R2n = an        R3n   = Sum f1.fn  .....................     Rn2+2n  = Sum fn.fn     Rn2+3n   = Sum y.fn

Flags:  F01 is cleared by this program

Subroutines:   "LS"  synthetic version    ( cf linear and non-linear systems for the HP-41 )
and the n functions  f1, f2, ... ,  fn  ( assuming x is in X-register upon entry )

Note:  No coefficient of determination is calculated by this program.

001  LBL A
002  STO M
003  X<>Y
004  STO N
005  RCL 00
006  STO O
007  LBL 00
008  RCL IND O
009  RCL M
010  XEQ IND Y
011  RCL O
012  RCL 00
013  +
014  X<>Y
015  STO IND Y
016  LASTX
017  X^2
018  LASTX
019  +
020  ST+ Z
021  CLX
022  RCL N
023  *
024  ST+ IND Y
025  DSE O
026  GTO 00
027  RCL 00
028  X^2
029  LASTX
030  ST+ X
031  ST+ Y
032  RCL 00
033   E3
034  /
035  +
036  STO O
037  LBL 01
038  RCL IND X
039  RCL IND O
040  *
041  ST+ IND Z
042  RDN
043  DSE Y
044  DSE X
045  GTO 01
046  RCL 00
047  +
048  DSE O
049  GTO 01
050  RCL M
051  RTN
052  LBL "LSF"
053  3
054  RCL 00
055  ST+ Y
056  ST* Y
057   E2
058  /
059  +
060   E3
061  /
062  RCL 00
063  ST+ X
064  +
065  1
066  +
067  0
068  CF 01
069  XEQ "LS"
070  R^
071  STO 00
072  1
073  +
074  X^2
075  RCL 00                     If you don't have an X-function module, replace lines 76 to 91 with:
076  .1                               ST+ X
077  %                                E3
078  +                                /
079  1                                RCL 00
080  +                                +
081  STO Z                        1
082   E3                             +
083  /                                 STO Z
084  +                                SIGN
085  REGMOVE               LBL 03
086  CLX                          CLX
087  RCL 00                     RCL IND Y
088   E3                            STO IND Z
089  /                                 ISG Y
090  +                               CLX
091  RTN                          ISG Z
092  LBL F                       GTO 03
093  CLA                          LASTX
094  STO M                      RTN
095  RCL 00
096  ST+ 00
097  STO N
098  LBL 02
099  RCL IND N
100  RCL M
101  XEQ IND Y
102  RCL IND 00
103  *
104  ST+ O
105  DSE 00
106  DSE N
107  GTO 02
108  X<> M
109  SIGN
110  X<> O
111  CLA
112  RTN
113  GTO F
114  END

( 178 bytes / SIZE n2+ 3n +1 )

Example:  Find 4 real numbers a1 , a2 , a3 , a4  such that    y = a1 + a2 sin 10.x + a3.x + a4 ln x   best fits the following data ( x in degrees ):

 y 5.22 6.69 7.28 7.01 3.06 -0.14 x 2 3 4 6 10 12

Then, evaluate y(5) and y(7)

-SIZE 029  or greater.
-Set your HP-41 in DEG mode.
-Clear registers R09 thru R28:                        9.028  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-Load the 4 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  1
03  RTN
04  LBL "U"        ( any global label, maximum of 6 characters )
05  10
06  *
07  SIN
08  RTN
09  LBL "V"        ( any global label, maximum of 6 characters )
10  RTN
11  LBL "W"       ( any global label, maximum of 6 characters )
12  LN
13  RTN

-There are 4 functions:      4  STO 00
-Store their names into registers R01 thru R04     "T" ASTO 01    "U"  ASTO 02    "V"  ASTO 03   "W"  ASTO 04
-and store the data ( GTO "LSF" )

( enter y-values first )

5.22    ENTER^    2     [A]      ( the Sigma+ key in user mode )      x = 2  ends up in X-register within  11 seconds
6.69    ENTER^    3     [A]                                                            x = 3   ----------------------------------------
7.28    ENTER^    4     [A]
7.01    ENTER^    6     [A]                                                                                 ......  etc  ......
3.06    ENTER^   10    [A]
-0.14    ENTER^   12     [A]

-Then  R/S  or  XEQ "LSF"  produces ( in 31 seconds )        5.008  =  the control number of the solution:
a1 , a2 , a3 , a4  are in registers   R05 , R06 , R07 , R08
RCL 05  gives       a1    =  2.9957
RCL 06  -----       a2    =  4.0011
RCL 07  -----       a3    = -2.0014
RCL 08  -----       a4    =  7.0089

Thus,    y = 2.9957 + 4.0011 sin 10.x - 2.0014 x + 7.0089 ln x    is the best fitting curve for these data.

-To obtain y(5)       5   [F]  or   R/S   yields   7.3339     ( in 3.5 seconds )              x = 5  is saved in L-register
--------- y(7)        7   [F]  or   R/S   -----    6.3841                                             x = 7  --------------------

Remark:   Execution time is approximately proportional to n2

Linear Combination of n functions ( 2 dimensional problem )

-Given a set of N data points ( xi ; yi ; zi ) in space, we seek a best fitting surface z = a1 f1(x,y) + a2 f2(x,y) + ...... + an fn(x,y)           ( n <= N )
where  f1 ,  f2 , ...... ,  fn      are n  known functions of 2 variables and  a1 , a2 , ...... , an   are unknown.
-The following program is quite similar to "LSF"

Data Registers:     •  Registers    R00 R01 R02  ......  Rnn   are to be initialized and  R2n+1 thru Rn2+3n  are to be cleared

•  R00 = n
•  R01 = f1 name    Rnn+1 = a1    R2n+1 = Sum f1.f1  .....................  Rn2+n+1 = Sum fn.f1    Rn2+2n+1  = Sum z.f1
•  R02 = f2 name    Rnn+2 = a2    R2n+2 = Sum f1.f2  .....................  Rn2+n+2  = Sum fn.f2   Rn2+2n+2  = Sum z.f2
...................................................................................................................................

•  Rnn = fn name      R2n = an        R3n   = Sum f1.fn  .....................     Rn2+2n  = Sum fn.fn     Rn2+3n   = Sum z.fn

Flags:  F01 is cleared by this program

Subroutines:   "LS"  synthetic version        ( cf linear and non-linear systems for the HP-41 )
and the n functions  f1, f2, ... ,  fn  ( assuming x is in X-register and y is in Y-register upon entry )

Notes:  -No coefficient of determination is calculated by this program.
-Status register P is used in LBL A and LBL F, therefore the fi suroutines must have no digit entry line,
and unfortunately, register a wouldn't work because XEQ IND Z instructions clear this register.

001  LBL A
002  STO M
003  RDN
004  STO N
005  X<>Y
006  STO O
007  RCL 00
008  STO P
009  LBL 00
010  RCL IND P
011  RCL N
012  RCL M
013  XEQ IND Z
014  RCL P
015  RCL 00
016  +
017  X<>Y
018  STO IND Y
019  LASTX
020  X^2
021  LASTX
022  +
023  ST+ Z
024  CLX
025  RCL O
026  *
027  ST+ IND Y
028  DSE P
029  GTO 00
030  RCL 00
031  X^2
032  LASTX
033  ST+ X
034  ST+ Y
035  RCL 00
036   E3
037  /
038  +
039  STO O
040  LBL 01
041  RCL IND X
042  RCL IND O
043  *
044  ST+ IND Z
045  RDN
046  DSE Y
047  DSE X
048  GTO 01
049  RCL 00
050  +
051  DSE O
052  GTO 01
053  RCL N
054  RCL M
055  RTN
056  LBL "LSFXY"
057  3
058  RCL 00
059  ST+ Y
060  ST* Y
061   E2
062  /
063  +
064   E3
065  /
066  RCL 00
067  ST+ X
068  +
069  1
070  +
071  0
072  CF 01
073  XEQ "LS"
074  R^
075  STO 00
076  1
077  +
078  X^2
079  RCL 00                     If you don't have an X-function module, replace lines 80 to 95 with:
080  .1                               ST+ X
081  %                                E3
082  +                                /
083  1                                RCL 00
084  +                                +
085  STO Z                        1
086   E3                             +
087  /                                 STO Z
088  +                                SIGN
089  REGMOVE               LBL 03
090  CLX                          CLX
091  RCL 00                     RCL IND Y
092   E3                            STO IND Z
093  /                                 ISG Y
094  +                               CLX
095  RTN                          ISG Z
096  LBL F                       GTO 03
097  CLA                          LASTX
098  STO M                      RTN
099  X<>Y
100  STO N
101  RCL 00
102  ST+ 00
103  STO P
104  LBL 02
105  RCL IND P
106  RCL N
107  RCL M
108  XEQ IND Z
109  RCL IND 00
110  *
111  ST+ O
112  DSE 00
113  DSE P
114  GTO 02
115  RCL N
116  RCL M
117  SIGN
118  X<> O
119  CLA
120  RTN
121  GTO F
122  END

( 195 bytes / SIZE n2+ 3n +1 )

Example1:  Find 3 real numbers   a ; b ; c  such that the plane y = a x + b y + c   best fits the following data:

 z 1.1 0.1 -7.2 9.75 4.97 y 2 1 0 3 2 x 1 2 0 4 3

Then, evaluate z(1;3) and z(-2;4)

-SIZE 019  or greater.
-Clear registers R07 thru R18:                        7.018  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-Load the 3 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  RTN
03  LBL "U"        ( any global label, maximum of 6 characters )
04  X<>Y
05  RTN
06  LBL "V"        ( any global label, maximum of 6 characters )
07  CLX
08  SIGN            ( 1 would alter register P )
09  RTN

-There are 3 functions:      3  STO 00
-Store their names into registers R01 thru R03     "T" ASTO 01    "U"  ASTO 02    "V"  ASTO 03
-and store the data ( GTO "LSFXY" )

( enter z-values first , then y-values and x-values at last )

1.1    ENTER^    2    ENTER^    1    [A]      ( the Sigma+ key in user mode )      y = 2 and x = 1 end up in Y and X-registers
0.1    ENTER^    1    ENTER^    2    [A]                                                            y = 1 and x = 2 ----------------------------
-7.2   ENTER^    0    ENTER^           [A]
9.75  ENTER^    3    ENTER^    4    [A]                                                                                 ......  etc  ......
4.97  ENTER^    2    ENTER^    3    [A]

-Then  R/S  or  XEQ "LSFXY"  yields        4.006  =  the control number of the solution:
a , b , c   are in registers   R04 , R05 , R06:
RCL 04  gives       a    =  1.9485
RCL 05  -----       b    =  3.0475
RCL 06  -----       c    = -7.0290

Thus, the least-squares plane is    z = 1.9485 x + 3.0475 y - 7.0290

3  ENTER^   1    R/S   ( or [F]  )   --------->    z(1;3) = 4.0620                 ( y = 3 ends up in Y-register and x =  1 ends up in L-register )
4  ENTER^  -2   R/S    ( or [F] )   --------->    z(-2;4) = 1.2640                ( y = 4 ends up in Y-register and x = -2 ends up in L-register )

Example2:  Find 4 real numbers   a1 , a2 , a3 , a4  such that    y = a1 sin ( x + y ) + a2 ex / y + a3.x.y + a4 ln ( x.y )   best fits the following data ( x , y in radians ):

 z -0.647 -6.301 -22.679 -57.46 -2.336 y 2 2 3 1 1 x 1 2 4 3 1

Then, evaluate z(1;4)

-SIZE 029  or greater.
-Clear registers R09 thru R28:                        9.028  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-Load the 4 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  +
03  SIN
04  RTN
05  LBL "U"        ( any global label, maximum of 6 characters )
06  E^X
07  X<>Y
08  /
09  RTN
10  LBL "V"        ( any global label, maximum of 6 characters )
11  *
12  RTN
13  LBL "W"       ( any global label, maximum of 6 characters )
14  *
15  LN
16  RTN

-There are 4 functions:      4  STO 00
-Store their names into registers R01 thru R04     "T" ASTO 01    "U"  ASTO 02    "V"  ASTO 03    "W"   ASTO 04
-and store the data ( GTO "LSFXY" )

( enter z first , then y and x at last )

-0.647    ENTER^    2    ENTER^    1    [A]      ( the Sigma+ key in user mode )      y = 2 and x = 1 end up in Y and X-registers
-6.301    ENTER^    2    ENTER^          [A]                                                            y = 2 and x = 2 ----------------------------
-22.679   ENTER^    3    ENTER^    4    [A]
-57.460   ENTER^    1    ENTER^    3    [A]                                                                                 ......  etc  ......
-2.336    ENTER^    1    ENTER^          [A]

-Then  R/S  or  XEQ "LSFXY"    yields        5.008  =  the control number of the solution:         a1 , a2 , a3 , a4  are in registers   R05 , R06 , R07 , R08

RCL 05  gives       a1    =  2.0005
RCL 06  -----       a2    = -3.0000
RCL 07  -----       a3    =  3.9996
RCL 08  -----       a4    = -6.9985

Thus, the least-squares surface is    z = 2.0005 sin ( x + y ) - 3.0000 ex / y + 3.9996 x.y - 6.9985 ln ( x.y )

Then:   4  ENTER^ 1    R/S   ( or [F]  )   --------->    z(1;4) = 2.3393         ( y = 4 ends up in Y-register and x = 1 is saved in L-register )
... and so on ...

Remark:   If one function is, for instance,   sin [12.345 ( 2x + y )] , store 12.345  into an unused regiter ( say R99 ) and define this function by

ST+ X
+
RCL 99
*
SIN
RTN

Reference:   F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN  07-055197-9