/Users/mourrain/Devel/mmx/algebramix/include/algebramix/algebraic_extension.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : algebraic_extension.hpp
* DESCRIPTION: Algebraic extensions based on a primitive element
* COPYRIGHT  : (C) 2011  Joris van der Hoeven
*******************************************************************************
* 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_ALGEBRAIC_EXTENSION_HPP
#define __MMX_ALGEBRAIC_EXTENSION_HPP
#include <algebramix/polynomial.hpp>
#include <algebramix/matrix.hpp>
#include <basix/symbol.hpp>
namespace mmx {
#define TMPL_DEF template<typename C>
#define TMPL template<typename C>
#define Polynomial polynomial<C>
#define Extension algebraic_extension<C>
#define Element typename algebraic_extension<C>::El
TMPL class algebraic_extension;

/******************************************************************************
* The algebraic_extension type
******************************************************************************/

TMPL_DEF
class algebraic_extension {
MMX_ALLOCATORS;
public:
  Polynomial mp;          // minimal polynomial
  typedef Polynomial El;  // representation of elements

public:
  inline algebraic_extension ():
    mp () {}
  inline algebraic_extension (const format<C>& fm):
    mp (vec<C> (promote (0, fm), promote (1, fm))) {}
  inline algebraic_extension (const Polynomial& mp2):
    mp (mp2) {}
  template<typename C2> inline
  algebraic_extension (const algebraic_extension<C2>& ext2):
    mp (as<Polynomial > (ext2.mp)) {}
  inline Polynomial operator * () const { return mp; }
};

/******************************************************************************
* Basic operations
******************************************************************************/

TMPL inline format<C> CF (const Extension& x) { return CF(*x); }
TMPL inline nat hash (const Extension& x) { return hash (*x); }
TMPL inline nat exact_hash (const Extension& x) { return exact_hash (*x); }
TMPL inline nat hard_hash (const Extension& x) { return hard_hash (*x); }
TMPL inline bool operator == (const Extension& x, const Extension& y) {
  return (*x) == (*y); }
TMPL inline bool operator != (const Extension& x, const Extension& y) {
  return (*x) != (*y); }
TMPL inline bool exact_eq (const Extension& x, const Extension& y) {
  return exact_eq (*x, *y); }
TMPL inline bool exact_neq (const Extension& x, const Extension& y) {
  return exact_neq (*x, *y); }
TMPL inline bool hard_eq (const Extension& x, const Extension& y) {
  return hard_eq (*x, *y); }
TMPL inline bool hard_neq (const Extension& x, const Extension& y) {
  return hard_neq (*x, *y); }

TMPL inline syntactic flatten (const Extension& x) {
  return syn ("Extension", flatten (x.mp)); }

TMPL
struct binary_helper<Extension >: public void_binary_helper<Extension > {
  typedef Polynomial R;
  typedef binary_helper<R> H;
  static inline string short_type_name () {
    return "Algext" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Algebraic_extension", Full_type_name (C)); }
  static inline generic disassemble (const Extension& x) {
    return binary_disassemble<R> (*x); }
  static inline Extension assemble (const generic& x) {
    return Extension (binary_assemble<R> (x)); }
  static inline void write (const port& out, const Extension& x) {
    return H::write (out, *x); }
  static inline Extension read (const port& in) {
    return H::read (in); }
};

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Element
square (const Extension& ext, const Element& p1) {
  return rem (square (p1), ext.mp);
}

TMPL inline Element
mul (const Extension& ext, const Element& p1, const Element& p2) {
  return rem (p1 * p2, ext.mp);
}

TMPL Element
invert (const Extension& ext, const Element& p1) {
  Element a= promote (0, p1), b= promote (0, p1), c= gcd (p1, ext.mp, a, b);
  return a;
}

TMPL Element
div (const Extension& ext, const C& c1, const Element& p2) {
  Element a= promote (0, p2), b= promote (0, p2), c= gcd (p2, ext.mp, a, b);
  return c1 * a;
}

TMPL Element
div (const Extension& ext, const Element& p1, const Element& p2) {
  Element a= promote (0, p2), b= promote (0, p2), c= gcd (p2, ext.mp, a, b);
  return rem (p1 * a, ext.mp);
}

TMPL inline bool
is_zero (const Extension& ext, const Element& p1) {
  return deg (p1) < 0;
}

TMPL inline int
sign (const Extension& ext, const Element& p1) {
  if (deg (p1) <= 0) return sign (p1[0]);
  ERROR ("cannot compute sign");
}

/******************************************************************************
* Joining algebraic extensions
******************************************************************************/

TMPL vector<C>
shift1 (const Extension& ext1, const Extension& ext2, const vector<C>& v) {
  // multiply the bivariate polynomial represented by v with x1
  nat d1= deg (ext1.mp), d2= deg (ext2.mp);
  vector<C> r= fill<C> (promote (0, CF(ext1)), d1*d2);
  for (nat i2=0; i2<d2; i2++) {
    vector<C> c= fill<C> (promote (0, CF(ext1)), d1);
    for (nat i1=0; i1<d1; i1++)
      c[i1]= v[i1*d2 + i2];
    Element p= rem (lshiftz (Element (c), 1), ext1.mp);
    for (nat i1=0; i1<d1; i1++)
      r[i1*d2 + i2]= p[i1];
  }
  return r;
}

TMPL vector<C>
shift2 (const Extension& ext1, const Extension& ext2, const vector<C>& v) {
  // multiply the bivariate polynomial represented by v with x2
  nat d1= deg (ext1.mp), d2= deg (ext2.mp);
  vector<C> r= fill<C> (promote (0, CF(ext1)), d1*d2);
  for (nat i1=0; i1<d1; i1++) {
    vector<C> c= fill<C> (promote (0, CF(ext1)), d2);
    for (nat i2=0; i2<d2; i2++)
      c[i2]= v[i1*d2 + i2];
    Element p= rem (lshiftz (Element (c), 1), ext2.mp);
    for (nat i2=0; i2<d2; i2++)
      r[i1*d2 + i2]= p[i2];
  }
  return r;
}

TMPL matrix<C>
mul_matrix (const Extension& ext1, const Extension& ext2, const vector<C>& v) {
  // let p (x1, x2) be the bivariate polynomial represented by v
  // return the matrix whose rows represent polynomials p (x1, x2) x1^i x2^j
  nat d1= deg (ext1.mp), d2= deg (ext2.mp);
  matrix<C> r (promote (0, CF(ext1)), d1*d2, d1*d2);
  vector<C> aux= fill<C> (promote (0, CF(ext1)), d1*d2);
  for (nat i1=0; i1<d1; i1++) {
    if (i1 == 0)
      for (nat k=0; k<d1*d2; k++)
        r (0, k)= v[k];
    else {
      for (nat k=0; k<d1*d2; k++)
        aux[k]= r((i1-1)*d2, k);
      aux= shift1 (ext1, ext2, aux);
      for (nat k=0; k<d1*d2; k++)
        r(i1*d2, k)= aux[k];
    }
    for (nat i2=1; i2<d2; i2++) {
      for (nat k=0; k<d1*d2; k++)
        aux[k]= r(i1*d2 + i2-1, k);
      aux= shift2 (ext1, ext2, aux);
      for (nat k=0; k<d1*d2; k++)
        r(i1*d2 + i2, k)= aux[k];
    }
  }
  return r;
}

TMPL matrix<C>
pow_matrix (const Extension& ext1, const Extension& ext2, const vector<C>& v) {
  // let p (x1, x2) be the bivariate polynomial represented by v
  // return the matrix whose rows represent the powers p (x1, x2)^i
  nat d1= deg (ext1.mp), d2= deg (ext2.mp);
  matrix<C> m= mul_matrix (ext1, ext2, v);
  matrix<C> r (promote (0, CF(ext1)), d1*d2+1, d1*d2);
  vector<C> aux= fill<C> (promote (0, CF(ext1)), d1*d2);
  aux[0]= promote (1, CF(ext1));
  for (nat i1=0; i1<=d1*d2; i1++) {
    for (nat i2=0; i2<d1*d2; i2++)
      r (i1, i2)= aux[i2];
    aux= aux * m;
  }
  return r;
}

TMPL matrix<C>
pow_matrix (const Extension& ext1, const Extension& ext2) {
  // return matrix whose rows represent the powers of a primitive element
  nat d1= deg (ext1.mp), d2= deg (ext2.mp);
  vector<C> v= fill<C> (promote (0, CF(ext1)), d1*d2);
  for (nat i=1; i<1000; i++) {
    v[1]= promote (1, CF(ext1)); v[d2]= promote ((int) i, CF(ext1));
    matrix<C> m= pow_matrix (ext1, ext2, v);
    if (rank (m) == d1*d2) return m;
  }
  ERROR ("unexpected situation");
}

template<typename C> polynomial<C>
square_free (const polynomial<C>& p) {
  polynomial<C> g= gcd (p, derive (p));
  if (deg (g) == 1) return p;
  else return quo (p, g);
}

TMPL Extension
join (const Extension& ext1, const Extension& ext2, Element& z1, Element& z2) {
  // Return the smallest common extension ext of ext1 and ext2
  // On exit, z1 and z2 contains the primitive els of ext1 and ext2 inside ext
  ASSERT (N(ext1.mp) > 0, "uninitialized algebraic extension");
  ASSERT (N(ext2.mp) > 0, "uninitialized algebraic extension");
  if (deg (ext1.mp) == 1) {
    z1= Polynomial (-ext1.mp[0]/ext1.mp[1]);
    z2= Polynomial (promote (1, CF(ext1)), (nat) 1);
    return ext2;
  }
  else if (deg (ext2.mp) == 1) {
    z1= Polynomial (promote (1, CF(ext1)), (nat) 1);
    z2= Polynomial (-ext2.mp[0]/ext2.mp[1]);
    return ext1;
  }
  else if (hard_eq (ext1, ext2)) {
    z1= Polynomial (promote (1, CF(ext1)), (nat) 1);
    z2= Polynomial (promote (1, CF(ext1)), (nat) 1);
    return ext1;
  }
  else {
    nat n= deg (ext1.mp) * deg (ext2.mp);
    matrix<C> m= transpose (pow_matrix (ext1, ext2));
    matrix<C> u= invert (range (m, 0, 0, n, n));
    vector<C> v= - (u * column (m, n));
    v << promote (1, CF(ext1));
    Polynomial mp= Polynomial (v);
    Polynomial sf= square_free (mp);
    Extension ext= Extension (sf);
    v= fill<C> (promote (0, CF(ext1)), n);
    v[deg(ext2.mp)]= promote (1, CF(ext1));
    z1= Element (u * v);
    if (deg (sf) < deg (mp)) z1= rem (z1, sf);
    v= fill<C> (promote (0, CF(ext1)), n);
    v[1]= promote (1, CF(ext1));
    z2= Element (u * v);
    if (deg (sf) < deg (mp)) z2= rem (z2, sf);
    return ext;
  }
}

/******************************************************************************
* Upgrading to larger algebraic extensions
******************************************************************************/

TMPL Element
compose (const Extension& ext, const Polynomial& p, const Element& q) {
  // Evaluate p at q in a given extension
  Element sum= promote (0, q);
  for (int i= deg (p); i>=0; i--) {
    if (i != deg (p)) sum= rem (sum * q, ext.mp);
    sum += p[i];
  }
  return sum;
}

TMPL Extension
upgrade (const Extension& ext1, const Extension& ext2,
         Element& p1, Element& p2)
{
  if (hard_eq (ext1, ext2)) return ext1;
  Element z1= promote (0, CF(ext1)), z2= promote (0, CF(ext1));
  Extension ext= join (ext1, ext2, z1, z2);
  if (!hard_eq (ext1, ext)) p1= compose (ext, p1, z1);
  if (!hard_eq (ext2, ext)) p2= compose (ext, p2, z2);
  return ext;
}

/******************************************************************************
* Computing the minimal annihilator of an element
******************************************************************************/

TMPL Polynomial
annihilator (const Extension& ext, const Element& p) {
  nat n= deg (ext.mp);
  matrix<C> m (promote (0, CF(ext)), n+1, n);
  Element pp= promote (1, CF(ext));
  for (nat i=0; i<=n; i++) {
    for (nat j=0; j<N(pp); j++) m (i, j)= pp[j];
    pp= rem (p * pp, ext.mp);
  }
  matrix<C> e, k;
  e= row_echelon (m, k);
  for (nat i=0; i<=n; i++) {
    bool ok= true;
    for (nat j=0; j<n; j++)
      if (e (i, j) != 0) ok= false;
    if (ok) return square_free (Polynomial (row (k, i)));
  }
  ERROR ("unexpected situation");
}

TMPL Extension
normalize (const Extension& ext, const Element& p) {
  return Extension (annihilator (ext, p));
}

TMPL Extension
ramify (const Extension& ext, nat p) {
  return Extension (dilate (ext.mp, p));
}

#undef TMPL_DEF
#undef TMPL
#undef Polynomial
#undef Extension
#undef Element
} // namespace mmx
#endif // __MMX_ALGEBRAIC_EXTENSION_HPP