两个算法是等效的：
{$IFDEF M1}
const TIMES = 2;
procedure Xiao(const X: array of Integer; const Y: array of Double; N: Integer;
var d: array of Double);
var
dX1, dX2, dX3, dX4, dY, dYX, dYX2, dYX3: Double;
dX, Ax, Bx: Double;
i: Integer;
begin
dX1 := 0.0;
dX2 := 0.0;
dX3 := 0.0;
dX4 := 0.0;
dY := 0.0;
dYX := 0.0;
dYX2 := 0.0;
dYX3 := 0.0;
for i := 0 to N - 1 do
begin
dX := X;
dX1 := dX1 + dX;
dX2 := dX2 + dX*dX;
dX3 := dX3 + dX*dX*dX;
dX4 := dX4 + dX*dX*dX*dX;
dY := dY + Y;
dYX := dYX + dX * Y;
dYX2 := dYX2 + dX*dX * Y;
dYX3 := dYX3 + dX*dX*dX * Y;
end;
dX1 := dX1 / N;
dX2 := dX2 / N;
dX3 := dX3 / N;
dX4 := dX4 / N;
dY := dY / N;
dYX := dYX / N;
dYX2 := dYX2 / N;
//dYX3 := dYX3 / N;
Ax := (dX4 - dX2 * dX2) - (dX3 - dX1 * dX2) / (dX2 - dX1 * dX1)*(dX3 - dX1 * dX2);
Bx := (dYX2 - dY*dX2) - (dX3 - dX1 * dX2) / (dX2 - dX1 * dX1) * (dYX - dY * dX1);
d[2] := Bx / Ax;
d[1] := ((dYX - dY * dX1) - (dX3 - dX1 * dX2) * d[2]) / (dX2 - dX1 * dX1);
d[0] := dY - d[2] * dX2 - d[1] * dX1;
end;
{$ENDIF}
{$IFDEF M2}
const TIMES = 3;
procedure Xiao(const X: array of Integer; const Y: array of Double; N: Integer;
var d: array of Double);
var
a1, a2, a3, b1, b2, c0, c1, c2, c3: Double;
P1, P12, EP12, EXP12, EYP1: Double;
P2, P22, EP22, EXP22, EYP2: Double;
P3, P32, EP32, EYP3: Double;
i: Integer;
begin
a1 := 0;
c0 := 0;
for i := 0 to N - 1 do
begin
a1 := a1 + X;
c0 := c0 + Y;
end;
a1 := a1 / N;
c0 := c0 / N;
EP12 := 0;
EXP12 := 0;
EYP1 := 0;
for i := 0 to N - 1 do
begin
P1 := X - a1;
P12 := P1 * P1;
EP12 := EP12 + P12;
EXP12 := EXP12 + X * P12;
EYP1 := EYP1 + Y * P1;
end;
a2 := EXP12 / EP12;
b1 := EP12 / N;
c1 := EYP1 / EP12;
EP22 := 0;
EXP22 := 0;
EYP2 := 0;
for i := 0 to N - 1 do
begin
P1 := X - a1;
P2 := (X - a2) * P1 - b1;
P22 := P2 * P2;
EP22 := EP22 + P22;
EXP22 := EXP22 + X * P22;
EYP2 := EYP2 + Y * P2;
end;
a3 := EXP22 / EP22;
b2 := EP22 / EP12;
c2 := EYP2 / EP22;
EP32 := 0;
EYP3 := 0;
for i := 0 to N - 1 do
begin
P1 := X - a1;
P2 := (X - a2) * P1 - b1;
P3 := (X - a3) * P2 - b2 * P1;
P32 := P3 * P3;
EP32 := EP32 + P32;
EYP3 := EYP3 + Y * P3;
end;
c3 := EYP3 / EP32;
d[0] := c0 - c1 * a1 + c2 * (a1 * a2 - b1);
d[1] := c1 - c2 * (a1 + a2);
d[2] := c2;
if TIMES = 3 then
begin
d[0] := d[0] - c3 * (a3 * (a1 * a2 - b1) - b2 * a1);
d[1] := d[1] + c3 * (a1 * a2 - b1 + a3 * (a1 + a2) - b2);
d[2] := d[2] - c3 * (a1 + a2 + a3);
d[3] := c3;
end;
end;
{$ENDIF}
