ginac/polynomial/primes_factory.h

   1 #ifndef GINAC_CHINREM_GCD_PRIMES_FACTORY_H
   2 #define GINAC_CHINREM_GCD_PRIMES_FACTORY_H
   3 #include <cln/integer.h>
   4 #include <cln/numtheory.h>
   5 #include <limits>
   6 #include "smod_helpers.h"
   7 #include "debug.h"
   8
   9 namespace GiNaC
  10 {
  11
  12 /**
  13  * Find a `big' prime p such that lc mod p != 0. Helper class used by modular
  14  * GCD algorithm.
  15  */
  16 class primes_factory
  17 {
  18 private:
  19         // These primes need to be large enough, so that the number of images
  20         // we need to reconstruct the GCD (in Z) is reasonable. On the other
  21         // hand, they should be as small as possible, so that operations on
  22         // coefficients are efficient. Practically this means we coefficients
  23         // should be native integers. (N.B.: as of now chinrem_gcd uses cl_I
  24         // or even numeric. Eventually this will be fixed).
  25         cln::cl_I last;
  26         // This ensures coefficients are immediate.
  27         static const int immediate_bits = 8*sizeof(void *) - __alignof__(void *);
  28         static const long opt_hint = (1L << (immediate_bits >> 1)) - 1;
  29 public:
  30         primes_factory()
  31         {
  32                 last = cln::nextprobprime(cln::cl_I(opt_hint));
  33         }
  34
  35         bool operator()(long& p, const cln::cl_I& lc)
  36         {
  37                 static const cln::cl_I maxval(std::numeric_limits<long>::max());
  38                 while (last < maxval) {
  39                         long p_ = cln::cl_I_to_long(last);
  40                         last = cln::nextprobprime(last + 1);
  41
  42                         if (!zerop(smod(lc, p_))) {
  43                                 p = p_;
  44                                 return true;
  45                         }
  46                 }
  47                 return false;
  48         }
  49
  50         bool has_primes() const
  51         {
  52                 static const cln::cl_I maxval(std::numeric_limits<long>::max());
  53                 return last < maxval;
  54         }
  55 };
  56
  57 } // namespace GiNaC
  58
  59 #endif /* GINAC_CHINREM_GCD_PRIMES_FACTORY_H */
  60