1 /** @file cra_garner.cpp
3 * Garner's algorithm. */
6 * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
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.
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.
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
23 #include "cra_garner.h"
26 #include <cln/integer.h>
27 #include <cln/modinteger.h>
37 retract_symm(const cl_MI& x, const cl_modint_ring& R,
40 cl_I result = R->retract(x);
41 if (result > (modulus >> 1))
42 result = result - modulus;
47 compute_recips(vector<cl_MI>& dst,
48 const vector<cl_I>& moduli)
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);
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)
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]);
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]);
78 for (size_t j = k - 1 /* NOT k - 2 */; j-- != 0; )
79 tmp = tmp*moduli[j] + R->canonhom(dst[j]);
81 cl_MI next = (R->canonhom(residues[k]) - tmp)*recips[k-1];
82 dst[k] = retract_symm(next, R, moduli[k]);
87 mixed_radix_2_ordinary(const vector<cl_I>& mixed_radix_coeffs,
88 const vector<cl_I>& moduli)
90 size_t k = mixed_radix_coeffs.size() - 1;
91 cl_I u = mixed_radix_coeffs[k];
93 u = u*moduli[k] + mixed_radix_coeffs[k];
97 cl_I integer_cra(const vector<cl_I>& residues,
98 const vector<cl_I>& moduli)
100 if (unlikely(moduli.size() < 2))
101 throw std::invalid_argument("integer_cra: need at least 2 moduli");
103 vector<cl_MI> recips(moduli.size() - 1);
104 compute_recips(recips, moduli);
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);