ginac/polynomial/smod_helpers.h

   1 /** @file smod_helpers.h
   2  *
   3  *  Utility functions for modular arithmetic. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #ifndef GINAC_POLYNOMIAL_SMOD_HELPERS_H
  24 #define GINAC_POLYNOMIAL_SMOD_HELPERS_H
  25
  26 #include "ex.h"
  27 #include "numeric.h"
  28 #include "debug.h"
  29
  30 #include <cln/integer.h>
  31 #include <cln/integer_io.h>
  32
  33 namespace GiNaC {
  34
  35 /// Z -> Z_p (in the symmetric representation)
  36 static inline cln::cl_I smod(const cln::cl_I& a, long p)
  37 {
  38         const cln::cl_I p2 = cln::cl_I(p >> 1);
  39         const cln::cl_I m = cln::mod(a, p);
  40         const cln::cl_I m_p = m - cln::cl_I(p);
  41         const cln::cl_I ret = m > p2 ? m_p : m;
  42         return ret;
  43 }
  44
  45 static inline cln::cl_I recip(const cln::cl_I& a, long p_)
  46 {
  47         cln::cl_I p(p_);
  48         cln::cl_I u, v;
  49         const cln::cl_I g = xgcd(a, p, &u, &v);
  50         cln::cl_I ret = smod(u, p_);
  51         cln::cl_I chck = smod(a*ret, p_);
  52         bug_on(chck != 1, "miscomputed recip(" << a << " (mod " << p_ << "))");
  53         return ret;
  54
  55 }
  56
  57 static inline numeric recip(const numeric& a_, long p)
  58 {
  59         const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
  60         const cln::cl_I ret = recip(a, p);
  61         return numeric(ret);
  62 }
  63
  64 static inline cln::cl_I to_cl_I(const ex& e)
  65 {
  66         bug_on(!is_a<numeric>(e), "argument should be an integer");
  67         bug_on(!e.info(info_flags::integer),
  68                 "argument should be an integer");
  69         return cln::the<cln::cl_I>(ex_to<numeric>(e).to_cl_N());
  70 }
  71
  72 struct random_modint
  73 {
  74         typedef long value_type;
  75         const value_type p;
  76         const value_type p_2;
  77
  78         random_modint(const value_type& p_) : p(p_), p_2((p >> 1))
  79         { }
  80         value_type operator()() const
  81         {
  82                 do {
  83                         cln::cl_I tmp_ = cln::random_I(p);
  84                         value_type tmp = cln::cl_I_to_long(tmp_);
  85                         if (tmp > p_2)
  86                                 tmp -= p;
  87                         return tmp;
  88                 } while (true);
  89         }
  90
  91 };
  92
  93 } // namespace GiNaC
  94
  95 #endif // GINAC_POLYNOMIAL_SMOD_HELPERS_H