Roots of Legendre Polynomials c++

Question

I'm writing a program to find the roots of nth order Legendre Polynomials using c++; my code is attached below:

double* legRoots(int n)
{
 double myRoots[n];
 double x, dx, Pi = atan2(1,1)*4;
 int iters = 0;
 double tolerance = 1e-20;
 double error = 10*tolerance;
 int maxIterations = 1000;

 for(int i = 1; i<=n; i++)
 {
  x = cos(Pi*(i-.25)/(n+.5));
  do
  {
   dx -= legDir(n,x)/legDif(n,x);
   x += dx;
   iters += 1;
   error = abs(dx);
  } while (error>tolerance && iters<maxIterations);
  myRoots[i-1] = x;
 }
 return myRoots;
}

Assuming the existence of functioning Legendre Polynomial and Legendre Polynomial derivative generating functions, which I do have but I thought that would make for unreadable walls of code text. This function is functioning in the sense that it's returning an array calculated values, but they're wildly off, outputting the following:

3.95253e-323
6.94492e-310
6.95268e-310
6.42285e-323
4.94066e-323
2.07355e-317

where an equivalent function I've written in Python gives the following:

[-0.90617985 -0.54064082  0.          0.54064082  0.90617985]

I was hoping another set of eyes could help me see what the issue in my C++ code is that's causing the values to be wildly off. I'm not doing anything different in my Python code that I'm doing in C++, so any help anyone could give on this is greatly appreciated, thanks. For reference, I'm mostly trying to emulate the method found on Rosetta code in regards to Gaussian Quadrature: http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature.

this "double myRoots[n];" is not legal C++. n must be const. Also, the for loop indexing goes from 1 to n and arrays in C and C++ are indexed from zero so the last loop would access unallocated data. — doug, Nov 23 '15 at 00:49
Ah I see, what I wanted to do there was initialize an array of size n, a similar iteration worked when I implemented a version with std::vector but I'll double check. — chuckFURD56, Nov 23 '15 at 00:56
A vector does initialize the doubles to 0.0. A double myRoots[n] would not unless it was static and, of course, it requires that n be a const unless the compiler has extensions to C++. — doug, Nov 23 '15 at 02:32

Mikael Kuisma · Answer 1 · 2015-11-23T00:58:21.370

2

You are returning an address to a temporary variable in stack

{
    double myRoots[n];
    ...
    return myRoots; // Not a safe thing to do
}

I suggest changing your function definition to

void legRoots(int n, double *myRoots)

omitting the return statement, and defining myroots before calling the function

double myRoots[10];
legRoots(10, myRoots);

Option 2 would be to allocate myRoots dynamically with new or malloc.

edited Nov 23 '15 at 00:58

answered Nov 23 '15 at 00:53

Mikael Kuisma

302
1
2
9

Well that did the trick, defining the array outside of the loop according to a # defined n made it work. I'll see if I can return the same array by changing the function type to void as per your suggestions, but thanks! – chuckFURD56 Nov 23 '15 at 01:00

Roots of Legendre Polynomials c++

1 Answers1