ginac/polynomial/smod_helpers.h

   1 #ifndef GINAC_POLYNOMIAL_SMOD_HELPERS_H
   2 #define GINAC_POLYNOMIAL_SMOD_HELPERS_H
   3 #include <cln/integer.h>
   4 #include <cln/integer_io.h>
   5 #include "ex.h"
   6 #include "numeric.h"
   7 #include "debug.hpp"
   8
   9 namespace GiNaC
  10 {
  11
  12 /// Z -> Z_p (in the symmetric representation)
  13 static inline cln::cl_I smod(const cln::cl_I& a, long p)
  14 {
  15         const cln::cl_I p2 = cln::cl_I(p >> 1);
  16         const cln::cl_I m = cln::mod(a, p);
  17         const cln::cl_I m_p = m - cln::cl_I(p);
  18         const cln::cl_I ret = m > p2 ? m_p : m;
  19         return ret;
  20 }
  21
  22 static inline cln::cl_I recip(const cln::cl_I& a, long p_)
  23 {
  24         cln::cl_I p(p_);
  25         cln::cl_I u, v;
  26         const cln::cl_I g = xgcd(a, p, &u, &v);
  27         cln::cl_I ret = smod(u, p_);
  28         cln::cl_I chck = smod(a*ret, p_);
  29         bug_on(chck != 1, "miscomputed recip(" << a << " (mod " << p_ << "))");
  30         return ret;
  31
  32 }
  33
  34 static inline numeric recip(const numeric& a_, long p)
  35 {
  36         const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
  37         const cln::cl_I ret = recip(a, p);
  38         return numeric(ret);
  39 }
  40
  41 static inline cln::cl_I to_cl_I(const ex& e)
  42 {
  43         bug_on(!is_a<numeric>(e), "argument should be an integer");
  44         bug_on(!e.info(info_flags::integer),
  45                 "argument should be an integer");
  46         return cln::the<cln::cl_I>(ex_to<numeric>(e).to_cl_N());
  47 }
  48
  49 struct random_modint
  50 {
  51         typedef long value_type;
  52         const value_type p;
  53         const value_type p_2;
  54
  55         random_modint(const value_type& p_) : p(p_), p_2((p >> 1))
  56         { }
  57         value_type operator()() const
  58         {
  59                 do {
  60                         cln::cl_I tmp_ = cln::random_I(p);
  61                         value_type tmp = cln::cl_I_to_long(tmp_);
  62                         if (tmp > p_2)
  63                                 tmp -= p;
  64                         return tmp;
  65                 } while (true);
  66         }
  67
  68 };
  69
  70 } // namespace GiNaC
  71
  72 #endif // GINAC_POLYNOMIAL_SMOD_HELPERS_H