HP 15C and INT(1/√(1-x),0,1)

[TheKaneB]

Back to the original topic:

I used a C implementation of the Romberg algorithm to find the integral of
f(x) = 1 / sqrt(1 - x) over the interval 0 - 1

It runs for a relatively long time (several seconds on my 3 GHz Intel i5) and it came up with this result:

result = 2.000052194143781214563660
function evaluations = 1073741825 (1 billion of function evaluations!)

I set an accuracy of 0.0001 with max steps = 100. I also used the interval from 0 to 0.99999999999 or else it would evaluate to NAN.

If i set an accuracy of 0.01 I get this:
result = 2.005444550190635499831160
function evaluations = 16777217

Code:


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

void
dump_row(size_t i, double *R){
   printf("R[%2zu] = ", i);
   for(size_t j = 0; j <= i; ++j){
      printf("%.12f ", R[j]);
   }
   printf("\n");
}

double
romberg(double (*f/* function to integrate */)(double), double /*lower limit*/ a, double /*upper limit*/ b, size_t max_steps, double /*desired accuracy*/ acc){
   double R1[max_steps], R2[max_steps]; //buffers
   double *Rp = &R1[0], *Rc = &R2[0]; //Rp is previous row, Rc is current row
   double h = (b-a); //step size
   Rp[0] = (f(a) + f(b))*h*.5; //first trapezoidal step

   dump_row(0, Rp);

   for(size_t i = 1; i < max_steps; ++i){
      h /= 2.;
      double c = 0;
      size_t ep = 1 << (i-1); //2^(n-1)
      for(size_t j = 1; j <= ep; ++j){
         c += f(a+(2*j-1)*h);
      }
      Rc[0] = h*c + .5*Rp[0]; //R(i,0)

      for(size_t j = 1; j <= i; ++j){
         double n_k = pow(4, j);
         Rc[j] = (n_k*Rc[j-1] - Rp[j-1])/(n_k-1); //compute R(i,j)
      }

      //Dump ith column of R, R[i,i] is the best estimate so far
      dump_row(i, Rc);

      if(i > 1 && fabs(Rp[i-1]-Rc[i]) < acc){
         return Rc[i-1];
      }

      //swap Rn and Rc as we only need the last row
      double *rt = Rp;
      Rp = Rc;
      Rc = rt;
   }
   return Rp[max_steps-1]; //return our best guess
}

int functionCounter = 0;

double myfunc(double x) {
    functionCounter++;
    double l = 1 / sqrt( 1 - x );
    return l;
}

int main() {
    printf("result = %.24f\nfunction evaluations = %d\n", romberg(&myfunc, 0, 0.99999999999, 100, 0.0001), functionCounter);
}

The code source is from wikipedia: https://en.wikipedia.org/wiki/Romberg%27s_method