a concerned citizen - 1 year ago 86

C++ Question

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`

`long long`

Now, I made this program as a test, which is why the

`main()`

`gmpxx.h`

`f()`

`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`

`f()`

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=`

Answer Source

Gcc can point out one 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 is always true [-Wtype-limits]
for (unsigned int k=coeff.size(); k>=0; --k)
~^~~
```

Unsigned types should be avoided, unless you are implementing a bigint or a bitfield and know what you are doing.