Polynomiale Regression

program PolynomialeRegression(input, output);

        {Programmsequenz zur polynomialen Regression von n Datenpunkten}

    uses
        stdfile;

    const
        MaxArray = 32; {<= 32767}

    var
        i, j, k, l, n: integer;
        order:         integer;
        x, y, c:       array [0..MaxArray] of double;
        a:             array [0..MaxArray,0..MaxArray] of double;
        sum:           double;

    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Berechnung der Matrix A und des Vektors c, mit Ar = c, sodass die r_i');
            writeln('die Koeffizienten der polynomialen Regression von n Datenpunkten sind.');
            writeln;
            writeln('In welcher Datei liegt die Funktion? (In pol_reg.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                        writeln('Die Datei "', fName, '" hat zu viele Punkte (mehr als ', MaxArray, ').');
                        halt
                    end; {if}
                until not StdRead2(f, x[n], y[n]);
                n := n - 1
            end else begin {if}
                writeln('Die Datei "', fName, '" kann nicht geoeffnet werden.');
                halt
            end; {if}
            close(f);
            writeln('Insgesamt wurden ', n + 1, ' Werte eingelesen.');
            n := n + 1;
            writeln('Geben Sie die Ordnung des Polynoms an');
            readln(order);
        end; {init}

    procedure result; begin
        for i:=1 to order+1 do begin
            for j:=1 to order+1 do begin
                write(a[i, j]:8:3, ' ')
            end; {for}
            writeln(':', c[i]:8:3)
        end {for}
    end; {result}

    begin {Polynomiale Regression}
        init;
        for i:=1 to order+1 do begin
            for j:=1 to i do begin
                k := i + j - 2;
                sum := 0.0;
                for l:=0 to n-1 do begin
                    sum := sum + exp(k*ln(x[l]))
                end; {for}
                a [i, j] := sum;
                a [j, i] := sum
            end; {for}
            sum := 0.0;
            for l:=0 to n-1 do begin
                sum := sum + y[l]*exp((i - 1)*ln(x[l]))
            end; {for}
            c[i] := sum
        end; {for}
        result
    end. {Polynomiale Regression}