#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define NMAX 10000

double xx[NMAX+1],yy[NMAX+1];
double xx2[NMAX+1],ww2[NMAX+1];
double p[NMAX+1];
double i1,i2,i3,pi;

double func(double x)
   {
   	double y;
   	y=8.0/3.0*sqrt(pow((1.0-x*x),3));
   	return y;
   }

void legendre(int n,double x)
   {
    int j;
    p[0]=1.0;
    p[1]=x;
    for(j=1;j<=n-1;j++){p[j+1]=((double)(2*j+1)*x*p[j]-(double)j*p[j-1])/(double)(j+1);}
    return;
   }

void zeroten(int n)
   {
    double delta,eps,a,b,c,fa,fb,fc,diff;
    int i,m;
    delta=1.0e-10;
    eps=1.0e-10;
    m=n/2;
    for(i=1;i<=m;i++)
       {
        a=sin(pi*(double)(n-2*i)/(double)(2*n));
        b=sin(pi*(double)(n-2*i+2)/(double)(2*n));
        legendre(n,a);
        fa=p[n];
        legendre(n,b);
        fb=p[n];
        if(fa*fb>0.0){printf("error a,b,i==%d \n");}
        fc=1.0;
        while(fabs(fc)>delta&&fabs(a-b)>eps)
           {
            diff=fabs(a-b);
            c=(a+b)/2.0;
            legendre(n,c);
            fc=p[n];
            if(fc*fa>0.0){a=c;fa=fc;}
            else{b=c;fb=fc;}
           }
        xx2[n-i+1]=c;
        xx2[i]=-c;
       }
    return;
   }

int main(void)
   {
   	int ns,n,m,i,k;
   	double h,a;
   	pi=4.0*atan(1.0);
   	printf("N     Trapezoidal rule     Simpson's rule     Gauss-Legendre rule \n");
   	for(ns=1;ns<=10;ns++)
   	{
     if(ns==1){n=2;}
     if(ns==2){n=3;}
     if(ns==3){n=4;}
     if(ns==4){n=5;}
     if(ns==5){n=6;}
     if(ns==6){n=10;}
     if(ns==7){n=20;}
     if(ns==8){n=100;}
     if(ns==9){n=1000;}
     if(ns==10){n=10000;}
     h=2.0/(double)n;
     xx[0]=-1.0;
     xx[n]=1.0;
     yy[0]=0.0;
     yy[n]=0.0;
     for(i=1;i<=n-1;i++)
     {
      xx[i]=xx[0]+(double)i*h;
      yy[i]=func(xx[i]);
     }
     i1=yy[0]+yy[n];
     for(i=1;i<=n-1;i++){i1=i1+2.0*yy[i];}
     i1=i1*h/2.0;
     m=n/2;
     i2=0.0;
     for(i=0;i<=m-1;i++){i2=i2+yy[2*i]+4.0*yy[2*i+1]+yy[2*i+2];}
     i2=i2*h/3.0;
     zeroten(n);
     for(i=1;i<=n;i++)
     {
      a=0.0;
      legendre(n,xx2[i]);
      for(k=0;k<=n-1;k++){a=a+(double)(2*k+1)*p[k]*p[k];}
      ww2[i]=2.0/a;
     }
     i3=0.0;
     for(i=1;i<=n;i++){i3=i3+ww2[i]*func(xx2[i]);}
     printf("N= %d %18.14lf %18.14lf %18.14lf \n",n,i1,i2,i3);
   	}
    return 0;
   }