#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double f(double x, double a, double b)
{
    double drob_1 = (1.4*x)/(31.7*a+x*x);
    double drob_2 = (78*x*exp(-pow(x-1,2)/(2*a*a)))/(56.2+a);
    double drob_3 = (14*sin(pi*x/b))/(x*x+3.5*a);
    double drob_4 = (cos(pi*x/2*b))/(x*x*x+3);

    return (drob_1 + drob_2 + drob_3 + drob_4);
}

double rectangle_integral(pointFunc f, double a, double b, int n) {
      const double h = (b - a) / n;
      double sum = 0.0;
      #pragma omp parallel for reduction(+:sum)
      for (int i = 0; i < n; i++) {
        sum += f(a + i*h);
      }
      return (sum * h);
    }

int main(){
  double a, b, eps;
  double s1, s;
  int n = 1; //начальное число шагов

  a =0;
  b= 1000;
  eps = 0.001;

  s1 = rectangle_integral(f, a, b, n); //первое приближение для интеграла
  do {
    s = s1;     //второе приближение
    n = 2 * n;  //увеличение числа шагов в два раза,
                //т.е. уменьшение значения шага в два раза
    s1 = rectangle_integral(f, a, b, n);
  }
  while (fabs(s1 - s) > eps);  //сравнение приближений с заданной точностью
  cout << "\nИнтеграл = " << s1 << endl;
}