include/algebramix/root_modular.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : root_modular.hpp
* DESCRIPTION: Roots for modulars over integers
* COPYRIGHT  : (C) 2010  Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_ROOT_MODULAR_HPP
#define __MMX_ROOT_MODULAR_HPP

#include <basix/list.hpp>
#include <numerix/modular_int.hpp>
#include <numerix/modular_integer.hpp>
#include <algebramix/modular_polynomial.hpp>


namespace mmx {

#define C          typename scalar_type_helper<Polynomial>::val
#define TMPL       template<typename Polynomial>
#define Modulus_polynomial modulus<Polynomial>
#define Modular_polynomial modular<Modulus_polynomial>

struct root_modular_naive {

TMPL static Polynomial
pow_mod (const Polynomial& a, const integer& n, const Polynomial& f) {
  Modulus_polynomial _f= Modular_polynomial::get_modulus ();
  Modular_polynomial::set_modulus (f);
  Polynomial g= * binpow (Modular_polynomial (a), n);
  Modular_polynomial::set_modulus (_f);
  return g;
}

TMPL static Polynomial
degree_one_factorization (const Polynomial& f, const integer& p) {
  // return the product of the linear factors
  Polynomial x (C(1), 1);
  Polynomial g= pow_mod (x, p, f) - x;
  return gcd (g, f);
}

TMPL static void
linear_splitting (Polynomial& g, const Polynomial& f, const integer& p) {
  VERIFY (f != 0, "bug");
  nat n= deg (f); g= f;
  if (n == 1) return;
  Polynomial x (C(1), 1);
  nat counter= 0;
  while (deg (g) <= 0 || (nat) deg (g) == n) {
    vector<C> _a (C (0), n);
    for (nat i= 0; i < n; i++) _a[i]= C(rand ()); // TODO << improve
    Polynomial a (_a), b;
    g= gcd (a, f);
    if (deg (g) > 0) return;
    g= p == 2 ? gcd (a - 1, f)
              : gcd (pow_mod (a, (p-1) / 2, f) - 1, f);
    counter++;
  }
  //mmout << "Number of equal degree splitting failures: " << counter << "\n";
  ASSERT (deg (g) > 0 && (nat) deg (g) < n, "cannot find a linear factor"); }

TMPL static vector<Polynomial>
linear_factorization (const Polynomial& f, const integer& p) {
  VERIFY (f != 0, "bug");
  list<Polynomial> todo (f);
  vector<Polynomial> factors= vec<Polynomial> ();
  if (deg (f) <= 0) return factors;
  while (!is_nil (todo)) {
    Polynomial h= car(todo), g; todo= cdr(todo);
    linear_splitting (g, h, p); 
    if ((nat) deg (g) == 1) factors << g;
    else todo= cons (g, todo);
    if (deg (g) != deg (h)) todo= cons (h / g, todo);
  }
  return factors; }

TMPL static vector<C>
roots (const Polynomial& f, const integer& p) {
  if (deg (f) <= 0) {
    return vector<C> ();
  }
  //mmout << f << "\n";
  Polynomial g= degree_one_factorization (f, p);
  //mmout << g << "\n";
  vector<Polynomial> lin_facts= linear_factorization (g, p);
  vector<C> b (C(), N(lin_facts));
  for (nat i= 0; i < N(lin_facts); i++)
    b[i]= - lin_facts[i][0] / lin_facts[i][1];
  return b; }
}; // struct modular_root_naive

#undef C
#undef TMPL
#undef Modulus_polynomial
#undef Modular_polynomial

#define Modulus modulus<C,V>
#define Modular modular<Modulus,W>

template<typename C, typename V, typename W> vector<Modular>
separable_roots (const Modular& a, nat r) {
  // return all the r th roots of a
  C p= * get_modulus (a);
  ASSERT (Modular (r) != 0, "wrong argument");
  polynomial<Modular> x_r (C(1),r);
  return root_modular_naive::roots (x_r - a, p);
}

template<typename C, typename V, typename W> Modular
separable_root (const Modular& a, nat r) {
  // return one r th roots of a
  // FIXME << to optimize
  if (a == 1) return a;
  vector<Modular> ans= separable_roots (a, r); 
  ASSERT (N (ans) != 0, "no root");
  return ans [0];
}

template<typename C, typename V, typename W> vector<Modular>
polynomial_modular_roots (const polynomial<Modular>& P) {
  // return all the roots of P
  ASSERT (N (P) != 0, "wrong argument");
  return root_modular_naive::roots (P, * get_modulus (P[0]));
}

template<typename C, typename V, typename W> Modular
polynomial_modular_root (const polynomial<Modular>& P) {
  // return one of the regular roots of P
  // FIXME << to optimize
  if (P[0] == Modular (0)) return Modular (0);
  //FIXME: should always compute one root
  vector<Modular> ans= polynomial_modular_roots (P);
  polynomial<Modular> dP = derive(P);
  ASSERT (N (ans) != 0, "no root");
  for (nat i= 0; i < N(ans); i++)
    if (eval (dP,ans[i]) != Modular (0)) return ans[i];
  ERROR ("no regular root");
}

template<typename C, typename V, typename W> inline Modular
pth_root (const Modular& a) {
  // p-th root of a
  return a;
}

// general case

template<typename C, typename V, typename W> vector<Modular>
nth_roots (const Modular& a, nat r) {
  // return all the r th roots of a
  C p= * get_modulus (a);
  if (Modular (r) == 0) return nth_roots (a, r / as<nat> (p));
  return separable_roots (a, r);
}

template<typename C, typename V, typename W> Modular
nth_root (const Modular& a, nat r) {
  if (a == 1) return a;
  vector<Modular> ans= nth_roots (a, r); 
  ASSERT (N (ans) != 0, "no root");
  return ans [0];
}
#undef Modulus
#undef Modular
} // namespace mmx    
#endif // __MMX_ROOT_MODULAR_HPP

---