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

#define NMAX 200
double p[NMAX+1];
double xx2[NMAX+1];
double pi;

double func(double x)
   {
    double y;
    y=sqrt(1.0-x*x);
    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)
   {
    double xx[NMAX+1],yy[NMAX+1];
    double ww2[NMAX+1];
    double i1,i2,i3,a,e1,e2,e3;
    int i,k,m,n;
    double h;
    pi=4.0*atan(1.0);
    printf("N     Trapezoidal rule     Simpson's rule     Gauss-Legendre rule \n");
    for(n=8;n<=20;n=n+4)
       {
        h=2.0/(double)n;
        xx[0]=-1.0;
        xx[n]=+1.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]);}
        e1=fabs(i1-pi/2.0);
        e2=fabs(i2-pi/2.0);
        e3=fabs(i3-pi/2.0);
        printf("N=%2d %15.10lf %15.10lf %15.10lf \n",n,e1,e2,e3);
       }
    return 0;
   }