Happy New Year!
[ginac.git] / ginac / polynomial / cra_garner.cpp
1 /** @file cra_garner.cpp
2  *
3  *  Garner's algorithm. */
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 #include "cra_garner.h"
24 #include "compiler.h"
25
26 #include <cln/integer.h>
27 #include <cln/modinteger.h>
28 #include <cstddef>
29 #include <vector>
30
31 namespace cln {
32
33 using std::vector;
34 using std::size_t;
35
36 static cl_I 
37 retract_symm(const cl_MI& x, const cl_modint_ring& R,
38              const cl_I& modulus)
39 {
40         cl_I result = R->retract(x);
41         if (result > (modulus >> 1))
42                 result = result - modulus;
43         return result;
44 }
45
46 static void
47 compute_recips(vector<cl_MI>& dst, 
48                const vector<cl_I>& moduli)
49 {
50         for (size_t k = 1; k < moduli.size(); ++k) {
51                 cl_modint_ring R = find_modint_ring(moduli[k]);
52                 cl_MI product = R->canonhom(moduli[0]);
53                 for (size_t i = 1; i < k; ++i)
54                         product = product*moduli[i];
55                 dst[k-1] = recip(product);
56         }
57 }
58
59 static void
60 compute_mix_radix_coeffs(vector<cl_I>& dst,
61                          const vector<cl_I>& residues,
62                          const vector<cl_I>& moduli,
63                          const vector<cl_MI>& recips)
64 {
65         dst[0] = residues[0];
66
67         do {
68                 cl_modint_ring R = find_modint_ring(moduli[1]);
69                 cl_MI tmp = R->canonhom(residues[0]);
70                 cl_MI next = (R->canonhom(residues[1]) - tmp)*recips[0];
71                 dst[1] = retract_symm(next, R, moduli[1]);
72         } while (0);
73
74         for (size_t k = 2; k < residues.size(); ++k) {
75                 cl_modint_ring R = find_modint_ring(moduli[k]);
76                 cl_MI tmp = R->canonhom(dst[k-1]);
77
78                 for (size_t j = k - 1 /* NOT k - 2 */; j-- != 0; )
79                         tmp = tmp*moduli[j] + R->canonhom(dst[j]);
80
81                 cl_MI next = (R->canonhom(residues[k]) - tmp)*recips[k-1];
82                 dst[k] = retract_symm(next, R, moduli[k]);
83         }
84 }
85
86 static cl_I
87 mixed_radix_2_ordinary(const vector<cl_I>& mixed_radix_coeffs,
88                        const vector<cl_I>& moduli)
89 {
90         size_t k = mixed_radix_coeffs.size() - 1;
91         cl_I u = mixed_radix_coeffs[k];
92         for (; k-- != 0; )
93                 u = u*moduli[k] + mixed_radix_coeffs[k];
94         return u;
95 }
96
97 cl_I integer_cra(const vector<cl_I>& residues,
98                  const vector<cl_I>& moduli)
99 {
100         if (unlikely(moduli.size() < 2))
101                 throw std::invalid_argument("integer_cra: need at least 2 moduli");
102
103         vector<cl_MI> recips(moduli.size() - 1);
104         compute_recips(recips, moduli);
105
106         vector<cl_I> coeffs(moduli.size());
107         compute_mix_radix_coeffs(coeffs, residues, moduli, recips);
108         cl_I result = mixed_radix_2_ordinary(coeffs, moduli);
109
110         return result;
111 }
112
113 } // namespace cln