Simpsonregel

program Simpson(input, output);

    {Simpson-Integration einer Funktion mit n + 1 gegebenen Stuetzstellen}
    
    uses
       stdfile;
       
    const
        MaxArray = 1024;
        
    var
        fx:         array [0..MaxArray] of double;
        a, b:       double;
        m, n, i:    integer;
        h, t:       double;
        res:        double;
        sum1, sum2: double;
        
    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
        	writeln('Simpson-Integration einer Funktion f');
        	writeln;
        	writeln('Welche Datei enthaelt Intervallgrenzen und Stuetzstellen? ');
                writeln('(In simpson.dat liegt ein Beispiel.)');
        	readln(fName);
        	if StdOpen(f, fName) then begin
        	    if not StdRead(f, a) then begin
                    writeln('Erste Intervallgrenze nicht gefunden.');
                    halt
                end; {if}
                if not StdRead(f, b) then begin
                    writeln('Zweite Intervallgrenze nicht gefunden.');
                    halt
                end; {if}
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                        writeln('Zu viele Stuetzstellen (mehr als ', MaxArray,
                            ') vorhanden.');
                        halt
                    end {if}
                until not StdRead(f, fx[n]);
                n := n - 1;
                close(f);
                writeln('Insgesamt ', n + 1, ' Werte gelesen.');
                writeln
            end else begin {if}
                writeln('Kann Datei "', fName, '" nicht oeffnen.');
                halt
            end {if}
         end; {init}

    procedure result;
        begin {result}
            writeln('Das Integral hat den Wert ', res:7:3)
        end; {result}

    begin {Simpson}
        init;
        h := (b - a) / n;
        if (n mod 2 <> 0) and (n > 1) then begin
            t := fx[n - 3] + 3 * (fx[n - 2] + fx[n - 1]) + fx[n];
            res := 3 * h * t / 8;
            m := n - 3;
        end else begin {if}
            res := 0;
            m := n
        end; {if}
        if m > 1 then begin
            sum1 := 0;
            sum2 := 0;
            i := 1;
            while i <= m - 1 do begin
                sum1 := sum1 + fx[i];
                i := i + 2
            end; {while}
            i := 2;
            while i <= m - 2 do begin
                sum2 := sum2 + fx[i];
                i := i + 2
            end; {while}
            res := res + h * (fx[0] + 4 * sum1 + 2 * sum2 + fx[m]) / 3
        end; {if}
        result     
    end. {Simpson}