Math.lab6

#include <stdio.h>
#include <math.h>
#define alpha1 2
#define beta1 0.5
#define gamma1 1.8765
#define alpha2 0.3
#define beta2 2
#define gamma2 3.4567

float S(const float *X, const float *F, float x, int N);
int find_i(const float *X, float x, int N);
void down_up(float a, float b, float c, float f, int N, float C[]);

int main(int argc, char **argv)
{
    const int N = 5;
    const float X[] = {0.2, 0.26, 0.28, 0.31, 0.32, 0.38};
    const float F[] = {1.2214, 1.2765, 1.3071, 1.3456, 1.3775, 1.4568};
    const float x = 0.32;
    int i = find_i(X,x,N);
    printf("\nS%d = %f\n\n", i, S(X,F,x,N));

    return 0;
}

float S(const float *X, const float *F, float x, int N)
{
    int i = find_i(X,x,N);
    float h0 = X[i] - X[i-1];
    float h1 = X[i+1] - X[i];
    float f = 6*( (F[i+1] - F[i])/h1 - (F[i] - F[i-1])/h0 );
    float C[N+1];
    down_up(h1, 2*(h0+h1), h0, f, N, C);

    float a = F[i];
    float d = (C[i] - C[i-1])/h0;
    float b = (h0/2)*C[i] - (pow(h0,2)/6)*d + (F[i]-F[i-1])/h0;
    float Si = a + b*(x-X[i]) + (C[i]/2)*pow((x-X[i]),2) + (d/6)*pow((x-X[i]),3);


    return Si;
}

void down_up(float a, float b, float c, float f, int N, float Y[])
{
    float X[N+1], Z[N+1];
    int i;

    X[0] = -1*beta1/alpha1;
    Z[0] = gamma1/alpha1;

    for (i = 0; i < N; i++) {
        X[i+1] = -1*a/(b+c*X[i]);
        Z[i+1] = (f-c*Z[i])/(b+c*X[i]);
    }

    Y[N] = (gamma2 - alpha2*Z[N]) / (beta2 + alpha2*X[N]);

    for (i = N; i > 0; i--) {
        Y[i-1] = X[i] * Y[i] + Z[i];
    }
}

int find_i(const float *X, float x, int N)
{
    int i;

    for (i = 0; i < N+1; i++)
        if (x <= X[i])
            break;
    return i;
}