#include <iostream>
#include <cmath>
#include <vector>

using namespace std;


vector<vector<double>> romber(double (*function)(double), double a, double b, int M){
    double h = b - a;
    vector<vector<double>> R (M+1, vector<double> (M+1, 0));
    R[0][0] = h * (function(a) + function(b)) / 2.0;
    for (int n = 0; n < M; n++){
        h = h / 2.0;
        double inner_sum = 0;
        for (int i = 1; i <= pow(2, n); i++){
            inner_sum += function(a + (2*i - 1) * h);
        }
        R[n+1][0]= 1.0 / 2.0 * R[n][0] + h * inner_sum;
        for(int m = 0; m < n; m++){
            R[n+1][m+1] = R[n+1][m] + (R[n+1][m]-R[n][m]) / (pow(4, m+1) - 1);
        }
    }
    //print
    for (int i = 0; i<R.size(); i++){
        for(int j = 0; j<R[i].size(); j++){
            cout << fixed << R[i][j] << ' ';
        }
        cout << endl;
    }
    return R;
}


int main(){
    vector<vector<double>> result;
    romber(&sin, 0, M_PI, 4);
    return 0;
}