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

i trying root-finding algorithm based on durand-kerner , need polynomials of orders maybe 50, coefficients go beyond long double or long long, trying (for 1st time) gmp.

now, made program test, why main() (plus not advanced programmer), , works if don't involve gmp, think it's in implementation makes crash. codeblock's debugger points lines inside gmpxx.h , @ function f(), line sum +=.

here's code, minus polynomial coefficient generation (here 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};     (unsigned int i=1; i<=x; ++i)         z *= i;     return z; }  // 2^n mpz_class pwr2(const unsigned int &n) {     mpz_class z {1};     (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) {     (unsigned int i=0; i<=n; ++i)         a[i] = fact(n + i) / fact(n - i) / fact(i) / pwr2(i);     //return *a; }  // definition f(x), called every iteration cplx f(const std::vector<mpz_class> &coeff, const cplx &xterm) {     cplx sum {0, 0};     cplx factor {1, 0};     (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 current root's value cplx prod(const std::vector<cplx> &roots, const cplx &root) {     cplx product {1.0, 0};     (unsigned int k=0; k<roots.size(); ++k)     {         // skip if element of 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 coefficients, initial , 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 , maximum numbers of iterations.     unsigned int maxiter {200};     unsigned int iters {0};      (unsigned int k=0; k<arraysize; ++k)     {         coeffs[k] = k+1;         std::cout << coeffs[k] << " ";     }     std::cout << '\n';      // initialize seed roots     (unsigned int k=0; k<n; ++k)         initial[k] = pow(cplx(0.6, 0.8), k);     // temporary array hold next iteration's values     bool delta[n];     mpf_class tmp[n];     while (iters <= maxiter)     {         (unsigned int k=0; k<n; ++k)         {             roots[k] = initial[k] - f(coeffs, initial[k]) / prod(initial, initial[k]);         }         // calculate abs() average of temporary roots         bool test {1};         (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         (unsigned int k=0; k<n; ++k)             initial[k] = roots[k];         ++iters;     }     /// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     (unsigned short k=0; k<n; ++k)         std::cout << tmp[k] << " ";     std::cout << "\n\n";     (unsigned short k=0; k<n; ++k)         std::cout << initial[k] << "\n";     std::cout << iters << "\n";      return 0; } 

the std::cout line in f() function gets printed if it's inside loop, else won't, error must have connection, i'm afraid goes on head. can me error? gets me whole algorithm works fine standard types.

---

here's 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, example.

gcc can point out 1 important issue:

a.cc: in function ‘cplx f(const std::vector<__gmp_expr<__mpz_struct [1], __mpz_struct [1]> >&, const cplx&)’: a.cc:40:40: warning: comparison of unsigned expression >= 0 true [-wtype-limits]      (unsigned int k=coeff.size(); k>=0; --k)                                        ~^~~ 

unsigned types should avoided, unless creating own bigint type (or possibly using bitfield) , know doing.