/*
 * integrate.c - adaptive quadrature of the sqrt function.
 *
 * Alan J Rosenthal, December 2001.
 */

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


#define MAXINTERVALS 40

/* each interval is already really broken up into the two subtrapezoids */
struct {
    double left, mid, right;  /* values of x */
    double fleft, fmid, fright;  /* values of f(x) */
    double area;  /* area of the whole trapezoid from 'left' to 'right' */
    double leftarea, rightarea;  /* the two subtrapezoids (divided at mid) */
    double errest;  /* total error estimate for this interval */
} interval[MAXINTERVALS];
int nintervals;

extern void calc(int i);
extern int worst_error();


int main()
{
    int i;
    double h, result;

    /* start with ten intervals */
    /* we write into the eleventh struct in one place; that then gets ignored */
    nintervals = 20;
    h = 1.0 / nintervals;
    interval[0].left = 0;
    interval[0].fleft = sqrt(0);
    for (i = 0; i < nintervals; i++) {
        interval[i].mid = interval[i].left + h / 2;
        interval[i].fmid = sqrt(interval[i].mid);
        interval[i+1].left = interval[i].right = interval[i].left + h;
        interval[i+1].fleft = interval[i].fright = sqrt(interval[i].right);
        calc(i);
    }

    /* now keep subdividing the worst interval */
    for (; nintervals < MAXINTERVALS; nintervals++) {

	int worstinterval = worst_error();

        /*
         * for worstinterval's "left", "mid", "right":
         * make it into the interval [left,mid];
         * new interval is [mid,right] (in new slot in interval array).
         */
        interval[nintervals].left = interval[worstinterval].mid;
        interval[nintervals].right = interval[worstinterval].right;
        interval[nintervals].mid =
                (interval[nintervals].left + interval[nintervals].right) / 2;
        interval[nintervals].fleft = interval[worstinterval].fmid;
        interval[nintervals].fright = interval[worstinterval].fright;
        interval[nintervals].fmid = sqrt(interval[nintervals].mid);

        /* interval[worstinterval].left, fleft are ok */
        interval[worstinterval].right = interval[worstinterval].mid;
        interval[worstinterval].fright = interval[worstinterval].fmid;
        interval[worstinterval].mid =
                (interval[worstinterval].left + interval[worstinterval].right)
                / 2;
        interval[worstinterval].fmid = sqrt(interval[worstinterval].mid);

        /* we could avoid recalculating the overall areas of the new intervals.. */
        calc(worstinterval);
        calc(nintervals);
    }

/*
    for (i = 0; i < nintervals; i++)
        printf("interval %d goes from %lg to %lg, error estimate %lg\n", i,
                interval[i].left, interval[i].right, interval[i].errest);
 */

    /* add subtrapezoidal areas to get best value for overall integral */
    result = 0;
    for (i = 0; i < nintervals; i++)
        result += interval[i].leftarea + interval[i].rightarea;
    
    printf("result %g\n", result);
    return 0;
}


/* set area, leftarea, rightarea, errest for interval[i] */
void calc(int i)
{
    interval[i].area = (interval[i].right - interval[i].left)
            * ((interval[i].fleft + interval[i].fright) / 2);
    interval[i].leftarea = (interval[i].mid - interval[i].left)
            * ((interval[i].fleft + interval[i].fmid) / 2);
    interval[i].rightarea = (interval[i].right - interval[i].mid)
            * ((interval[i].fmid + interval[i].fright) / 2);
    interval[i].errest = (interval[i].area -
                        (interval[i].leftarea + interval[i].rightarea))
                      / 3;
}


/*
 * The following is a very inefficient algorithm, which the lab handout
 * said was ok.  We would normally keep this structure in order so that the
 * worst error was at the top of a data structure called a "heap".
 */
int worst_error()
{
    int i;
    double worst = 0;
    int worstinterval = 0;

    /* basic "max" kind of loop */
    for (i = 0; i < nintervals; i++) {
        if (fabs(interval[i].errest) > worst) {
            worst = fabs(interval[i].errest);
            worstinterval = i;
        }
    }

    return worstinterval;
}