sweeneychris / rpolyplusplus Goto Github PK

A modern c++ root finding algorithm based on the original Jenkins-Traub RPOLY software.

CMake 0.80% Shell 8.00% C++ 83.73% Makefile 2.05% C 0.05% Python 4.87% M4 0.51%

polynomials algorithm zeros roots root-finding jenkins-traub

rpolyplusplus's Introduction

RPOLY - A Polynomial Root-finding library

A three-stage algorithm for finding roots of polynomials with real coefficients as outlined in: "A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration" by Jenkins and Traub, SIAM 1970. Please note that this variant is different than the complex-coefficient version, and is estimated to be up to 4 times faster.

The algorithm works by computing shifts in so-called "K-polynomials" that reveal the roots. These shifts are applied in three stages: Zero-shifts, Fixed-shifts, and Variable-shift iterations. Roots are revealed as real roots or as a pair of complex conjugate roots. After a root (or pair of roots) is found, it is divided from the polynomial and the process is repeated.

Dependencies

Eigen3 library: http://eigen.tuxfamily.org/

This library is header-only so the installation is simple.

Building

Run the following commands from the root directory of RpolyPlusPlus.

mkdir build
cd build
cmake ..
make

This should build the library. Note that the unit tests are enabled by default. To build without the unit tests change the cmake line to:

cmake .. -DBUILD_TESTING=Off

If testing is enabled, you can run the unit test from the build directory with:

./bin/find_polynomial_roots_jenkins_traub_test

All unit tests should pass.

Questions

Contact Chris Sweeney at [email protected]

rpolyplusplus's People

Contributors

Stargazers

Watchers

Forkers

caomw davidcl fmschleif bellout sytang velurimithun petroleumcyberneticsgroup gillesdegottex luisf1997cr georgbrown forkbackups see3p00 advanpix baldvin-kovacs icodein oktengfei cumbenic relaxar

rpolyplusplus's Issues

Incorrect roots found

I have been testing the library on some polynomials with integer coefficients. On many of them incorrect roots are found. One of the examples:

Coefficients: [1,0,-1,0,1,1,1,1,1,1,1,0,-1,0,1]

Roots found:

  1.1910 + 0.6054i
  1.1910 - 0.6054i
 -1.2766 + 0.3311i
 -1.2766 - 0.3311i
 -0.9114 + 0.8715i
 -0.9114 - 0.8715i
 -0.2781 + 1.1229i
 -0.2781 - 1.1229i
  0.4173 + 0.9860i
  0.4173 - 0.9860i
  0.6739 + 0.3241i
  0.6739 - 0.3241i
  0.1838 + 0.2115i
  0.1838 - 0.2115i

True roots:

  1.2051 + 0.5795i
  1.2051 - 0.5795i
 -0.5464 + 0.8375i
 -0.5464 - 0.8375i
 -0.9793 + 0.2023i
 -0.9793 - 0.2023i
 -0.8508 + 0.5254i
 -0.8508 - 0.5254i
  0.5333 + 0.8459i
  0.5333 - 0.8459i
 -0.0358 + 0.9994i
 -0.0358 - 0.9994i
  0.6739 + 0.3241i
  0.6739 - 0.3241i

Roots are sorted by magnitude.

Not correctly solving certain quartics

Under some situations, I have found that the solver is outputting all zeros, below are some coefficients that caused it for me. I have tried putting some of these into wolfram alpha, which gives non-zero solutions.

I was using the solver to render a torus in a ray-tracing application for an assignment.

I am currently working around this issue by detecting the all-zeros condition, and retrying with slightly different numbers. As this is relatively rare, it hasn't been causing a noticable performance impact.

          1
   -198.995
    15006.6
    -508127
6.52025e+06

          1
   -200.246
    15131.5
    -511320
6.52025e+06

          1
   -202.193
    15327.2
    -516281
6.52025e+06

          1
   -202.205
    15328.4
    -516313
6.52025e+06

          1
   -202.215
    15329.4
    -516336
6.52025e+06

          1
   -202.268
    15334.8
    -516476
6.52025e+06

          1
   -202.289
    15336.8
    -516525
6.52025e+06

          1
    -202.35
      15343
    -516680
6.52025e+06

          1
   -202.389
      15347
    -516781
6.52025e+06

          1
   -202.388
    15346.8
    -516775
6.52025e+06

          1
   -202.427
    15350.8
    -516876
6.52025e+06

          1
   -202.427
    15350.8
    -516876
6.52025e+06

code producing the unwanted behaviour:

//find closest intersection with a ray originating at posn, and traveling in unit direction d
float Torus::intersect(glm::vec3 posn, glm::vec3 d)
{

	glm::vec3 p = posn - this->center;

        //checking for intersection with bounding sphere for optimization reasons
	float len = glm::length(posn);
	glm::vec3 finalPos = posn + d*len;

	if (glm::length(finalPos) > (r1 + d1))
	{
		return -1.0;
	}


	float a = glm::dot(d, d);
	float b = 2 * glm::dot(p, d);
	float c = glm::dot(p, p) + powf(d1, 2) - powf(r1, 2);
	float alpha = powf(d.x, 2) + powf(d.z, 2);
	float beta = 2 * (p.x*d.x + p.z*d.z);
	float gamma = powf(p.x, 2) + powf(p.z, 2);
	Eigen::VectorXd polynomial(5);
	polynomial[0] = powf(a, 2);
	polynomial[1] = 2 * a*b;
	polynomial[2] = 2 * a*c + powf(b, 2) - 4 * powf(d1, 2) *alpha;
	polynomial[3] = 2 * b*c - 4 * powf(d1, 2) *beta;
	polynomial[4] = powf(c, 2) - 4 * powf(d1, 2) *gamma;

	Eigen::VectorXd *realRoots = new Eigen::VectorXd();
	Eigen::VectorXd *compRoots = new Eigen::VectorXd();
	rpoly_plus_plus::FindPolynomialRootsJenkinsTraub(polynomial, realRoots, compRoots);

	float closest = -1.0;
	bool nonzero = 0;
	for (int i = 0; (*realRoots).size() > i; i++)
	{
		double d = (*realRoots).data()[i];
		double c = (*compRoots).data()[i];
		nonzero = nonzero || (d != 0 || c != 0);
		if (fabs(c) < 1e-7 && d > 1e-7 && (d < closest || closest == -1.0))
		{
			closest = d;
		}
	}
	if (!nonzero)
	{
		std::cout << polynomial << "\n" << "\n";
		//return(this->intersect(glm::vec3(posn.x+0.0001, posn.y, posn.z), d));
	}
	//std::cout<<closest<<"\n";
	return closest;
}

program output(set to display surface normals for debugging purposes), the pixel sized holes represent the locations where the solver returned zeros:

Clang Reports Use of Uninitialized Variable

Clang complains that the uninitialized variable 'polynomial_at_root' at line 603 of find_polynomial_roots_jenkins_traub.cc is assigned to 'prev_polynomial_at_root' at line 618, which is subsequently used uninitialized.

Segfault for specific case

Creating the 6th degree polynomials with coefficients as given at the end
of this message ends with a segmentation fault error. I don't remember
which of the 4 cases given crashes, but I think it was number 2 or 3.
Running RPolyPlusPlus with the commit from Aug 11, 2015 will
not result in the crash. This means the error only happens with the
newest commit (on Mar 9, 2016).

Hilmar

coeffs =
1
-359951
3.78335e+10
-1.37262e+15
5.26491e+18
0
0

coeffs =
1
-611771
8.23017e+10
-5.62821e+14
-4.18369e+18
0
0

coeffs =
1
-691197
1.27678e+11
-6.9045e+15
8.80928e+19
0
0

coeffs =
1
-280525
2.53443e+10
-7.45354e+14
-7.27967e+17
0
0

sweeneychris / rpolyplusplus Goto Github PK

rpolyplusplus's Introduction

RPOLY - A Polynomial Root-finding library

Dependencies

Building

Questions

rpolyplusplus's People

Contributors

Stargazers

Watchers

Forkers

rpolyplusplus's Issues

Incorrect roots found

Not correctly solving certain quartics

Clang Reports Use of Uninitialized Variable

Segfault for specific case

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent