GauLegExample/code
Jump to navigation
Jump to search
//// can be generated with the code below. // No any external funcitons are required, // but your C++ compiler 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
//