Całkowanie numeryczne

#include <iostream>
#include <ctime>
#include <cmath>
#include <cstdlib>

using namespace std;
double z;
double f1(double), f2(double), f3(double), f4(double);
typedef double(*funct)(double);
void prostokat(funct, double, double);
void trapez(funct, double, double);
void monte(funct, double, double);
int main(int argc, char** argv) {

    double pi=M_PI;

    cout<<"1) f(x)=sinx dla x=<0, pi>: \n\na) Metoda prostokatow (usredniona): ";
    prostokat(f1, 0, pi);
    cout<<"\nb) Metoda trapezow: ";
    trapez(f1, 0, pi);
    cout<<"\nc) Metoda Monte Carlo: ";
    monte(f1, 0, pi);

    cout<<"\n\n\n2) f(x)=1/(x^2+1)  dla x=<0, 1>: \n\na) Metoda prostokatow (usredniona): ";
    prostokat(f2, 0, 1);
    cout<<"\nb) Metoda trapezow: ";
    trapez(f2, 0, 1);
    cout<<"\nc) Metoda Monte Carlo: ";
    monte(f2, 0, 1);

    cout<<"\n\n\n3) f(x)=sqrt(2x)   dla x=<0, 8>: \n\na) Metoda prostokatow (usredniona): ";
    prostokat(f3, 0, 8);
    cout<<"\nb) Metoda trapezow: ";
    trapez(f3, 0, 8);
    cout<<"\nc) Metoda Monte Carlo: ";
    monte(f3, 0, 8);

    cout<<"\n\n\n4) f(x)=x*sin4x    dla x=<0, pi/4>: \n\na) Metoda prostokatow (usredniona): ";
    prostokat(f4, 0, pi/4);
    cout<<"\nb) Metoda trapezow: ";
    trapez(f4, 0, pi/4);
    cout<<"\nc) Metoda Monte Carlo: ";
    monte(f4, 0, pi/4);












    return 0;
}
double f1(double x) {
    return sin(x);
}
double f2(double x) {
    double a=pow(x,2);
    return 1.0/(a+1);
}
double f3(double x) {
    return sqrt(2*x);
}
double f4(double x) {
    return x*sin(4*x);
}
void prostokat(funct f, double prze, double dzial) {
    double x=(dzial-prze)*1.0/1000000;
    double a=prze;
    double b=a+x;
    double p1, p2, wynik;
    while(b<=dzial) {
        p1+=fabs(f(b))*x;
        p2+=fabs(f(a))*x;
        wynik+=((p1+p2)*1.0/2);
        a+=x;
        b+=x;
    }
    cout<<wynik;
}
void trapez(funct f, double prze, double dzial) {
    double x=(dzial-prze)/1000000;
    double a=prze, b=a+x, p, c;
    while(b<=dzial) {
        c=(fabs(f(a))+fabs(f(b)));
        p+=c*x*(1.0/2);
        a+=x;
        b+=x;
    }
    cout<<p;
}
void monte(funct f, double prze, double dzial) {
    void (srand(time(NULL)));
    double a, x, wart;
    for(int i=0; i<10000000; i++) {
        a=(rand()*1.0)/RAND_MAX;
        x=a*(dzial-prze)+prze;
        wart+=fabs(f(x));
    }
    cout<<((wart*1.0)/10000000)*(dzial-prze);


}