GauLegExample/code

// Code of the example of calculation of nodes and weights of the
// Gaussian-Legendre 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!PS-Adobe-2.0 EPSF-2.0\n",'%');
	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
//