Recursion calculating Legendre Polynom not working (segmentation fault) [C++]

Question

I made a function in C++ that calculates the N-th Legendre polynomial and its derivative based on recursive formulas. It returns either zeros, bogus numbers or the "segmentation fault" error. I can't see what I did wrong. Care to help me out? Compiler is g++.

#include <iostream>
#include <cmath>
#include <iomanip>
#include <vector>
#include <functional>

using namespace std;

struct LegendreValues
{
  double function_value;
  double derivative_value;
};

LegendreValues calculateLegendrePolynomAndDerivative(const int&, const double&);
double calculateLegendrePolynom(const int&, const double&);
double calculateLegendreDerivative(const int&, const double&);

int main()
{
  int N;
  LegendreValues values;
  
  cout << "Unesite red Legendreovog polinoma: " << endl;
  cin >> N;
     
  cout << calculateLegendrePolynom(N, 0.5) << "\n" << calculateLegendreDerivative(N, 0.5) << endl;
}

LegendreValues calculateLegendrePolynomAndDerivative(const int& N, const double& x)
{
  LegendreValues results;
  if (N==0)
  {
    results.function_value = 1;
    results.derivative_value = 0;

  }
  
  else if (N==1)
  {
    results.function_value = x;
    results.derivative_value = 1;

  }

  else if (N>1)
  {
    results.function_value = (1/N) * ((2*N-1) * x * calculateLegendrePolynomAndDerivative(N-1, x).function_value - (N-1) * calculateLegendrePolynomAndDerivative(N-2, x).function_value);
    
    results.derivative_value = N/(1-x*x) * ( calculateLegendrePolynomAndDerivative(N-1, x).function_value - x * calculateLegendrePolynomAndDerivative(N, x).function_value);
    
  }

  return results;
}

double calculateLegendrePolynom(const int& N, const double& x)
{
  double fval;
  if (N==0)
  {
    fval =  1;
  }
  
  else if (N==1)
  {
    fval =  x;
  }
  
  else if (N>1)
  {
    fval = (1/N) * ((2*N-1) * x * calculateLegendrePolynom(N-1, x) - (N-1) * calculateLegendrePolynom(N-2, x));
  }
  return fval;
}

double calculateLegendreDerivative(const int& N, const double& x)
{
  double dfval;
  if (N==0)
  {
    dfval =  0;
  }
  
  else if (N==1)
  {
    dfval =  1;
  }
  
  else if (N>1)
  {
    dfval = N/(1-x*x) * ( calculateLegendreDerivative(N-1, x) - x * calculateLegendreDerivative(N, x));
  }
  return dfval;
}

The recursive formulas used: Legendre recursions

As you can see, I tried splitting the two things apart, and it still doesn't work.

Answer 1

The (1/N) versus (1.0/N) problem is easily fixed. Less so is the runaway expansion in function evaluations if you actually compute them by recursion without memoisation (storing previously-calculated values for immediate recall).

Consider N=6, say. This spawns calls for N=5 and N=4. The first spawns calls for N=4 (which spawns further lower calls) and N=3 (which will already have been calculated); the latter spawns calls for N=3 (already calculated several times over) and 2 (ditto).

So you make many extra unnecessary function evaluations and can run out of recursion stack.

It is simpler just to iteratively (rather than recursively) calculate for i=0,1,...,N. After all, you are going to have to work these out anyway. The code below uses a vector, but you could improve that further because you only actually need to store the previous two levels at each stage - a bit like a Fibonacci sequence.

#include <iostream>
#include <vector>
using namespace std;

struct LegendreValues
{
  double function_value;
  double derivative_value;
};

LegendreValues calculate( const int& N, const double& x );



int main()
{
  int N;
  double x;
  cout << "Enter the degree of polynomial (N >= 0): ";
  cin >> N;
  cout << "Enter the value of x ( 0 <= x < 1 ): ";
  cin >> x;
     
  LegendreValues values = calculate( N, x );
  cout << "Function:   " << values.function_value << '\n'
       << "Derivative: " << values.derivative_value << '\n';
}

LegendreValues calculate( const int& N, const double& x )
{
  vector<LegendreValues> store( 1 + N );
  for ( int i = 0; i <= N; i++ )
  {
     if ( i == 0 )
     {
        store[i] = { 1.0, 0.0 };
     }
     else if ( i == 1 )
     {
        store[i] = { x, 1.0 };
     }
     else
     {
        store[i].function_value   = ( 1.0 / i ) * ( ( 2 * i - 1 ) * x * store[i-1].function_value - ( i - 1 ) * store[i-2].function_value );
        store[i].derivative_value = i / ( 1.0 - x * x ) * ( store[i-1].function_value - x * store[i].function_value );
     }
  }
  return store[N];
}

Output:

Enter the degree of polynomial (N >= 0): 100
Enter the value of x ( 0 <= x < 1 ): 0.5
Function:   -0.060518
Derivative: -7.03317