CLN/ginac: factoring announce, questions
Hello, the following is concerned with the numtheory part of CLN. I have appended a patch to numtheory/cl_IF.h file and three new files for the same directory that would provide a factoring engine. The code is a much improved version of what I sent to Richy Kreckel a few days ago, and it explains itself. I have serious design questions on the interface, though. Shortly summed, you can factor any integer now within a minute whose second largest prime factor is less than ~10^11. I am not an expert, and I will not even think about implementing one of the more exotic algorithms (squfof, elliptic curves, all kind of sieves). I did this because I need ginac with a divisors() function. The people familiar with ginac/CLN design are asked to please answer QUESTION 1. The function cl_factor_rec() is the recursive main engine and will be called by factor() that is the interface to it. What should factor() return to the outside? We have - two lists (prime factors, exponents) - a list of pairs (p1,e1) (p2,e2) ... pari does return data like this - a simple list (p1, e1, p2, e2, p3, ...) QUESTION 2. How should the pairs/lists be implemented? We have - vectors provided by CLN (GV/SV where's the difference?) - STL containers An alternative would be to implement factor() outside CLN, ie, within ginac using lst, but I think there is interest having it in CLN. I just ask because I do not think the STL should be duplicated by GV/SV which, supposedly, was started when gcc did not have a working STL. I am no expert there, either, but I think CLN's memory scheme can be adapted to STL by simply(?) using a tailored allocator object. QUESTION 3. Should divisors() go into CLN or ginac? Many thanks for your time reading through to here. Please help---I am just starting C++ again after 8 years of hiatus, and design thinking is a bit unfamiliar at this time. Anyway, here is the beef: Sincerely, ralf --- cl_IF.h.orig Wed Jun 16 16:16:08 2004 +++ cl_IF.h Tue Jun 15 19:50:14 2004 @@ -3,6 +3,7 @@ #ifndef _CL_IF_H #define _CL_IF_H +#include <vector> #include "cln/number.h" #include "cln/integer.h" @@ -49,6 +50,15 @@ // Returns false if n is definitely composite, and then sets factor = some // nontrivial factor or 0. extern cl_boolean cl_miller_rabin_test (const cl_I& n, int count, cl_I* factor); + +// Pollard rho algorithm +extern cl_I cl_pollard_rho (const cl_I& n, long iters); + +// Pollard p-1 algorithm +extern cl_I cl_pollard_pmone (const cl_I& n); + +// Core recursive factoring routine, needs container +extern int cl_factor_rec (const cl_I& n, std::vector<cl_I>& v); } // namespace cln --- cl_IF_factor_rec.cc.orig Wed Jun 16 16:18:44 2004 +++ cl_IF_factor_rec.cc Tue Jun 15 19:10:12 2004 @@ -0,0 +1,109 @@ +// cl_factor_rec(). + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_IF.h" + + +// Implementation. + +#include "cln/integer.h" + +namespace cln { + +int cl_factor_rec(const cl_I& n, std::vector<cl_I>& v) +// recursive buildup of prime factors vector, part of code stolen from +// isprobprime.cc. +// prime factors of n are collected in v, 0 returned if factorization complete. +{ + // With a Miller-Rabin count = 50 the final error probability is + // 4^-50 < 10^-30. + int count = 50; + // Trial division (rules out 87% of all numbers quickly). + const uintL l = integer_length(n); + if (l <= 32) + { + uint32 nn = cl_I_to_UL(n); + if (nn <= cl_small_prime_table_limit) + { + // Table lookup. + uintL i = cl_small_prime_table_search(nn); + if (i < cl_small_prime_table_size + && ((unsigned int) cl_small_prime_table[i] == nn + || nn == 2)) + { v.push_back (n); return 0; } // end leaf of recursion + } + if ((nn % 2) == 0) + { v.push_back (2); return cl_factor_rec (n>>1, v); } + if ((nn = cl_trialdivision(nn,1,cl_small_prime_table_limit)) >0) + { v.push_back (nn); return cl_factor_rec (exquo (n,nn), v); } + + // For small n, only few Miller-Rabin tests are needed. + if (nn < 2000U) count = 1; // {2} + else if (nn < 1300000U) count = 2; // {2,3} + else if (nn < 25000000U) count = 3; // {2,3,5} + else if (nn < 3200000000U) count = 4; // {2,3,5,7} + } + else if (l <= 64) + { + uint32 nhi = cl_I_to_UL(ldb(n,cl_byte(32,32))); + uint32 nlo = cl_I_to_UL(ldb(n,cl_byte(32,0))); + if ((nlo % 2) == 0) + { v.push_back (2); return cl_factor_rec (n>>1, v); } + uint32 t; + if ((t = cl_trialdivision (nhi,nlo,1,cl_small_prime_table_limit)) >0) + { v.push_back (t); return cl_factor_rec (exquo (n,t), v); } + } + else + { + if (evenp (n)) + { v.push_back (2); return cl_factor_rec (n>>1, v); } + uint32 t; + if ((t = cl_trialdivision (n,1,cl_small_prime_table_limit)) >0) + { v.push_back (t); return cl_factor_rec (exquo (n,t), v); } + } + + // Handle perfect powers. + int m; + cl_I r; + if (rootp (n, m, &r)) + { + std::vector<cl_I> vv; + int ret = cl_factor_rec (r, vv); + for (int l=0; l<m; ++l) + copy (vv.begin(), vv.end(), v.begin()); + return ret; + } + + // After factoring out small factors: is remaining number prime? + if (cl_miller_rabin_test (n, count, NULL)) + { + v.push_back (n); + return 0; + } // end leaf of recursion + + // 30k iterations of Pollard rho finds most 8-digit factors. + // 150k iterations of Pollard rho finds most 10-digit factors. + cl_I f; + f = cl_pollard_rho (n, 150000L); + if (f != 1 && f != n) + { v.push_back (f); return cl_factor_rec (exquo (n,f), v); } + + // Pollard p-1 using cl_small_prime_table. Does not often find a factor + // but does not cost much, either. If it finds one, it may have + // up to 18-20 digits, occasionally. + f = cl_pollard_pmone (n); + if (f != 1 && f != n) + { v.push_back (f); return cl_factor_rec (exquo (n,f), v); } + + // Last step (for now): let Pollard rho run 'forever'. + // If you don't want this, return 1 without calling rho() again. + v.push_back (cl_pollard_rho (n, 0)); + return 0; // end leaf of recursion +} + + +} // namespace cln + --- cl_IF_pollard_rho.cc.orig Wed Jun 16 16:21:23 2004 +++ cl_IF_pollard_rho.cc Tue Jun 15 17:04:20 2004 @@ -0,0 +1,42 @@ +// cl_pollard_rho(). + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_IF.h" + + +// Implementation. + +#include "cln/integer.h" + +namespace cln { + +cl_I cl_pollard_rho (const cl_I& n, long iters) +// If N is composite, returns a factor of N, for iters>0 not necessarily <>N. +// Adaptation of an implementation by Richard Crandall (giantint). +// This gets slow when factors become larger than 10^10. iters=0 runs until +// factor is found. Using modinteger slows function down by 10 percent. +{ + cl_I g=1, t1=3, t2=3; + + for (long l=0; !iters || l<iters; ++l) + { + if (l%100 == 0) // removing this line seems to make no big difference + { + if ((g = gcd (g, n)) != 1) + break; + } + t1 = mod (square (t1) + 2, n); + t2 = mod (square (t2) + 2, n); + t2 = mod (square (t2) + 2, n); + g = mod (g * abs (t2-t1), n); + } + g = gcd (g, n); + return g; +} + +} // namespace cln + + --- cl_IF_pollard_pmone.cc.orig Wed Jun 16 16:22:45 2004 +++ cl_IF_pollard_pmone.cc Tue Jun 15 17:48:57 2004 @@ -0,0 +1,48 @@ +// cl_pollard_pmone(). + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_IF.h" + + +// Implementation. + +#include "cln/modinteger.h" + +namespace cln { + +cl_I cl_pollard_pmone (const cl_I& n) +// Pollard p-1 using cl_small_prime_table. Does not often find a factor +// but does not cost much, either. If it finds one, it may have +// up to 18-20 digits, occasionally. +// Returns a factor of N, not necessarily <>n. +// Adaptation of an implementation by Richard Crandall (giantint). +{ + cl_modint_ring R = find_modint_ring (n); + cl_MI g(R,cl_I(1)), t1(R,cl_I(3)); + cl_I t; + + for(long l=0; l<cl_small_prime_table_size; ++l) + { + long cnt = long((8.0*log(2.0))/log(1.0*cl_small_prime_table[l])); + if (cnt < 2) + cnt = 1; + for (long k=0; k<cnt; ++k) + t1 = expt_pos (t1, cl_small_prime_table[l]); + + g = g * (t1-1); + if (l%100==0) + { + if ((t = gcd (R->retract(g), n)) != 1) + break; + g = R->canonhom (t); + } + } + t = gcd (R->retract(g), n); + return t; +} + +} // namespace cln +
Dear Ralf, On Wed, 16 Jun 2004, Ralf Stephan wrote:
the following is concerned with the numtheory part of CLN. I have appended a patch to numtheory/cl_IF.h file and three new files for the same directory that would provide a factoring engine. The code is a much improved version of what I sent to Richy Kreckel a few days ago, and it explains itself. I have serious design questions on the interface, though.
Shortly summed, you can factor any integer now within a minute whose second largest prime factor is less than ~10^11. I am not an expert,
Oh, c'mon!
and I will not even think about implementing one of the more exotic algorithms (squfof, elliptic curves, all kind of sieves). I did this because I need ginac with a divisors() function.
The people familiar with ginac/CLN design are asked to please answer QUESTION 1. The function cl_factor_rec() is the recursive main engine and will be called by factor() that is the interface to it. What should factor() return to the outside? We have - two lists (prime factors, exponents) - a list of pairs (p1,e1) (p2,e2) ... pari does return data like this - a simple list (p1, e1, p2, e2, p3, ...)
QUESTION 2. How should the pairs/lists be implemented? We have - vectors provided by CLN (GV/SV where's the difference?) - STL containers An alternative would be to implement factor() outside CLN, ie, within ginac using lst, but I think there is interest having it in CLN. I just ask because I do not think the STL should be duplicated by GV/SV which, supposedly, was started when gcc did not have a working STL. I am no expert there, either, but I think CLN's memory scheme can be adapted to STL by simply(?) using a tailored allocator object.
QUESTION 3. Should divisors() go into CLN or ginac?
Many thanks for your time reading through to here. Please help---I am just starting C++ again after 8 years of hiatus, and design thinking is a bit unfamiliar at this time. Anyway, here is the beef:
Thank you very much for your keen resubmission. I was still pondering over your first patch which did have some, err, rough edges -- like that concat() function poking on the heap. There are still issues that we need to iron out, like const correctness and how to setup tests, but that can wait. I generally have a quite restrictive attitude towards extensions to CLN, admitting only state-of-the-art implementations. But in the case of integer factorization I do agree that code that finds all factors of order 10^10 is a very good start and should be included in CLN proper. (Yes, that answers your third question.) It may be helful for many and it may be improved in the future if people feel inclined. However, it is important to get the interface right, as you have noticed already. So, to question 1, what should factor() return? I'ld hate to see two lists. (I would rather have two vectors, since these are somewhat more rigidly coupled by nature.) I also dislike the alternating list of factors and exponents. It bears semantics invisible to the language/compiler. The list of pairs of factors and exponents is the most reasonable return value. And, I think, the most common. Not only Pari uses this notation, I think NTL does as well (for polynomial factors, that is). Question 2 gives me headaches. You've already noticed that CLN does not use std::list<>. I guess this is because the univariate polynomial code was written at a time when it was en vogue to write one's own vector/list/map/... Since a list of integer factors does not necessarily need to have anything to do with the coefficients of a univariate polynomial, I tend to suggest not to focus on what is already present it CLN. It is more important to focus on readabiliby and to keep extensibility in mind and this leads to STL containers. Now I do have two questions: What do you mean by adapting CLN's memory scheme to STL? The user-visible side of the bridge constitute perfectly copyable objects with value semantics and are hence suited for STL, right? Just have a look at the Bernoulli number cache in GiNaC. Second, is std::list< factor_exponent_pair_or_whatever > really a wise decision? Why not use a std::vector<> instead? After all it doesn't make too much sense removing/inserting elements therein at random locations. Methinks NTL returns vectors as well, but I need to check again. Best wishes -richy. PS: Besides throwing away lots of emails because I can't wade through all my fu^Hine SPAM I sometimes don't have enough time to answer email instantly. This is especilly true if the email is quite non-trivial. Well, sorry for the late answer. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Richard, thanks for the detailed reply. Const correctness is noted for my next submission (which will be to you), as well as testing. I understand from your comments that, while factor() will go into CLN, divisors() won't.
Now I do have two questions: What do you mean by adapting CLN's memory scheme to STL?
As far as I understand, CLN doesn't use new/delete and its objects are allocated using malloc_hook. OTOH, it may not be sufficient to malloc_hook sizeof(std::vector<...>) and give that to the caller since STL might itself allocate space using its own (default) allocator in the process. If that is really of concern, the solution could be to write an allocator (STL object) using malloc_hook for that and other cases. That was what I meant. Meanwhile, I favor moving the burden of space allocation to the caller via working on a referenced vector of pairs of cl_Is.
The user-visible side of the bridge constitute perfectly copyable objects with value semantics and are hence suited for STL, right?
CLN's objects are suited for STL, of course, since STL uses the object's interface. I was sceptical about the reverse, i.e., if STL deletes everything(!) it creates with its normal malloc. The user relies on garbage collection of an STL container with CLN objects that is returned by CLN. But the point is moot, anyway, with factor having the interface int factor (const cl_I& n, std::vector<std::pair<cl_I,cl_I> >& pv);
Just have a look at the Bernoulli number cache in GiNaC. Second, is std::list< factor_exponent_pair_or_whatever > really a wise decision?
I used the word 'list' not careful enough, sorry. I meant 'container'.
Well, sorry for the late answer.
I would have waited for it but nothing in CLN indicates who is its maintainer, so I was unsure if I mailed to the right address. Regards, ralf
Dear Ralf, On Thu, 17 Jun 2004, Ralf Stephan wrote:
thanks for the detailed reply. Const correctness is noted for my next submission (which will be to you), as well as testing. I understand from your comments that, while factor() will go into CLN, divisors() won't.
No, I didn't say that. What's the implementation of divisors anyways? You factorize, and then build all possible permutations of the factors, taking into account their multiplicities, right? Is there anything else one can do that I'm not aware of? Anyways, I see no reason why it shouldn't go into CLN, even if it is only a minor wrapper function for factor(). Is there any?
Now I do have two questions: What do you mean by adapting CLN's memory scheme to STL?
As far as I understand, CLN doesn't use new/delete and its objects are allocated using malloc_hook. OTOH, it may not be sufficient to malloc_hook sizeof(std::vector<...>) and give that to the caller since STL might itself allocate space using its own (default) allocator in the process. If that is really of concern, the solution could be to write an allocator (STL object) using malloc_hook for that and other cases. That was what I meant.
Meanwhile, I favor moving the burden of space allocation to the caller via working on a referenced vector of pairs of cl_Is.
Sorry, I still don't understand why you'ld have to fiddle with malloc_hook if all you're dealing with are the cl_I objects that have value semantics.
The user-visible side of the bridge constitute perfectly copyable objects with value semantics and are hence suited for STL, right?
CLN's objects are suited for STL, of course, since STL uses the object's interface. I was sceptical about the reverse, i.e., if STL deletes everything(!) it creates with its normal malloc. The user relies on garbage collection of an STL container with CLN objects that is returned by CLN.
I'm still confused. Are you worried about heap storage allocated by CLN and not freed early enough by the STL container's dtors? Even if the STL container dtors don't return everything to the system, they have to invoke their elements' dtors right away, thus ensuring that CLN doesn't leak, right? Please correct me if I misunderstand.
But the point is moot, anyway, with factor having the interface
int factor (const cl_I& n, std::vector<std::pair<cl_I,cl_I> >& pv);
We don't need a cl_I for the multiplicities - machine precision integers are more than enough. Also, I would recommend to declare a type representing factors and their multiplicity instead of using pair<cl_I,uintC>. Also, what's the integer return value, anyway? Wouldn't the prototype std::vector<somethingContaining_cl_I_and_unsigned> factor(const cl_I&); give us more natural syntax? I can't see much sense in C-like argument shuttling.
Just have a look at the Bernoulli number cache in GiNaC. Second, is std::list< factor_exponent_pair_or_whatever > really a wise decision?
I used the word 'list' not careful enough, sorry. I meant 'container'.
Ah, okay. Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
your comments that, while factor() will go into CLN, divisors() won't.
No, I didn't say that.
Ah, noted.
What's the implementation of divisors anyways? You factorize, and then build all possible permutations of the factors, taking into account their multiplicities, right? Is there anything else one can do that I'm not aware of?
Well, I don't claim to be able to understand pari innards but they build divisors up while factoring. As speed is dependent on the number of distinct prime factors, there is no big difference IMHO. And their code is a mess, though fast.
Anyways, I see no reason why it shouldn't go into CLN, even if it is only a minor wrapper function for factor(). Is there any?
I understood you were cautious to put anything new into CLN but you clarified that now.
int factor (const cl_I& n, std::vector<std::pair<cl_I,cl_I> >& pv);
We don't need a cl_I for the multiplicities - machine precision integers are more than enough. Also, I would recommend to declare a type representing factors and their multiplicity instead of using pair<cl_I,uintC>.
These natural changes were done already, I use now long for the exponent, and the type becomes (in cln/numtheory.h) // prime factor-exponent pair. A vector of these is written by factor(). typedef struct cl_factexp { cl_I fact; long exp; cl_factexp (const cl_I& i, long l) { fact = i; exp = l; } bool operator== (const struct cl_factexp& other) { return fact==other.fact && exp==other.exp; } bool operator!= (const struct cl_factexp& other) { return !(*this == other); } void *operator new (size_t size) { return malloc_hook (size); } void operator delete (void *p) { free_hook (p); } } cl_IFE; typedef std::vector<cl_IFE*> cl_IFEV;
Also, what's the integer return value, anyway?
I had meanwhile changed it to bool: false if factoring is incomplete---which can only happen if one changes the code---but pending the following discussion this will change again.
Wouldn't the prototype
std::vector<somethingContaining_cl_I_and_unsigned> factor(const cl_I&);
give us more natural syntax? I can't see much sense in C-like argument shuttling.
First, you want to have signed exponents, it costs only one bit. Then, who is going to delete the vector? The cl_Is are garbage collected but not the STL vector itself. Thanks, ralf
participants (2)
-
Ralf Stephan
-
Richard B. Kreckel