extern float f(float);
void romberg(float r[10][10], float, float, int);
static int int_pow(int, int);

void romberg(float r[10][10], float a, float b, int M)
{
  int n, m, i;
  float h, s;
  
  h = b - a;
  r[0][0] = (f(a) + f(b)) * h / 2.0;

  for (n = 1; n <= M; n++) {
    h = h / 2.0;
    s = 0.0;

     for (i = 1; i <= int_pow(2, n - 1); i++) { 
   /* for (i = 1; i <= (1 << (n - 1)); i++) {*/
      s = s + f(a + (float)(2.0 * i - 1) * h);
    }

    r[n][0] = r[n-1][0]/2.0 + h * s;

    for(m = 1; m <= n; m++) {
      r[n][m] = r[n][m-1] + 
	 (float)(1.0/(int_pow(4,m)-1)) * (r[n][m-1] - r[n-1][m-1]); 
/* 	(float)(1.0/((1 << (2*m))-1)) * (r[n][m-1] - r[n-1][m-1]); */

      /* printf("<r[%d][%d] = %4.4f>\n", n, m, r[n][m]); */

    }
  }
}

static int int_pow(int base, int exp)
{
  int acc = 1;

  while (exp > 0) {
    acc = acc * base;
    exp = exp -1;
  }

  return acc;
}