辛普生公式和龙贝格递推的实现

如题，数值分析编程作业

 

代码：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<stack>
#include<queue>
#include<map>
#include<vector>
using namespace std;
typedef long long LL;
#define INF 1e18
#define endl '\n'
const int N = 5e2 + 10;
#define PI acos(-1) 
double eps;
double f(double x)
{
    if (x == 0)return 1;
    return sin(x) / x;
}

double Xinpson(double a, double b, int n)
{
    double h = (b - a) / n;
    double S = (f(b) - f(a));
    double x = a;
    for (int k = 0; k < n; k++)
    {
        S += 4 * f(x + h / 2) + 2 * f(x);
        x += h;
    }
    return S*(h/6);
}
double longbeige(double a,double b)
{
    double T1, T2, S1, S2, C1, C2, R1, R2;
    double h = (b - a);
    T1 = h/2 * (f(a) + f(b));
    int k = 1;
    while (1)
    {
        double x = a + h / 2;
        double S = 0;
        while (x < b)
        {
            S += f(x);
            x += h;
        }
        T2 = T1 / 2 + S * h / 2;
        S2 = T2 + (T2 - T1) / 3;
        if (k == 1)
        {
            k++;
            h /= 2;
            T1 = T2;
            S1 = S2;
            continue;
        }
        C2 = S2 + (S2 - S1) / 15;
        if (k == 2)
        {
            k++;
            h /= 2;
            C1 = C2;
            S1 = S2;
            T1 = T2;
            continue;
        }
        R2 = C2 + (C2 - C1) / 63;

        if (k == 3)
        {
            k++;
            h /= 2;
            C1 = C2;
            S1 = S2;
            T1 = T2;
            R1 = R2;
            continue;

        }
        if (fabs(R2 - R1) < eps)
            break;
        k++;
        h /= 2;
        C1 = C2;
        S1 = S2;
        T1 = T2;
        R1 = R2;

    }
    return R2;
}
int main()
{

    
    printf("%.8lf\n",Xinpson(0, 1, 4));
    eps = 1e-9;
    printf("%.10lf", longbeige(0, 1));
    

    return 0;
}