unit2.pas

Numeric Programs

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
procedure TOEPLZ(R:array of double;var X:array of real; Y:array of real; N:integer);
implementation
procedure TOEPLZ(R:array of double;var X:array of real;
                                   Y:array of real; N:integer);
label 99;
var
    G, H:array[0..100] of double;
    M,M1,J,K,M2:integer;
    SXN,SD,SGN,SGD,SHN,PT1,PT2,QT1,QT2,PP,QQ:double;
begin
    If R[N] = 0 Then GoTo 99;
    X[1]:=Y[1] / R[N];
    If N = 1 Then Exit;
    G[1]:=R[N - 1] / R[N];
    H[1]:=R[N + 1] / R[N];
    For M:=1 To N do
    begin
        M1:=M + 1;
        SXN:=-Y[M1];
        SD:=-R[N];
        For J:=1 To M do
        begin
            SXN:=SXN + R[N + M1 - J] * X[J];
            SD:=SD + R[N + M1 - J] * G[M - J + 1];
        end;
        If SD = 0 Then GoTo 99;
        X[M1]:=SXN / SD;
        For J:=1 To M do
            X[J]:=X[J] - X[M1] * G[M - J + 1];
        If M1 = N Then Exit;
        SGN:=-R[N - M1];
        SHN:=-R[N + M1];
        SGD:=-R[N];
        For J:=1 To M do
        begin
            SGN:=SGN + R[N + J - M1] * G[J];
            SHN:=SHN + R[N + M1 - J] * H[J];
            SGD:=SGD + R[N + J - M1] * H[M - J + 1];
        end;
        If (SD = 0) Or (SGD = 0) Then GoTo 99;
        G[M1]:=SGN / SGD;
        H[M1]:=SHN / SD;
        K:=M;
        M2:=(M + 1) div 2;
        PP:=G[M1];
        QQ:=H[M1];
        For J:=1 To M2 do
        begin
            PT1:=G[J];
            PT2:=G[K];
            QT1:=H[J];
            QT2:=H[K];
            G[J]:=PT1 - PP * QT2;
            G[K]:=PT2 - PP * QT1;
            H[J]:=QT1 - QQ * PT2;
            H[K]:=QT2 - QQ * PT1;
            K:=K - 1;
        end;
    end;
    ShowMessage('Never get here');
99: ShowMessage('Levinson method fails: singular principal minor');
end;
end.