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

/*
   This is a test case for the adaptive, recursive numerical integration
   routine integral.  Note that the integrand is singular at x = 0, but
   that the integrator indeed finds the Riemann integral, although the
   integrator itself is not told of the location of the singularity.  To
   prevent numerical problems the exact singularity is excluded from the
   integration, but the absolute error introduced is less than 2.0E-19.
   This test case must be linked with s=10000 to prevent stack overflow.

     cc tstintegral s=10000

   The stack size depends upon the depth of recursion demanded.
   Compared to most practical problems, this test case is very severe.
*/

static int calls;

main()
{  double
      a,	    /* lower limit of integral */
      b,	    /* upper limit of integral */
      eps,	    /* tolerance on integration */
      ans,	    /* result of numerical integration */
      integral(),   /* integration routine */
      f1();	    /* function to be integrated */
   static int
      int_fail;     /* number of times integration failed to converge */

   a = -9.0;
   b = 1.0e4;
   /* Exact answer is 2 * (sqrt(b) + sqrt(|a|)) = 206. */

   for (eps = 0.1; eps >= 1.0e-10; eps = eps * 0.10)
   {  calls = 0;  /* For fun count the number of function calls. */
      ans = integral(f1, a, b, eps, &int_fail);
      printf("a= %8.3e b=%8.3e eps= %8.2e ans= %e int_fail=%3d calls=%5d\n",
	     a, b, eps, ans, int_fail, calls);
   }
} /* end main */

double f1(x)  /* function to be integrated */
   double x;  /* argement of function */
{  double
      f1,     /* function value to be returned */
      ax;     /* temporary variable = abs(x) */

   ++calls;	   /* This is only for display purposes; not needed. */
   ax = fabs(x);   /* Avoid calling fabs twice. */

   /* The following test avoids a possible zero divide; with the choice of
      limits of integration no truncation ever occurs.	The truncation
      error introduced in principle is 2.0E-19 in the integral.
   */

   if (ax > 1.0e-38) f1 = 1.0 / sqrt(ax);
   else 	     f1 = 1.0e19;
   return (f1);
} /* end f1 */

/* The output of the above test case follows:

a= -9.000E+0 b=1.000E+4 eps=  1.00E-1 ans= 2.002225E+2 int_fail=  0 calls=   55
a= -9.000E+0 b=1.000E+4 eps=  1.00E-2 ans= 2.060024E+2 int_fail=  3 calls=  427
a= -9.000E+0 b=1.000E+4 eps=  1.00E-3 ans= 2.060011E+2 int_fail=  6 calls=  607
a= -9.000E+0 b=1.000E+4 eps=  1.00E-4 ans= 2.060000E+2 int_fail=  9 calls=  979
a= -9.000E+0 b=1.000E+4 eps=  1.00E-5 ans= 2.060000E+2 int_fail= 15 calls= 1555
a= -9.000E+0 b=1.000E+4 eps=  1.00E-6 ans= 2.059999E+2 int_fail= 27 calls= 2563
a= -9.000E+0 b=1.000E+4 eps=  1.00E-7 ans= 2.059999E+2 int_fail= 42 calls= 4171
a= -9.000E+0 b=1.000E+4 eps=  1.00E-8 ans= 2.059999E+2 int_fail= 69 calls= 6763
a= -9.000E+0 b=1.000E+4 eps=  1.00E-9 ans= 2.059999E+2 int_fail=117 calls=11011
a= -9.000E+0 b=1.000E+4 eps= 1.00E-10 ans= 2.059999E+2 int_fail=195 calls=17971

   The correct answer is 206.
*/