gaussian quadrature

here is an implantation of gaussian quadrature using lagrange polynomials. rather than generating the sample points and weights separately, this routine generates them as it goes along integrating. this is good if you dont want to store them and bad if you want to repeat the process on a different function, but with the same number of samples.

the coefficients depend on the number of samples so in practice you always want to vary the number of samples until you get enough accuracy.

i thought originally, this might make a good candidate for a calculator implementation, but unfortunately, its way too slow having to do the work of gauss-legendre weights and coefficients as well (even though it doesnt require the memory to store them).

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

typedef double (*Fn)(double);

double gauss(int n, double a, double b, Fn f)
{    
    int     m, i, j; 
    double  t, t1, pp, p1, p2, p3;
    double  pi = 3.1415926535897932385;
    double  eps = DBL_EPSILON;
    double ba1, ba2;
    double w, q;
    double ca, sa, cd, sd;

    ba1 = (b-a)/2.0;
    ba2 = (b+a)/2.0;
    q = 0;

    m = (n+1)/2;

    ca = cos(pi*0.75/(n + 0.5));
    sa = sin(pi*0.75/(n + 0.5));
    cd = cos(pi/(n + 0.5));
    sd = sin(pi/(n + 0.5));

    for(i=0; i<m; ++i) {  
        t  = ca;
        t1 = ca*cd - sa*sd;
        sa = sa*cd + ca*sd;
        ca = t1;
        t1 = 1;
        while((fabs(t-t1))>eps) { 
            p1 = 1.0;
            p2 = 0.0;
            for(j=1; j<=n; j++) {
                p3 = p2;
                p2 = p1;
                p1 = ((2*j-1)*t*p2-(j-1)*p3)/j;
            }
            pp = n*(t*p1-p2)/(t*t-1);
            t1 = t;
            t  = t1 - p1/pp;
        }

        w = 2.0*ba1/((1-t*t)*pp*pp);
        q += w*((*f)(-t*ba1+ba2)+ (*f)(t*ba1+ba2));
    }
    return q;
}