ginac/polynomial/cra_garner.cpp

   1 /** @file cra_garner.cpp
   2  *
   3  *  Garner's algorithm. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2022 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