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.
 * 1) include 
 * 2) include 
 * 3) include
 * 4) define DO(x,y) for(x=0;x 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");

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");
 * 1) define M(x,y) fprintf(o,"%6.2f %6.2f M\n",0.+x,0.+y);
 * 2) define L(x,y) fprintf(o,"%6.2f %6.2f L\n",0.+x,0.+y);
 * 3) define o(x,y) fprintf(o,"%6.2f %6.2f o\n",0.+x,0.+y);
 * 4) define O(x,y) fprintf(o,"%6.2f %6.2f O\n",0.+x,0.+y);
 * 5) define Q(x,y) fprintf(o,"%6.2f %6.2f O\n",0.+x,0.+y);

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