ginac/polynomial/cra_garner.cpp

   1 #include <cln/integer.h>
   2 #include <cln/modinteger.h>
   3 #include <vector>
   4 #include <cstddef>
   5 #include "cra_garner.hpp"
   6 #include "compiler.h"
   7
   8 namespace cln
   9 {
  10 using std::vector;
  11 using std::size_t;
  12
  13 static cl_I
  14 retract_symm(const cl_MI& x, const cl_modint_ring& R,
  15              const cl_I& modulus)
  16 {
  17         cl_I result = R->retract(x);
  18         if (result > (modulus >> 1))
  19                 result = result - modulus;
  20         return result;
  21 }
  22
  23 static void
  24 compute_recips(vector<cl_MI>& dst,
  25                const vector<cl_I>& moduli)
  26 {
  27         for (size_t k = 1; k < moduli.size(); ++k) {
  28                 cl_modint_ring R = find_modint_ring(moduli[k]);
  29                 cl_MI product = R->canonhom(moduli[0]);
  30                 for (size_t i = 1; i < k; ++i)
  31                         product = product*moduli[i];
  32                 dst[k-1] = recip(product);
  33         }
  34 }
  35
  36 static void
  37 compute_mix_radix_coeffs(vector<cl_I>& dst,
  38                          const vector<cl_I>& residues,
  39                          const vector<cl_I>& moduli,
  40                          const vector<cl_MI>& recips)
  41 {
  42         dst[0] = residues[0];
  43
  44         do {
  45                 cl_modint_ring R = find_modint_ring(moduli[1]);
  46                 cl_MI tmp = R->canonhom(residues[0]);
  47                 cl_MI next = (R->canonhom(residues[1]) - tmp)*recips[0];
  48                 dst[1] = retract_symm(next, R, moduli[1]);
  49         } while (0);
  50
  51         for (size_t k = 2; k < residues.size(); ++k) {
  52                 cl_modint_ring R = find_modint_ring(moduli[k]);
  53                 cl_MI tmp = R->canonhom(dst[k-1]);
  54
  55                 for (size_t j = k - 1 /* NOT k - 2 */; j-- != 0; )
  56                         tmp = tmp*moduli[j] + R->canonhom(dst[j]);
  57
  58                 cl_MI next = (R->canonhom(residues[k]) - tmp)*recips[k-1];
  59                 dst[k] = retract_symm(next, R, moduli[k]);
  60         }
  61 }
  62
  63 static cl_I
  64 mixed_radix_2_ordinary(const vector<cl_I>& mixed_radix_coeffs,
  65                        const vector<cl_I>& moduli)
  66 {
  67         size_t k = mixed_radix_coeffs.size() - 1;
  68         cl_I u = mixed_radix_coeffs[k];
  69         for (; k-- != 0; )
  70                 u = u*moduli[k] + mixed_radix_coeffs[k];
  71         return u;
  72 }
  73
  74 cl_I integer_cra(const vector<cl_I>& residues,
  75                  const vector<cl_I>& moduli)
  76 {
  77         if (unlikely(moduli.size() < 2))
  78                 throw std::invalid_argument("integer_cra: need at least 2 moduli");
  79
  80         vector<cl_MI> recips(moduli.size() - 1);
  81         compute_recips(recips, moduli);
  82
  83         vector<cl_I> coeffs(moduli.size());
  84         compute_mix_radix_coeffs(coeffs, residues, moduli, recips);
  85         cl_I result = mixed_radix_2_ordinary(coeffs, moduli);
  86
  87         return result;
  88 }
  89
  90 } // namespace cln
  91