a concerned citizen - 1 year ago 112
C++ Question

# C++ Durand-Kerner root finding algorithm crashes with GMP, but not with long double

I am trying a root-finding algorithm based on Durand-Kerner and I need it for polynomials of orders maybe up to 50, with coefficients that go beyond

`long double`
or
`long long`
, so I am trying (for the 1st time) GMP.

Now, I made this program as a test, which is why the
`main()`
is as it is (plus I am not an advanced programmer), and it works if I don't involve GMP, so I think it's something in my implementation that makes it crash. With Codeblock's debugger it points to some lines inside
`gmpxx.h`
and at the function
`f()`
, the line with the
`sum +=`
.

Here's the code, minus the polynomial coefficient generation (here a simple loop):

``````#include <iostream>
#include <vector>   // std::vector
#include <complex>
#include <cmath>
#include <iomanip>  // std::setprecision
#include <gmpxx.h>

typedef std::complex<mpf_class> cplx;

// x!
mpz_class fact(const unsigned int &x)
{
mpz_class z {1};
for (unsigned int i=1; i<=x; ++i)
z *= i;
return z;
}

// 2^n
mpz_class pwr2(const unsigned int &n)
{
mpz_class z {1};
for (unsigned int i=0; i<n; ++i)
z *= 2;
return (n == 0 ? 1 : z);
}

void bessel(std::vector<mpz_class> &a, const unsigned int &N)
{
for (unsigned int i=0; i<=N; ++i)
a[i] = fact(N + i) / fact(N - i) / fact(i) / pwr2(i);
//return *a;
}

// Definition for f(x), will be called for every iteration
cplx f(const std::vector<mpz_class> &coeff, const cplx &xterm)
{
cplx sum {0, 0};
cplx factor {1, 0};
for (unsigned int k=coeff.size(); k>=0; --k)
{
sum += static_cast<cplx>(coeff[k]) * factor;
factor *= xterm;
std::cout<<sum<<'\n';
}
return sum;
}

// Denominator, product of differences. root is the current root's value
cplx prod(const std::vector<cplx> &roots, const cplx &root)
{
cplx product {1.0, 0};
for (unsigned int k=0; k<roots.size(); ++k)
{
// Skip if an element of the matrix == current root
if (roots[k] == root)
product = product;
else
product *= root - roots[k];
}
return product;
}

int main()
{
std::cout << "N=";
unsigned int N;
double eps {1e-10}; // relative error
std::cin >> N;
std::cout << std::setprecision(16);

// Declaring arrays for coefficients, initial and final roots
unsigned int arraySize {N + 1};
std::vector<mpz_class> coeffs;
std::vector<cplx> initial;
std::vector<cplx> roots;
coeffs.resize(arraySize);
initial.resize(arraySize);
roots.resize(arraySize);
// Keep track of number and maximum numbers of iterations.
unsigned int maxIter {200};
unsigned int iters {0};

for (unsigned int k=0; k<arraySize; ++k)
{
coeffs[k] = k+1;
std::cout << coeffs[k] << " ";
}
std::cout << '\n';

// Initialize the seed roots
for (unsigned int k=0; k<N; ++k)
initial[k] = pow(cplx(0.6, 0.8), k);
// Temporary array to hold next iteration's values
bool delta[N];
mpf_class tmp[N];
while (iters <= maxIter)
{
for (unsigned int k=0; k<N; ++k)
{
roots[k] = initial[k] - f(coeffs, initial[k]) / prod(initial, initial[k]);
}
// Calculate the abs() average of the temporary roots
bool test {1};
for (unsigned int k=0; k<N; ++k)
{
tmp[k] = fabs(initial[k] - roots[k]);
delta[k] = tmp[k] <= eps;
test *= delta[k];
//test *= fabs(initial[k] - roots[k]) <= eps;
//std::cout << tmp[k] << " ";
}
//std::cout << '\n';
// Check if eps has been reached
if (test)
break;
// if not, initial roots take current roots' value
for (unsigned int k=0; k<N; ++k)
initial[k] = roots[k];
++iters;
}
/// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (unsigned short k=0; k<N; ++k)
std::cout << tmp[k] << " ";
std::cout << "\n\n";
for (unsigned short k=0; k<N; ++k)
std::cout << initial[k] << "\n";
std::cout << iters << "\n";

return 0;
}
``````

The
`std::cout`
line in the
`f()`
function gets printed only if it's inside the loop, else it won't, so the error must have some connection, but I'm afraid it goes over my head. Can someone help me with the error? What gets me is that the whole algorithm works just fine with standard types.

Here's the callstack:

``````#0 0x7ffff7951e00   __gmpf_set_z() (/usr/lib/libgmp.so.10:??)
#1 0x404562 __gmp_set_expr<__mpz_struct [1]>(f=0x7fffffffe2f0, expr=...) (/usr/include/gmpxx.h:2114)
#2 0x403661 __gmp_expr<__mpf_struct [1], __mpf_struct [1]>::__gmp_expr<__mpz_struct [1], __mpz_struct [1]>(this=0x7fffffffe2f0, expr=...) (/usr/include/gmpxx.h:1844)
#3 0x401d5c f(coeff=std::vector of length 5, capacity 5 = {...}, xterm=...) (/home/xxx/Documents/cpp/Durand_Kerner_gmp/main.cpp:51)
#4 0x40250f main() (/home/xxx/Documents/cpp/Durand_Kerner_gmp/main.cpp:112)
``````

The program starts showing
`N=`
, input 4, for example.

``````a.cc: In function ‘cplx f(const std::vector<__gmp_expr<__mpz_struct [1], __mpz_struct [1]> >&, const cplx&)’: