factorization and stuff
Hi, I am re-sending this email because it looks like the server problems this morning have swallowed the original posting. Hi, I recently checked in code for polynomial factorization. The multivariate case and maybe algebraic extensions are in the making. Some notes on the code: - The code is absolutely not optimized. Don't look at the code, please! - Alternative methods are missing yet. It is doing the factorization with only one combination of methods and there are a lot of (well, infinitely many) use cases where the code performs amazingly poor. - The code is debugged and it should give the right answer for a proper expressions. But handling of non-proper expressions is still missing. Proper ~ Polynomial with symbols. I want to have the code feature-complete first (because I need it for some other stuff), only then will I try to optimize the code and build in additional algorithms. But if you find a bug in the factorization now please tell me! I also recently updated the todo list on the GiNaC website. Well, factorization will be done, but for the other two new points I have done nothing yet and there is always the chance that I will do nothing in the future. If you have interest in working on them, just do so (but maybe tell me in advance). In case you are planning to add other new features, I could announce them on the todo list as well. @Alexei: if you want to repair that GCD and is_polynomial business go ahead! That would be great! Regards, Jens
Hi! Out of curiosity I had a short look at the code. One question: what is factor_sqrfree(const ex&) actually doing? As far as I can tell, it s not squarefree factorization of univariate polynomials. Is it? (I know it's static and inside an anonymous namespace so I'm not supposed to use it, but still.) Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi! Jens Vollinga wrote:
I want to have the code feature-complete first (because I need it for some other stuff), only then will I try to optimize the code and build in additional algorithms. But if you find a bug in the factorization now please tell me!
Well, sometimes it fails to factor. In ginsh:
factor(expand((x+1)*(x+4))); 4+x^2+5*x
-richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi, Richard B. Kreckel schrieb:
Well, sometimes it fails to factor. In ginsh:
factor(expand((x+1)*(x+4))); 4+x^2+5*x
thanks! I fixed this bug. Your bug report is probably connected to your question:
Out of curiosity I had a short look at the code. One question: what is factor_sqrfree(const ex&) actually doing? As far as I can tell, it s not squarefree factorization of univariate polynomials. Is it? (I know it's static and inside an anonymous namespace so I'm not supposed to use it, but still.)
It is the modular square-free factorization, but it had a bug ... (the counter for the exponent of a possible factor was not always increased, so the modular sqrfree factorization of your polynomial was (x+1)^1 and not (x+1)^2; also a check for that exponent was missing at another place). Since you already looked at the code I should maybe comment on it with respect to cln: there is quite some stuff in factor.cpp that actually could be part of cln. It is mainly the type/class for modular univariate polynomials. If you compare the features of my class with the equivalent type in cln you can see why I needed this additional type. My plan was/is to finish first the code (including the multivariate case and additional algorithms for the univariate case) in GiNaC and to see exactly what features such a type needs, before I then ask the friendly cln developers to consider upgrading their type. :-) Regards, Jens
Hi! Jens Vollinga wrote:
Richard B. Kreckel schrieb:
Well, sometimes it fails to factor. In ginsh:
factor(expand((x+1)*(x+4))); 4+x^2+5*x
thanks! I fixed this bug.
Indeed. Thanks! Now, again in ginsh notation:
factor(expand((x+4+x^2-x^3+43*x^4)*(x+1-x^2-3*x^3+4*x^4))); Segmentation fault
Since you already looked at the code I should maybe comment on it with respect to cln: there is quite some stuff in factor.cpp that actually could be part of cln. It is mainly the type/class for modular univariate polynomials. If you compare the features of my class with the equivalent type in cln you can see why I needed this additional type.
Well, I don't see. :-( What is your point? Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi, Richard B. Kreckel schrieb:
Now, again in ginsh notation:
factor(expand((x+4+x^2-x^3+43*x^4)*(x+1-x^2-3*x^3+4*x^4))); Segmentation fault
thanks again!
Since you already looked at the code I should maybe comment on it with respect to cln: there is quite some stuff in factor.cpp that actually could be part of cln. It is mainly the type/class for modular univariate polynomials. If you compare the features of my class with the equivalent type in cln you can see why I needed this additional type.
Well, I don't see. :-( What is your point?
Well, maybe I should have put my statement differently. Of course, there is no ultimate reason why the cln class could not have been used from the beginning on in the factor code. Obviously, everthing can be done with the cln class and several additional functions that do the job of the missing methods! But I started out missing several things and then, not sure what the future algorithms might need in addition, wrote a convenient class. If all is done, maybe some of the methods could go into the cln class, probably not things like to_ex() or operator<(), but other things like unit_normal() or divide(). Regards, Jens
Dear Jens, Jens Vollinga wrote:
Richard B. Kreckel schrieb:
Now, again in ginsh notation:
factor(expand((x+4+x^2-x^3+43*x^4)*(x+1-x^2-3*x^3+4*x^4))); Segmentation fault
thanks again!
factor(expand((1-x+x^2-x^3)*x^2)); Internal error: statement in file ./float/division/cl_F_ceil1.cc, line 21 has been reached!! Please send the authors of the program a description how you produced
And thank *you* for the patch which indeed fixes that problem! I hope I'm not being too pesky, but that square-free factorization keeps attracting my mind. And now I am getting: this error! Sorry for being such a killjoy. Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi, Richard B. Kreckel schrieb:
I hope I'm not being too pesky, but that square-free factorization keeps attracting my mind. And now I am getting:
what is behind this curiosity? New feature in cln maybe? :-)
factor(expand((1-x+x^2-x^3)*x^2)); Internal error: statement in file ./float/division/cl_F_ceil1.cc, line 21 has been reached!! Please send the authors of the program a description how you produced this error!
It is no problem of the sqrfree function. Actually, you can help me with this bug. The problematic line is cl_I normmc = ceiling1(the<cl_F>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N()))); maxcoeff is 4 in this case. The line cl_I normmc = ceiling1(the<cl_I>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N()))); does work in this case (cl_F -> cl_I in the cast). Why?!? It is just too late to read the manual ...
Sorry for being such a killjoy.
Well, I have to seriously warn you: any more bugs and the multivariate factorization will be delayed even more! ;-) Jens
Good morning Jens, Jens Vollinga wrote:
Richard B. Kreckel schrieb:
I hope I'm not being too pesky, but that square-free factorization keeps attracting my mind. And now I am getting:
what is behind this curiosity? New feature in cln maybe? :-)
No, there's no code in the pipeline.
factor(expand((1-x+x^2-x^3)*x^2)); Internal error: statement in file ./float/division/cl_F_ceil1.cc, line 21 has been reached!! Please send the authors of the program a description how you produced this error!
It is no problem of the sqrfree function. Actually, you can help me with this bug. The problematic line is
cl_I normmc = ceiling1(the<cl_F>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
maxcoeff is 4 in this case.
The line
cl_I normmc = ceiling1(the<cl_I>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
does work in this case (cl_F -> cl_I in the cast).
Why?!? It is just too late to read the manual ...
Each of the classes `cl_R', `cl_F', `cl_SF', `cl_FF', `cl_DF', `cl_LF' defines `TYPE sqrt (const TYPE& x)'. `x' must be >= 0. This function returns the square root of `x', normalized to be >= 0. If `x' is the square of a rational number, `sqrt(x)' will be a rational number, else it will return a floating-point approximation. So, cln::sqrt(cln::cl_N(4)) returns cl_N(2), an exact integer, but of static type cl_N. Casting that to cl_F, a float triggers the error. (You can't use the for constructing floating-point numbers, you should use the contructors.) But in your case, converting the<cl_R>(...) will convert to static type cl_R, which can at runtime be either integer or rational or floating-point, but not complex. Indeed, it seems to fix things.
Sorry for being such a killjoy.
Well, I have to seriously warn you: any more bugs and the multivariate factorization will be delayed even more! ;-)
I'm looking forward to testing the multivariate case, too! :-) nice day -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Jens, A couple of minuts ago, I wrote:
Sorry for being such a killjoy.
Well, I have to seriously warn you: any more bugs and the multivariate factorization will be delayed even more! ;-)
I'm looking forward to testing the multivariate case, too! :-)
nice day -richy.
Ah, I haven't had my morning coffee yet. So I forgot to include the bug of the day:
factor(x); Segmentation fault
:-) fun -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Jens, On 2008-08-21, I wrote:
Ah, I haven't had my morning coffee yet. So I forgot to include the bug of the day:
factor(x); Segmentation fault
That's fixed now. Thank you! Now, I'm getting
factor(expand((x+64*x^4)*(3*x+92*x^4+5*x^2))); Segmentation fault
Is this in any way related with this issue?
factor(x^37+1); Not a 32-bit integer: 137438953472
Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Jens, Richard B. Kreckel wrote:
On 2008-08-21, I wrote:
Ah, I haven't had my morning coffee yet. So I forgot to include the bug of the day:
factor(x); Segmentation fault
That's fixed now. Thank you!
Now, I'm getting
factor(expand((x+64*x^4)*(3*x+92*x^4+5*x^2))); Segmentation fault
Is this in any way related with this issue?
factor(x^37+1); Not a 32-bit integer: 137438953472
That's fixed now. Thank you! Now, I'm getting:
factor(expand((77+102*x^4+27*x+11*x^3+25*x^2)*(85+66*x^4+29*x+57*x^3+92*x^2))); Not a 32-bit integer: 2955677175 Woops! Sorry for being a nuisance. ;-) Well, it's easily fixed with this patch: diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 886819a..9081052 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -124,8 +124,8 @@ struct UniPoly // assert: poly is in Z[x] Term t; for ( int i=poly.degree(x); i>=poly.ldegree(x); --i ) { - int coeff = ex_to<numeric>(poly.coeff(x,i)).to_int(); - if ( coeff ) { + cl_I coeff = the<cl_I>(ex_to<numeric>(poly.coeff(x,i)).to_cl_N()); + if ( !zerop(coeff) ) { t.c = R->canonhom(coeff); if ( !zerop(t.c) ) { t.exp = i; This leaves me wondering whether other numeric -> int conversions should be treated likewise. What about gamma inside the hensel_univar function? I don't know what it is. Can we assume with reasonable certainty that it fits the target unsigned int format? Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Hi, Richard B. Kreckel schrieb:
Now, I'm getting:
factor(expand((77+102*x^4+27*x+11*x^3+25*x^2)*(85+66*x^4+29*x+57*x^3+92*x^2))); Not a 32-bit integer: 2955677175
Woops! Sorry for being a nuisance. ;-)
thanks for the report, and thanks for the patch!
This leaves me wondering whether other numeric -> int conversions should be treated likewise. What about gamma inside the hensel_univar function? I don't know what it is. Can we assume with reasonable certainty that it fits the target unsigned int format?
You are completely right. This is bad stuff and it is probably happening at more places. When I started programming the code I didn't fully understand the algorithm behind it. Programming it was actually the way to understand how the mathematics works. Anyway, now the code is feature complete. Next is the actual testing and debugging. Right, I didn't really test the code yet! If that second step is finished I am going to clean-up, rewrite, comment, and optimize the code. Currently, the multivariate code breaks for example for x+x*y-1. Have to fix that. Regards, Jens
participants (2)
-
Jens Vollinga
-
Richard B. Kreckel