# Legendre-Gauss Quadrature formula/code

```// Code of the example of calculation of nodes and weights of the
// Legendre-Gauss Quadrature formula.
//  // Click at the image to see the description.
// No any external functions are required, but
//  your C++ compiler (or computer) should support the long double floating point arithmetics.
// It it does not, you still can reproduce the upper part of the figure,
//replacing long double to double.
// Copyleft 2008 by Dmitrii Kouznetsov
// You may extract the function gaule that evaluates nodes x and weights w and use it as you need.

#include <stdio.h>
#include <math.h>
#include<stdlib.h>
#define DO(x,y) for(x=0;x<y;x++)
```
```void gaule(long double x[], long double w[], int n)
{int i,j,m; long double z,t, a,b,c,q;
m=(n+1)/2; // add unity; then the same formula also for odd n
DO(i,m){z=cos(M_PI*(i+.5)/(n));
do{ a=1.; b=0.;  DO(j,n){ c=b; b=a; a=((2*j+(long double)1)*z*b-j*c)/(j+1); }
q=n*(z*a-b)/(z*z-(long double)1.); t=z; z=t-a/q;
} while (fabs(z-t) > 1.e-32);
x[i]=-z;     x[n-1-i]=z;
w[i]=2./((1-z*z)*q*q); w[n-1-i]=w[i];
}
}

void ado(FILE *O, int X, int Y)
fprintf(O,"%c%cBoundingBox: 0 0 %d %d\n",'%','%',X,Y);
fprintf(O,"/M {moveto} bind def\n");
fprintf(O,"/L {lineto} bind def\n");
fprintf(O,"/S {stroke} bind def\n");
fprintf(O,"/s {show newpath} bind def\n");
fprintf(O,"/C {closepath} bind def\n");
fprintf(O,"/F {fill} bind def\n");
fprintf(O,"/W {setlinewidth} bind def\n");
fprintf(O,"/RGB {setrgbcolor} bind def\n");}

#define N 1024
int main(void){ int i,n; long double t,s; long double *x,*w;
x=(long double *)malloc((size_t)(N)*sizeof(long double));
w=(long double *)malloc((size_t)(N)*sizeof(long double));
long double s0,p2;
s0=3.14159265358979323846264338327950288419716939937510;
printf("%40.30lf\n", s0);
s0/=2;
printf("%40.30lf\n", s0);
//s0=1.570796326794896619231321691639751442098; // fails at my computer
s0=1.57079632679489661923;			// so I do this
printf("%40.30lf\n", s0); 			// and such a way
p2=s0+6.123233995736767e-17;	// requires this stupid correction.
printf("%40.30lf\n", p2); // Direct assignment to long double fails.
s0=p2;
FILE *o; o=fopen("gaulegExample.eps","w"); ado(o,522,390);
```
```fprintf(o,"/o {.25 0 360 arc C F} bind def\n");
fprintf(o,"/O {.36 0 360 arc C S} bind def\n");
fprintf(o,"50 350 translate\n");
fprintf(o,"10 10 scale\n");
```
```#define M(x,y) fprintf(o,"%6.2f %6.2f M\n",0.+x,0.+y);
#define L(x,y) fprintf(o,"%6.2f %6.2f L\n",0.+x,0.+y);
#define o(x,y) fprintf(o,"%6.2f %6.2f o\n",0.+x,0.+y);
#define O(x,y) fprintf(o,"%6.2f %6.2f O\n",0.+x,0.+y);
#define Q(x,y) fprintf(o,"%6.2f %6.2f O\n",0.+x,0.+y);
M(0,1)L(0,-30) M(0,0)L(46,0)  fprintf(o,".1 W S\n");
fprintf(o,"/adobe-Roman findfont 2.1 scalefont setfont\n");
```
```M(-4.8,1.5) fprintf(o,"(lg(|error|))s\n");
for(n=5;n<50;n+=5){M(n,0)L(n,-30)}
for(n=10;n<40;n+=10){M(0,-n)L(45,-n)}
fprintf(o,".03 W S\n");
```
```for(n=10;n<50;n+=10){M(n-1.3,.3) fprintf(o,"(%2d)s\n",n);}
for(n=10;n<31;n+=10){M(-4,-n-.5) fprintf(o,"(%2d)s\n",-n);}
```
```fprintf(o,"/Times-italic findfont 2.1 scalefont setfont\n");
M(45,.5) fprintf(o,"(N)s\n");
```
```fprintf(o,".05 W\n");
double u,v ;
```
```fprintf(o,"0 0 1 RGB .12 W\n");
fprintf(o,"/adobe-Roman findfont 1.6 scalefont setfont\n");
//M(-4.8,15) fprintf(o,"(lg(|Resi|))s\n");
M(19,-27) fprintf(o,"(f[x_]=1/(3+x))s\n");
for(n=1;n<46;n++)
{gaule(x,w,n);
s=0; DO(i,n){t=x[i];t+=3;t*=t; s += w[i]/t;}
s-=(long double)1./4.; //  printf("%4d %14.3le ", n,s);
u=(double)s;  // printf("u= %14.3e ", u);
if( u>0) printf("\n u=%6.2e ;log(u)=%6.2f  ???\n",u,log(u));
if(u>0) {v=log(u)/log(10.); o(n,v); printf("n=%3d %8.3le  v=%9.3e pos\n",n,s,v);}
if(u<0) {v=log(-u)/log(10.);Q(n,v); printf("n=%3d %8.3le  v=%9.3e neg\n",n,s,v);}
}
```
```fprintf(o,"0 .8 0 RGB\n");  M(30.8,-23) fprintf(o,"(f[x_]=1/(1+x^2))s\n");
for(n=1;n<46;n++)
{gaule(x,w,n);
s=0; DO(i,n){t=x[i]; s += w[i]/(1+t*t);}
s-=s0;      // printf("%4d %16.8le", n,s);
u=(double)s;  // printf("u= %14.3e ", u);
if( u>0) printf("\n u=%6.2e ;log(u)=%6.2f  ???\n",u,log(u));
if(u>0) {v=log(u)/log(10.); o(n,v); printf("n=%3d %8.3le  v=%9.3e pos\n",n,s,v);}
if(u<0) {v=log(-u)/log(10.);Q(n,v); printf("n=%3d %8.3le  v=%9.3e neg\n",n,s,v);}
}
```
```fprintf(o,"0 0 0 RGB\n");   M(3,-29.5) fprintf(o,"(f[x_]=x^32)s\n");
for(n=1;n<46;n++)
{gaule(x,w,n);
s=0; DO(i,n){t=x[i]; t*=t; t*=t; t*=t; t*=t; s += w[i]*t;}
printf(" %14.3le ",s-(long double)2./17.);
s-=(long double)2./17.;
u=(double)s;  printf("u= %14.3e ", u);
if( u>0) printf("\n u=%6.2e ;log(u)=%6.2f  ???",u,log(u));
if(u>0) {v=log(u)/log(10.); o(n,v); printf("n=%3d %8.3le  v=%9.3e pos\n",n,s,v);}
if(u<0) {v=log(-u)/log(10.);Q(n,v); printf("n=%3d %8.3le  v=%9.3e neg\n",n,s,v);}
}
```
```fprintf(o,"1 0 0 RGB\n");  M(25.4,-6.5) fprintf(o,"(f[x_]=Sqrt[1+x^2])s\n");
for(n=1;n<48;n++)
{gaule(x,w,n);
s=0; DO(i,n){t=x[i]; s += w[i]*sqrt((long double)1-t*t);}
printf(" %14.3le\n",s-M_PI_2);
s-=M_PI_2;
u=(double)s;  printf("u= %14.3e ", u);
if( u>0) printf("\n u=%6.2e ;log(u)=%6.2f  ???",u,log(u));
if(u>0) {v=log(u)/log(10.); o(n,v); printf("n=%3d %8.3le  v=%9.3e pos\n",n,s,v);}
if(u<0) {v=log(-u)/log(10.);Q(n,v); printf("n=%3d %8.3le  v=%9.3e neg\n",n,s,v);}
}
```
```fprintf(o,"showpage\n%cTrailer\n",'%');
fclose(o);
```
```       free((char*)(w));
free((char*)(x));
//system("open gaulegExample.eps");
//system("ps2pdf gaulegExample.eps");
//getchar(); system("killall Preview");
}
// End of copylefted source
```
 Some content on this page may previously have appeared on Citizendium.