void LU_decomposition (double *a, double *d, double *c, int n) {
    int i;

    c[0] = c[0]/a[0];

    for (i = 1; i < n; ++i) {
        a[i] -= d[i-1]*c[i-1];
        c[i] = c[i]/a[i];
    }
    a[n] = a[n] - d[n-1]*c[n-1];

    return;
}

void solve_line_system (double *a, double *d, double *c, double *b, double *x, int n) {
    int i;

    x[0] = b[0]/a[0];

    for (i = 1; i < n+1; ++i)
        x[3*i] = (b[i] - d[i-1]*x[3*(i-1)])/a[i];

    for (i = n-1; i > -1; --i)
        x[3*i] -= c[i]*x[3*(i+1)];

    return;
}