Лабы по МПО / lore / PS_TO_C2

#include<math.h>
#include<conio.h>
#include<stdlib.h>
#include<stdio.h>
//#include"Sampler.h"

/* integration by Simpson's method */

const double tol = 1.0E-6;

      char done;     //	: boolean;
      double sum, upper, lower, erf, twopi;

// SAMPLE;
double fx(double x)
{
// SAMPLE;
  return exp(-x*x);
}		/* function fx */

double dfx(double x)
{
 return (-2*x*exp(-x*x));
}


void simps(double lower,double upper,double tol)
/* numerical integration by Simpson's rule */
/* function is fx, limits are lower and upper */
/* with number of regions equal to pieces */
/* partition is delta_x, answer is sum */
{
      int i;
      double x, delta_x, even_sum, odd_sum, end_sum, end_cor, sum1;
      int pieces;

  pieces = 2;
  delta_x = (upper-lower)/pieces;
  odd_sum = fx(lower + delta_x);
  even_sum = 0.0;
  end_sum = fx(lower) + fx(upper);
  end_cor = dfx(lower) - dfx(upper);
  sum = (end_sum + 4.0 * odd_sum) * delta_x / 3.0;
 // printf("pieces: %-5d sum: %g\n", pieces, sum);
										  //writeln(pieces:5,sum);
										  ////?????????????????
  div_t tmp;
  int end_cicle;
// SAMPLE;
  do
  {
    pieces = pieces*2;
    sum1 = sum;
    delta_x = (upper-lower)/pieces;
    even_sum = even_sum + odd_sum;
    odd_sum = 0.0;
    tmp = div(pieces,2);
    end_cicle = tmp.quot;
//  SAMPLE;
    for (i = 1;i <= end_cicle;i++)
      {
	x = lower + delta_x * (2.0 * i - 1.0);
	odd_sum = odd_sum + fx(x);
      }
    sum = (7.0 * end_sum + 16.0 * odd_sum + 14.0 * even_sum + end_cor*delta_x) * delta_x/15.0;
   }
  while(!(sum!=sum1) && (abs(tol*sum1) <= abs(tol*sum) ) );
}	/* simps */


void main()  /* main program */
{
  clrscr();
  done = 0;	//false;
  twopi = 2.0 / sqrt(M_PI);
  lower = 0.0;

  do
  {
    printf("\nErf? \n");
    scanf("%lf",&upper);
//  SAMPLE;
    if (upper < 0.0)
      done = 1;//true
    else
//  SAMPLE;
      if (upper == 0.0)
       printf("Erf of 0.0 is 0.0\n");
      else	/* upper>0 */
	{
	    simps(lower, upper, tol);
	    erf = twopi * sum;
	    printf("Erf of %-7.2f, is %12.8f\n",upper,erf);
	}
  }while(done != 1);
}