Happy New Year!
[ginac.git] / 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