Berechnung von Lagrange-Polynomen

program Lagrange(input, output);

    {Berechnet ein Lagrange-Polynom n-ten Grades anhand von n + 1 
    vorgegebenen Stuetzstellen.}

    uses
        stdfile;

    const
        MaxArray = 1024; {<= 32767}

    var
        x, fx:   array [0..MaxArray] of double;
        n:       integer;
        x0:       double;   {an dieser Stelle wird das Polynom berechnet}

        y:       double;   {enthaelt das Ergebnis}
        produkt: double;   {Zwischenergebnis: der i-te Produktterm}
        i, j:    integer;  {Zaehler}


    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
             writeln('Berechnet ein Lagrange-Polynom n-ten Grades');
             writeln;
             write('In welcher Datei liegen die Stuetzstellen? (In lagrange.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 Dimension ueberschreitet die vorgegebene Arraygroesse.');
                         halt
                     end; {if}
                 until not StdRead2(f, x[n], fx[n]);
                 n := n - 1;
             end else begin {if}
                 writeln('Kann Datei "', fName, '" nicht oeffnen.');
                 halt
             end; {if}
             close(f);
             writeln('Insgesamt ', n + 1, ' Stuetzstellen gelesen.');
             writeln;
             write('An welcher Stelle soll das Polynom berechnet werden? ');
             readln(x0)
         end; {init}

    procedure result;
        begin {result}
            writeln('Das Lagrange-Polynom an der Stelle ', x0:6:2,
                ' hat den Wert ', y:6:2, '.')
        end; {result}

    begin {Lagrange}
        init;
        y := 0;        {Ergebnis initialisieren}
        for i := 0 to n do begin
            produkt := fx[i];
            for j := 0 to n do begin
                if i <> j then begin
                    produkt := produkt * (x0 - x[j]) /
                       (x[i] - x[j])
                end {if}
            end; {for}
            y := y + produkt         
        end; {for}
        result     
    end. {Lagrange}