bivpol.cc 3.31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
/* Common header file for the rrspace software

This file is part of the rrspace project.

This library is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free
Software Foundation; either version 2.1 of the License, or (at your option) any
later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
details.

You should have received a copy of the GNU Lesser General Public License along
with this library; if not, write to the Free Software Foundation, Inc., 51
Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

*/

#include "bivpol.h"

using namespace std;
using namespace NTL;

size_t 
BivPol::degree() const {
  size_t res = 0;
  for (size_t i = 0; i < coeffs.size(); ++i)
    res = max(res, i + max(deg(coeffs[i]), 0)); 
  return res;
}

bool 
BivPol::IsZero() const {
  for (size_t i =0; i < coeffs.size(); ++i)
    if (!NTL::IsZero(coeffs[i]))
      return false;
  return true;
}

ZZ_pX 
BivPol::mod_eval(const ZZ_pX& g, const ZZ_pX& modulus) const {
  assert(modulus != ZZ_pX(0));
  if (deg(modulus) == 0) 
   return ZZ_pX(0); 
  ZZ_pXModulus mod(modulus);
  ZZ_pX power_g(1);
  ZZ_pX res(0);
  for (size_t i = 0; i < coeffs.size(); ++i) {
    res += MulMod(coeffs[i] % mod, power_g, mod);
    MulMod(power_g, power_g, g, mod);
  }
  return res;
}
  
BivPol 
BivPol::diffY() const {
  std::vector<ZZ_pX> coeffs_res;
  for (size_t i = 1; i < coeffs.size(); ++i)
    coeffs_res.push_back(i*coeffs[i]);
  return BivPol(coeffs_res);
}

ZZ_pX
NewtonHenselStep(const BivPol& q, const ZZ_pX& s, const ZZ_pX& modulus) {
  assert(deg(modulus) > 0);
  assert(IsZero(q.mod_eval(s, modulus)));
  ZZ_pX mod2 = sqr(modulus);
  ZZ_pX denom    = q.diffY().mod_eval(s, mod2);
  ZZ_pX num      = q.mod_eval(s, mod2);
  assert(deg(GCD(denom, mod2)) == 0);
  ZZ_pX invdenom = InvMod(denom, mod2);
  ZZ_pX NH_step  = MulMod(num, invdenom, mod2);
  ZZ_pX res      = s - NH_step;
  assert(IsZero(q.mod_eval(res, mod2)));
  return res;
}

ZZ_pX
NewtonHensel(const BivPol& q, const ZZ_pX& s, const ZZ_pX& modulus, std::size_t k) {
  ZZ_pX R    = s;
  ZZ_pX mod2 = modulus;
  for (std::size_t i = 0; i < k; ++i) {
    R = NewtonHenselStep(q, R, mod2);
    mod2 = sqr(mod2);
  }
  return R;
}

ostream& 
operator<<(ostream& o, const BivPol& P) {
  o << "[ ";
  for (size_t i = 0; i <= P.degree_y(); ++i)
    o << P.coeff(i) << " ";
  o << "]";
  return o;
}

ostream& 
PrintMagma(ostream& o, const ZZ_pX& P) {
  o << "(";
  for (size_t i = 0; i < (size_t)deg(P); ++i)
    if (!IsZero(coeff(P, i)))
      o << coeff(P, i) << "*x^" << i << " + ";
  o << coeff(P, deg(P)) << "*x^" << deg(P);
  o << ")";
  return o;
}

ostream& 
PrintMagma(ostream& o, const BivPol& P) {
  for (size_t i = 0; i < P.degree_y(); ++i)
    if (!IsZero(P.coeff(i)))
      PrintMagma(o, P.coeff(i)) << "*y^" << i << " + ";
  PrintMagma(o, P.coeff(P.degree_y())) << "*y^" << P.degree_y();
  return o;
}


istream& 
operator>>(istream& i, BivPol& res) {
  string s = "";
  size_t deg;
  res.coeffs.clear();
  ZZ_pX tmp;
  i >> s;
  i >> deg;
  for (size_t j = 0; j <= deg; ++j) {
    i >> tmp;
    res.coeffs.push_back(tmp);
  }
  i >> s;
  return i;
}