Интерполяционный многочлен Ньютона: рекурсивная функция

#include <iostream>
#include <math.h>
#include <fstream>
#define a  0
#define b  3
#define m  4 
#define n  4
#define J  20
#define h 0.15
using namespace std;
float RR(int k, int i, float *Xi,float *Yi)
{
    if (k == 1) return (Yi[i] - Yi[i+1])/(Xi[i]-Xi[i+1]);
    if (k == n-1) return 0;
    //тут надо добавить условия на случай выхода из массива 
    //я запутался здесь
    return ((RR(k-1,i,Xi,Yi) - RR(k-1,i+2,Xi,Yi)))/(Xi[i]-Xi[i+2]);
}
int main()
{
    ofstream v("E:\\lr3.txt", ios_base::out);
    if (!v.is_open() == NULL) 
    {
    float *N = new float [J];   
    float *Xi = new float [m];
    float *Yi = new float [m];
    int kkk = 1;
    for (int i =0;i<m;i++)
    {
        Xi[i] = a+(kkk-1)*(b-a)/(m-1);
        kkk++;
        Yi[i] = Xi[i]*Xi[i] - 10*pow(sin(Xi[i]),2);
        cout<<Xi[i]<<endl;
        cout<<Yi[i]<<endl;
    }
    float iter = 0;
    float *Xj = new float [J];
    for (int i = 0; i<J+1; i++)
    {
    Xj[i] = a+iter*(b-a)/20;
    cout<<Xj[i]<<endl;
    iter++;
    }
    float *Yj = new float [J];
    for (int i =0;i<J;i++)
    {
    Yj[i] = Xj[i]*Xj[i] - 10*pow(sin(Xj[i]),2);
    }   
    float u;
    for(int j = 0;j<J;j++)
        {
    float den=1;
            N[j] = Yi[0];
            int k = 0;
            //int fact = 1;
            for (int i = 0; i<n; i++)
                {
                //cout<<fact<<endl;
                den*=(Xj[j]-Xi[i]);
                    
                    N[j]+=den*RR(k,i,Xi,Yi);
                    k++;
                    //cout << fact*pow(h,i+1) << endl;
                    //(fact*pow(h,i+1))
                    //cout << RR(k, i, Xi, Yi) << endl;
                    //fact *= (i+2);
                }
            u = N[j] -Yj[j];
            cout<<"d(X"<<1+j<<") = "<<u<<endl;
            v<<"d(X"<<1+j<<") = "<<u<<endl;
        }
    }
    system("pause");
    return 0;
}