#include <iostream>

using namespace std;

void natural_cubic_spline(double x[], double a[], int n, double b[], double c[], double d[]) {
    double h[n], alfa[n], l[n+1], myu[n], z[n+1];

    // Step 1:
    for (int i = 0; i < n; i++) {
        h[i] = x[i+1] - x[i];
    }

    // Step 2
    for (int i = 1; i < n; i++) {
        alfa[i] = (3/h[i])*(a[i+1]-a[i]) - (3/h[i-1])*(a[i]-a[i-1]);
    }

    // Step 3
    l[0] = 1;
    myu[0] = 0;
    z[0] = 0;

    // Step 4
    for (int i = 1; i < n; i++) {
        l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*myu[i-1];
        myu[i] = h[i] / l[i];
        z[i] = (alfa[i] - h[i-1]*z[i-1]) / l[i];
    }

    // Step 5
    l[n] = 1;
    z[n] = 0;
    c[n] = 0;

    // Step 6
    for (int j = n-1; j >= 0; j--) {
        c[j] = z[j] - myu[j]*c[j+1];
        b[j] = (a[j+1]-a[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3;
        d[j] = (c[j+1]-c[j])/(3*h[j]);
    }

    // Print the results
    cout << "\nj" << "\tx[j]" << "\ta[j]" << "\t\tb[j]" << "\t\tc[j]" << "\t\td[j]";
    for (int j = 0; j < n; j++) {
        cout << "\n" << j << "\t" << x[j] << "\t" << a[j] << "\t\t" << b[j] << "\t\t" << c[j] << "\t\t" << d[j];
    }
}

int main() {
  double x[100] = {0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5.0, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3};
  double y[100] = {1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.2};

  int n = 21; // panjang array diatas.

  double b[n], c[n+1], d[n];

  natural_cubic_spline(x, y, n, b, c, d);

  return 0;
}