🚀 go-pugleaf

// Since there have been a few requests for Gauss-Legendre
// nodes & weights, I thought I'd post this here. I am not
// including the routine "myclib_error"; all it does is
// print the error message to cerr and to cout, then stop.
// I have compared the output of this routine to the tables
// in Abramowitz & Stegun, and they seem to be ok for all
// the orders listed. This code is not copyright-protected,
// so you can hack it to your heart's content.
// (You will anyway, I know :-) ) Enjoy!

#include <math.h>

// This routine calculates the value of the Nth-order Legendre
// polynomial and its first and second derivatives at the point X.
// The Legendre polynomials satisfy the recursion relation
//	P_n+1   = A*X*P_n - B*P_n-1			(polynomial)
//	P_n+1'  = A*(X*P_n' + P_n) - B*P_n-1'		(first derivative)
//	P_n+1'' = A*(X*P_n'' + 2*P_n') - B*P_n-1''	(second derivative)
//	where A and B are
//	A = (2*n + 1)/(n+1)
//	B = n/(n+1)
// On output, PNP contains the value of the polynomial, PND contains
// the value of the first derivative, and PNA contains the value of
// the second derivative.

void plegpda(int n, double x, double& pnp, double& pnd, double& pna)
{
    int k;
    double a, b, pm, dm, am, pc, dc, ac, pp, dp, ap;

    if (n < 0) {
	myclib_error("bad input to plegpda");
    } else if (n == 0) {
	pnp = 1.0;
	pnd = 0.0;
	pna = 0.0;
    } else if (n == 1) {
	pnp = x;
	pnd = 1.0;
	pna = 0.0;
    } else {
	pm = 1.0;
	dm = 0.0;
	am = 0.0;
	pc = x;
	dc = 1.0;
	ac = 0.0;
	for (k = 1; k < n; ++k) {
	    a = ((double) (2*k + 1))/((double) (k + 1));
	    b = ((double) k)/((double) (k + 1));
	    pp = a*x*pc - b*pm;
	    dp = a*(x*dc + pc) - b*dm;
	    ap = a*(x*ac + 2*dc) - b*am;
	    pm = pc;
	    pc = pp;
	    dm = dc;
	    dc = dp;
	    am = ac;
	    ac = ap;
	}
	pnp = pc;
	pnd = dc;
	pna = ac;
    }
}

// This routine generates the Gauss-Legendre points, defined by
//
//	P_{N+1}(z_i) = 0		i = 0, ... , N
//
// and the weights for Gauss-Legendre integration, given by
//
//			    2
//	W_i = -----------------------------
//	      (1 - z_i^2)*(P'_{N+1}(z_i))^2

static const double PI = 3.14159265358979323846264338327950288;

void glzw(int n, double *z, double *w, double eps)
{
    int i, np;
    double pin, zz, pp, pd, pa;

    if (n >= 0) {
	np = n + 1;
	pin = PI/((double) np);
	for (i = 0; i <= np/2; ++i) {
	    zz = cos(pin*(0.5 + ((double) i)));
	    pa = 1.0;
	    while (abs(pa) > eps) {
		plegpda(np, zz, pp, pd, pa);
		pa = -pp/pd;
		zz += pa;
	    }
	    z[i] = -zz;
	    z[n-i] = zz;
	}
	for (i = 0; i <= np/2; ++i) {
	    plegpda(np, z[i], pp, pd, pa);
	    w[i] = 2.0/((1.0 - z[i]*z[i])*pd*pd);
	    w[n-i] = w[i];
	}
    } else {
	myclib_error("negative N in glzw");
    }
}
-- 
Uwe Hollerbach          uh@acm.princeton.edu  or  uh@alumni.caltech.edu
"In His infinite mercy, Allah does not subtract from one's alloted span those
hours which are spent in contemplation of net.news. All praise be to Allah."
-- an obscure commentator on the early work of Al-Khowariszmi