#include<stdio.h>

#include<math.h>

#define f(x) 1.0/(1.0+(x)*(x))   /*  ˾·ÌÞݶݽ³  */

#define   a  0.0                 /*  ¶ÀÝ         */

#define   b  1.0                 /*  ¼Þ®³ÀÝ      */

void main()

{

    long   n;                     /*  ÌÞÝ¶Â É ½³  */

    double h;                     /*  ¸¶Ý É ÊÊÞ   */

    double P21(long n,double h);

    double P22(long n,double h);

    double P23(long n,double h);

    printf("\n---- Definite integral ----\n\n");

    printf("f(x)=1/(1+x^2),(0<x<1) \n\n");

    printf("   ÌÞݶ    ¸¶Ý É ÊÊÞ :     ¸ÌÞÝ·­³¾·        ̸ºÞ³ÀÞ²¹²       T.Simpson \n\n");

    for(n=2;n<100000;n=n*2)

    {

        h=(b-a)/n;

        printf(" %7d   %9.6f  :",n,h);

        printf("  %14.10f   %14.10f   %14.10f \n",P21(n,h),P22(n,h),P23(n,h));

    }

    printf("\n");

}

double P21(long n,double h)

{

    long   i;

    double x=a,y=0.0;

    for(i=0;i<n;i++)

    {

        x=x+h;

        y=y+f(x);

    }

    return h*y;

}

double P22(long n,double h)

{

    long   i;

    double x=a,y=0.0;

    for(i=0;i<n;i++)

    {

        y=y+(f(x)+f(x+h))/2.0;

        x=x+h;

    }

    return h*y;

}

double P23(long n,double h)

{

    long   i;

    double x=a,y=0.0;

    for(i=0;i<n/2;i++)

    {

        y=y+(f(x)+4.0*(f(x+h))+f(x+h+h))/3.0;

        x=x+h+h;

    }

    return h*y;

}