Как точно рассчитать данные двойного типа для решения интерполяции кубического сплайна в C?

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SIZE 11

double cubicspline(double val);
void TDMA(void);

double x[SIZE] = { -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };
double y[SIZE] = { 0.038, 0.058, 0.1, 0.2, 0.5, 1.0, 0.5, 0.2, 0.1, 0.058, 0.038 };
double M[SIZE];

void main(void) {
    FILE *fp1, *fp2;
    double val;
    int i;
    fp2 = fopen("cubicSpline_output.txt", "w");
    val = -1;
    while(val<=1.0) 
    {
        printf("------------------------------\n");
        printf("val : %lf\n",val);
        printf("spline result : %.12lf\n", cubicspline(val));
        fprintf(fp2, "%lf\t%lf\n", val, cubicspline(val));
        printf("------------------------------\n");
        val += 0.1;
    }
    fclose(fp2);
    system("pause");
}

void TDMA() {
    double A[SIZE-1], B[SIZE-1], C[SIZE-1], R[SIZE-1];
    double f[SIZE], g[SIZE];

    int i;

    for (i = 1; i < SIZE - 1; i++) 
    {
        A[i] = (x[i] - x[i - 1])/6;
        B[i] = (x[i + 1] - x[i - 1])/3;
        C[i] = (x[i + 1] - x[i])/6;
        R[i] = (y[i + 1] - y[i]) / 6 * C[i] - (y[i] - y[i - 1]) / 6 * A[i];

    }

    for (i = 1; i < SIZE - 2; i++) 
    {
        f[i + 1] = B[i + 1] - (C[i] * A[i + 1]) / B[i];
        g[i + 1] = R[i + 1] - (R[i] * A[i + 1]) / B[i];
    }

    M[0] = 0;
    M[SIZE - 1] = 0;
    M[SIZE - 2] = g[SIZE - 2] / f[SIZE - 2];

    for (i = SIZE - 3; i > 1; i--) 
    {
        M[i] = (g[i + 1] / f[i + 1]) - (C[i + 1]*M[i + 1]);
    }  

    M[1] = (R[1] - C[1]*M[2])/B[1];
}

double cubicspline(double val)
{
    int i, j;
    double result=0;

    TDMA();

    for (i = 1; i < SIZE; i++)
    {
        if (x[i] >= val) {
            j = i;
            break;
        }
    }
    result = ((pow((x[j] - val), 3.0) * M[j - 1]) + (pow((val - x[j - 1]), 3.0)) * M[j])/ 6 * (x[j] - x[j - 1]) + ((x[j] - val)*y[j - 1] + (val - x[j - 1])*y[j]) / (x[j] - x[j - 1]) - ((x[j] - x[j - 1])*((x[j] - val)*M[j - 1] + (val - x[j - 1])*M[j])) / 6;
    return result;
}

Вот моя часть исходного кода, написанная в среде Visual Studio 2017. Глобально инициализированный массив x,y получает данные. Я должен сделать интерполяцию кубическим сплайном с этими данными. Но после выполнения результаты сплайна в заданных точках немного отличаются. Я думаю, что нет никакой логической проблемы в сплайне и функции TDMA. Возможно, эти проблемы связаны с хранением иррациональных чисел. Можете ли вы, ребята, объяснить, почему возникает эта ошибка и как исправить код?