proc PolyReg(XA, YA, N);
real XA, YA; int N;
int  I;
real XM, YM, X2M, X3M, X4M, XYM, X2YM;
real SXX, SXY, SXX2, SX2X2, SX2Y;
real A, B, C, Y, RN;
[XM:= 0.;  YM:= 0.;
X2M:= 0.;  X3M:= 0.;  X4M:= 0.;
XYM:= 0.;  X2YM:= 0.;
for I:= 0 to N-1 do
    [XM:= XM + XA(I);
     YM:= YM + YA(I);
     X2M:= X2M + XA(I)*XA(I);
     X3M:= X3M + XA(I)*XA(I)*XA(I);
     X4M:= X4M + XA(I)*XA(I)*XA(I)*XA(I);
     XYM:= XYM + XA(I)*YA(I);
     X2YM:= X2YM + XA(I)*XA(I)*YA(I);
    ];
RN:=  float(N);
XM:=  XM/RN;
YM:=  YM/RN;
X2M:= X2M/RN;
X3M:= X3M/RN;
X4M:= X4M/RN;
XYM:= XYM/RN;
X2YM:= X2YM/RN;

SXX:=   X2M - XM*XM;
SXY:=   XYM - XM*YM;
SXX2:=  X3M - XM*X2M;
SX2X2:= X4M - X2M*X2M;
SX2Y:= X2YM - X2M*YM;

B:= (SXY*SX2X2 - SX2Y*SXX2) / (SXX*SX2X2 - SXX2*SXX2);
C:= (SX2Y*SXX - SXY*SXX2) / (SXX*SX2X2 - SXX2*SXX2);
A:= YM - B*XM - C*X2M;

Format(1, 0);
Text(0, "y = ");  RlOut(0, A);
Text(0, " + ");   RlOut(0, B);
Text(0, "x + ");  RlOut(0, C);
Text(0, "x^^2");   CrLf(0);
Text(0, "  Input   Approx^m^j");
Text(0, "  x    y     y1^m^j");
for I:=0 to N-1 do
    [Format(3, 0);
    RlOut(0, XA(I));  Text(0, "  ");
    RlOut(0, YA(I));  Text(0, "  ");
    Y:= A + (B + C*XA(I)) * XA(I);
    Format(3, 1);
    RlOut(0, Y);
    CrLf(0);
    ];
];

real XA, YA;
[XA:= [0., 1.,  2.,  3.,  4.,  5.,   6.,   7.,   8.,   9.,  10.];
 YA:= [1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.];
PolyReg(XA, YA, 11);
]