/*-------------------------------------------------------
  8次のニュートン・コーツによる数値計算プログラム
  -------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

#define Tmin   0
#define Tmax   1000
#define CALC_N 10000
#define F(x)   (1.84*(1+0.39*pow(x/2200.,0.71)))   /* 非積分関数 */

///////////////////////////////////////////////////////
// 関数宣言
///////////////////////////////////////////////////////
void suuhyo(double y[],double a,double b,int n);
double nc8(double y[],double a,double b,int n);
///////////////////////////////////////////////////////

int main()
{
    double start_x=Tmin;
    double end_x  =Tmax;
    /* 非積分関数用の領域確保 */
    double *y=(double*)calloc(CALC_N,sizeof(double));
    /* 非積分用の関数を計算   */
    suuhyo(y,start_x,end_x,CALC_N);
    /* 積分値の計算・表示     */
    printf("square=%lf\n",nc8(y,start_x,end_x,CALC_N));

    /* メモリの解放 */
    free(y);
    return 0;
}

/* テストデータ作成 */
void suuhyo(double y[],double a,double b,int n)
{
    int i;
    double x,dx;
    dx=(b-a)/(double)n;
    for(i=0;i<=n;++i){
        x=a+(double)i*dx;
        y[i]=F(x);
    }
}

/* ８次のNC */
double nc8(double y[],double a,double b,int n)
{
    int        i;
    double  dx,s,s0,s1,s2,s3,s4;
    dx=(b-a)/(double)n;
    s0=s1=s2=s3=0.0;
    for(i=0;i<=n-8;i+=8){
        s0+=y[i  ]+y[i+8];        
        s1+=y[i+1]+y[i+7];        
        s2+=y[i+2]+y[i+6];
        s3+=y[i+3]+y[i+5];
        s4+=y[i+4];
    }
    s=989.0*s0+5888.0*s1-928.0*s2
        +10496.0*s3-4540.0*s4;
    return(s*dx*4.0/14175.0);
}