Współczynniki wielomianu interpolacji Lagrange'a - część wyników to -nan

Dobry wieczór. Jestem w trakcie pisania programu obliczającego wartości współczynników wielomianu interpolacji Lagrange'a. Napisałem kod, który dla dwóch (dość małych) zestawów danych działał prawidłowo. Niestety, w przypadku trzeciego, trochę większego zestawu, ostatnie 10 wartości wynosi -nan, a pozostałe to jedynki. Tablice z danymi wejściowymi są prawidłowe, nigdzie też nie wyłapałem dzielenia przez zero. Zmiana z double na long double skutkuje tym, że wszystkie wyniki wynoszą -nan. Skończyły mi się już pomysły, stąd też moja prośba o pomoc w znalezieniu błędu.

Oto kod programu:

#include <iostream>
#include <iomanip>

double lj(const double xn[], int n, int idx, int x)   //Oblicza wartość fj
{   //(x-x0)(x-x1)...(x-x(i-1))(x-x(i+1))...(x-xn)/(xj-x0)(xj-x1)...(xj-x(i-1))(xj-x(i+1))...(xj-xn)
    double wynik=1; //wartość funkcji

    //złożoność O(n)
    for(int t=0; t<n; ++t)
    {
        if(t!=idx)
        {
            wynik*=(x-xn[t])/(xn[idx]-xn[t]);
        }
    }

    return wynik; 
}

template<int N>
struct x_points
{
    constexpr x_points() : array()  //konstruktor constexpr
    {
        for(int i=0; i<N; ++i)  //węzły obliczane są w czasie kompilacji
        {
            array[i]=-1+i*(1/32.);
        }
    }
    double array[N];
};

double* valuesOfFunction(const double x[], int size)    //wartości w węzłach
{
    double *tab=new double[size];

    for(int i=0; i<size; ++i)
    {
        tab[i]=1/(1+5*x[i]*x[i]);
    }

    return tab;
}

int main()
{
    const int precision=4;
    const int zero_val=0;
    const int x_max=1;
    const int x_min=-1;
    constexpr int size=((x_max-x_min)/(1/32.))+1;   //rozmiar tablicy, liczba węzłów, wartości oraz współczynników

    constexpr auto xv=x_points<size>();
    const double *xn=xv.array;  //węzły
    
    /*for(int i=0; i<size; ++i)
    {
        std::cout<<xn[i]<<std::endl;
    }*/

    const double *fx=valuesOfFunction(xn, size);   //wartości w węzłach

    /*for(int i=0; i<size; ++i)
    {
        std::cout<<fx[i]<<std::endl;
    }*/

    //constexpr int size=std::size(xn);   //rozmiar tablicy, liczba węzłów, wartości oraz współczynników
    
    double temp[size];  //tablica dla wartości w węzłach kolejnych wielomianów
    double a[size]; //współczynniki wielomianu w postaci naturalnej

    int xn_zero_value_index=-1;   //indeks miejsca, w którym tablica xn zawiera element zerowy

    for(int i=0; i<size; ++i)
    {
        if(xn[i]==0)
        {
            xn_zero_value_index=i;
        }
        a[0]+=fx[i]*lj(xn, size, i, zero_val); //wyraz wolny
    }
    
    if(xn_zero_value_index!=-1)
    {
        temp[xn_zero_value_index]=a[0];
    }

    for(int j=0; j<size; ++j)   //oddzielnie ze względu na brak konieczności interpolacji wartości
    {   // (fj-a0)/xj
        
        if(xn[j]!=0)
        {
            temp[j]=(fx[j]-a[0])/xn[j]; //wartości y1(x) w kolejnych węzłach
        }
        a[1]+=temp[j]*lj(xn, size, j, zero_val);   //suma kolejnych wyrazów wielomianu interpolacyjnego w zerze
    }

    //temp[xn_zero_value_index]=a[1];

    //pozostałe współczynniki wg wzoru:
    // yi(xj)=(y{i-1}(xj)-a{i-1})/xj
    for(int k=2; k<size; ++k)   //pętla po kolejnych współczynnikach
    {   
        if(xn_zero_value_index!=-1)
        {
            temp[xn_zero_value_index]=a[k-1];
        }

        for(int l=0; l<size; ++l)   //pętla po kolejnych wartościach w węzłach
        {
            if(xn[l]!=0)
            {
                temp[l]=(temp[l]-a[k-1])/xn[l];
            }
            a[k]+=temp[l]*lj(xn, size, l, zero_val);
        }
    }

    for(int m=0; m<size; ++m)   //wyniki
    {
        std::cout<<std::setprecision(precision)<<std::fixed<<"a"<<m<<" = "<<a[m]<<std::endl;
    }

    return 0;
}

Ej no, Pan raczy sobie żartować, zaczyna się ładnie, wszystko w małych funkcjach, a potem ten wielgaśny main. Zrobiłbyś chociaż jakieś funkcje mnożące wektory przez wektor i dodające wektory do wektorów, to by się zrobiło czytelniej (a i może błędy by zniknęły).

Znalezienie problemu może być trudne jeśli patrzy się tylko na kod ponieważ problemem mogą być błędy numeryczne. Proponowałbym wprowadzić niezmienniki do kodu (invariants) w oparciu o assert. W szczególności testy czy rezultat wyrażenia nie jest NaN. Do tego możesz wykorzystać funkcję isnan z math.h. Potem debugger i jesteś w domu.
Kilka uwag przy okazji:

• tworzysz różne tablice na stosie, co prawda nie sprawdziłem ich rozmiaru ale jest to ryzykowna praktyka (okryty złą sławą stack overflow). Jeśli ich rozmiary są znaczne to alokuj je na stercie
• nie widzę zwalniania const double *fx. W tym szczególnym przypadku nie jest to wyciek pamięci (alokacja tylko jeden raz a po zakończeniu programu pamięć i tak będzie zwolniona przez OS) ale nie jestem pewny czy zrobiłeś to intencjonalnie, a tak czy siak jest to ryzykowna technika