Msieve 1.41 Feedback

[Brian Gladman]

Originally Posted by jasonp

10metreh: I was insufficiently precise; the first thing the NFS code tries is rating the polynomial

The 'integration failed' messages don't have anything to do with these crashes. The crash occurs because you have to find the (floating point) roots of the algebraic polynomial in order to set up the numerical integration, and unfortunately it's impossible to build a polynomial rootfinder that is simultaneously

- capable of finding complex roots
- bulletproof regardless of the polynomial coefficients
- less than 1000 lines of code that started life as a crappy Fortran program

The code in the library uses Laguerre's method from Numerical Recipes and is pretty simple, but fails once in a while and produces nonsensical roots, which cause the integrator to crash. It seems to fail more often for SNFS polynomials. For bulletproof rootfinders, the best choice is the Jenkins-Traub model, which is 1000 lines of horrible C.

I don't claim that it is pretty but it did start off as Fortran :-)

Would double precision Jenkins-Traub followed by multiple precision root refinement be likely to work? Or are the distinct (non-multiple) roots too close together to be resolved in double precision?

[jasonp]

Originally Posted by Brian Gladman

I don't claim that it is pretty but it did start off as Fortran :-)

Would double precision Jenkins-Traub followed by multiple precision root refinement be likely to work? Or are the distinct (non-multiple) roots too close together to be resolved in double precision?

The roots are well separated, it's just that the iteration diverges (very fast!) for some input polynomials.

I didn't mean to disparage the work you did porting that terrible Fortran code, and even cleaned it up a little bit before giving up trying to make it better. It's just that I don't want to resort to something I can't understand if I can possibly help it.

[FactorEyes]

Originally Posted by jasonp

I don't either, so the rootfinder needs a bulletproof method for approximating complex roots to enough accuracy that Laguerre's method is guaranteed to converge. Perhaps a continuation algorithm is needed.

I'll look into Jenkins-Traub, but I can't find my resource for it -- I moved last year, and some stuff is still buried deep in boxes.

But the problem still remains: if I muddle my way through an implementation, will it be understandable enough that it won't be a black box to other developers (i.e.: you)? Another black box would defeat the purpose of having maintainable, comprehensible code.

[Brian Gladman]

Originally Posted by jasonp

The roots are well separated, it's just that the iteration diverges (very fast!) for some input polynomials.

I didn't mean to disparage the work you did porting that terrible Fortran code, and even cleaned it up a little bit before giving up trying to make it better. It's just that I don't want to resort to something I can't understand if I can possibly help it.

I didn't take any offence Jason - I agree that the code is terrible. If it was in a better state making it multiple precision would be relatively easy.

I didn't publish this code but if anyone wants to use it as a starting point for something better I will be happy to lete them have it.

[schickel]

Originally Posted by Brian Gladman

I didn't take any offence Jason - I agree that the code is terrible. If it was in a better state making it multiple precision would be relatively easy.

I didn't publish this code but if anyone wants to use it as a starting point for something better I will be happy to lete them have it.

If it won't blind me just by looking at, I'd like to at least take a look at it. I can't guarantee anything will happen.....

[Brian Gladman]

Here it is for those who are brave enough to take a look at it. It is supposed to be a translation of CACM algorithm 493 from Fortran to C but it may well contain conversion errors and really needs checking against the original Fortran. I would appreciate reports if anyone does find any such errors.

There is one subroutine where the logic is difficult to translate into structured code so I have left the original Fortran gotos in place. This is a challenge for those who feel that all gotos should be replaced by structured, transparent and efficient code!

Brian Gladman

[miklin]

At the link below is posted a C++ translation of the FORTRAN routine RPOLY.FOR, which is posted off the NETLIB site as TOMS/493.

http://www.akiti.ca/rpoly_ak1_Intro.html

http://math.fullerton.edu/mathews/n2...Bib_lnk_1.html

http://www.hvks.com/Numerical/winsolve.htm

[Robert Holmes]

Originally Posted by miklin

At the link below is posted a C++ translation of the FORTRAN routine RPOLY.FOR, which is posted off the NETLIB site as TOMS/493.

http://www.akiti.ca/rpoly_ak1_Intro.html

http://math.fullerton.edu/mathews/n2...Bib_lnk_1.html

http://www.hvks.com/Numerical/winsolve.htm

For the sake of diversity, here's another:

http://www.crbond.com/download/misc/rpoly.cpp

[WraithX]

Question... Jason mentioned needing a polynomial solver that can handle polynomials with complex roots. But the algorithm that Brian converted is 493. All the references I've seen online say that 493 is for polys with real coefficients and algorithm 419 is for polys with complex coefficients. And the links recently posted point to the rpoly (real polynomial) fortran code, and not the cpoly (complex polynomial) fortran code. Am I mixing these up? Will either of these work for Jason, or does he need specifically one or the other?

[jasonp]

The specialization in these methods is for polynomials whose coefficients are real numbers, but all the methods (Jenkins-Traub included) can find the complex roots of these polynomials that the code needs.

Some of these rewritings of TOMS493 are actually quite nice.

[Brian Gladman]

Originally Posted by jasonp

The specialization in these methods is for polynomials whose coefficients are real numbers, but all the methods (Jenkins-Traub included) can find the complex roots of these polynomials that the code needs.

Some of these rewritings of TOMS493 are actually quite nice.

I have tried several of these (and my own) on both the original test data set and on ill conditioned data sets and compared the results with those given by the original Fortran code. So far the results only agree to five or six decimal digits and are all different beyond this point even when compiled and run on the same machine using the same underlying FP processor. This is disappointing to say the least

The original Fortran code used quad precision but all the conversions I have tried use only double precision. So I am now in the process of reconverting the original Fortran in a way that sllows MPFR to be used to obtain multiple precision FP operations in C.