{$N+}PROGRAM I2D;

CONST
	pi=3.141592654;
        k=80;

TYPE
        pole2d=ARRAY[1..k,1..k] OF double;
        pole1d=ARRAY[1..k] OF double;

FUNCTION f(x,y:double):double;
BEGIN
	f:=abs(sin(sqr(pi)*(sqr(x-0.5)+sqr(y-0.5))));
END;

VAR
        y:pole2d;
        z:pole1d;
        a,b,u,v,integral,krok:double;
        i,j:integer;
BEGIN
        a:=0;
        b:=1;
        krok:=(b-a)/(k-1);
        integral:=0;


        FOR j:=1 TO k DO
        BEGIN
                FOR i:=1 TO k DO
                BEGIN
                        u:=a+(i-1)*krok;
                        v:=a+(j-1)*krok;
                        y[j,i]:=f(u,v);
                END;
        END;

        FOR j:=1 TO k DO
        BEGIN
                FOR i:=1 TO k DO integral:=integral+y[j,i];
                integral:=integral-(y[1,j]+y[j,1]+y[k,j]+y[j,k])/2;
        END;
        integral:=integral+(y[1,1]+y[1,k]+y[k,1]+y[k,k])/4;

        integral:=integral*krok*krok;

        WRITELN('Vysledek integrace1 je ',integral:5:4);
        READLN;


        {Klasicky zpusob}

        integral:=0;

        FOR j:=1 TO k DO
        BEGIN
             v:=a+(j-1)*krok;
             z[j]:=0;
             FOR i:=1 TO k DO BEGIN
                 u:=a+(i-1)*krok;
                 z[j]:=z[j]+f(u,v);
             END;
             z[j]:=z[j]-(f(a,v)+f(b,v))/2;
        END;

        FOR j:=1 TO k DO integral:=integral+z[j];
        integral:=integral-(z[1]+z[k])/2;
        integral:=integral*krok*krok;

        WRITELN('Vysledek integrace2 je ',integral:5:4);
	READLN;
END.