include/algebramix/algebraic_number.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : algebraic_number.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_NUMBER_HPP
#define __MMX_ALGEBRAIC_NUMBER_HPP
#include <algebramix/algebraic.hpp>
#include <numerix/ball_complex.hpp>
namespace mmx {
#define TMPL_DEF template<typename C, typename Ball>
#define TMPL template<typename C, typename Ball>
#define Extension algebraic_extension<C>
#define Field algebraic_number_extension<C,Ball>
#define Polynomial polynomial<C>
#define Element typename algebraic_number_extension<C,Ball>::El
TMPL class algebraic_number_extension;
TMPL bool shrink (const Polynomial& p, Ball& x);
TMPL Polynomial annihilator (const Field& ext, const Element& p);

/******************************************************************************
* The algebraic number types
******************************************************************************/

TMPL_DEF
class algebraic_number_extension {
MMX_ALLOCATORS;
public:
  Extension ext;                      // algebraic extension
  Ball x;                             // ball containing unique root x
  typedef typename Extension::El El;  // representation of elements

public:
  inline algebraic_number_extension ():
    ext (), x (0) {}
  inline algebraic_number_extension (const format<C>& fm):
    ext (fm), x (0) {}
  inline algebraic_number_extension (const Extension& ext2, const Ball& x2):
    ext (ext2), x (x2) { shrink_check (ext.mp, x); }
  inline algebraic_number_extension (const Polynomial& mp, const Ball& x2):
    ext (mp), x (x2) { shrink_check (mp, x); }
  template<typename C2, typename Ball2> inline
  algebraic_number_extension (const algebraic_number_extension<C2,Ball2>& ext2):
    ext (ext2.ext), x (as<Ball> (ext2.x)) { shrink_check (ext.mp, x); }
  inline algebraic_number_extension (const pair<Extension,Ball>& p):
    ext (p.x1), x (p.x2) {}
  inline pair<Extension,Ball> operator * () const {
    return pair<Extension,Ball> (ext, x); }
};

typedef algebraic_number_extension<rational,ball<floating<> > >
        algebraic_real_extension;
typedef algebraic_number_extension<rational,ball<complex<floating<> > > >
        algebraic_complex_extension;
typedef algebraic<rational,algebraic_real_extension>
        algebraic_real;
typedef algebraic<rational,algebraic_complex_extension>
        algebraic_number;

UNARY_RETURN_TYPE(STMPL,Re_op,algebraic_number,algebraic_real);
UNARY_RETURN_TYPE(STMPL,abs_op,algebraic_number,algebraic_real);
BINARY_RETURN_TYPE(STMPL,gaussian_op,algebraic_real,algebraic_real,algebraic_number);

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

TMPL inline nat hash (const Field& x) { return hash (*x); }
TMPL inline nat exact_hash (const Field& x) { return exact_hash (*x); }
TMPL inline nat hard_hash (const Field& x) { return hard_hash (*x); }
TMPL inline bool operator == (const Field& x, const Field& y) {
  return (*x) == (*y); }
TMPL inline bool operator != (const Field& x, const Field& y) {
  return (*x) != (*y); }
TMPL inline bool exact_eq (const Field& x, const Field& y) {
  return exact_eq (*x, *y); }
TMPL inline bool exact_neq (const Field& x, const Field& y) {
  return exact_neq (*x, *y); }
TMPL inline bool hard_eq (const Field& x, const Field& y) {
  return hard_eq (*x, *y); }
TMPL inline bool hard_neq (const Field& x, const Field& y) {
  return hard_neq (*x, *y); }

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

TMPL
struct binary_helper<Field >: public void_binary_helper<Field > {
  typedef pair<Extension,Ball> 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 Field& x) {
    return binary_disassemble<R> (*x); }
  static inline Field assemble (const generic& x) {
    return Field (binary_assemble<R> (x)); }
  static inline void write (const port& out, const Field& x) {
    return H::write (out, *x); }
  static inline Field read (const port& in) {
    return H::read (in); }
};

/******************************************************************************
* Numeric approximation of algebraic numbers
******************************************************************************/

template<typename C, typename V, typename K> K
eval (const polynomial<C,V>& p, const K& x) {
  K sum= 0;
  for (int i=deg(p); i>=0; i--) {
    sum *= x;
    sum += as<K> (p[i]);
  }
  return sum;
}

TMPL bool
improve_zero (const Polynomial& p, Ball& z) {
  typedef Center_type(Ball) CBall;
  typedef Radius_type(Ball) RBall;
  if (p[0] == 0 && z == 0) { z= Ball (0); return true; }
  CBall x= center (copy (z));
  RBall r= radius (z);
  Polynomial dp= derive (p);
  while (true) {
    CBall nx= x - eval (p, x) / eval (dp, x);
    RBall nr= as<RBall> (abs (nx - x));
    while (true) {
      Ball nb= Ball (nx, nr);
      if (eval (dp, nb) == 0) return false;
      RBall rr= abs_up (eval (p, sharpen (nb)) / eval (dp, nb));
      if (rr <= nr) break;
      nr= max (2 * nr, rr);
    }
    bool ok= (2 * nr >= r);
    x= nx;
    r= nr;
    if (ok) break;
  }
  if (!included (Ball (x, r), z)) return false;
  z= Ball (x, r);
  return true;
}

TMPL bool
shrink (const Polynomial& p, Ball& x) {
  nat old_precision= mmx_bit_precision;
  mmx_bit_precision *= 2;
  Ball x2= copy (x);
  bool r= improve_zero (p, x2);
  mmx_bit_precision= old_precision;
  x= copy (x2);
  add_additive_error (x);
  return r;
}

TMPL void
shrink_check (const Polynomial& p, Ball& x) {
  if (!shrink (p, x)) {
    mmerr << "mp= " << p << "\n";
    mmerr << "x = " << x << "\n";
    ERROR ("root not uniquely specified");
  }
}

TMPL void
increase_precision (const Field& ext) {
  typedef Center_type(Ball) CBall;
  typedef Radius_type(Ball) RBall;
  if (precision (ext.x) >= mmx_bit_precision) return;
  Ball new_x= ext.x;
  if (!improve_zero (ext.ext.mp, new_x)) {
    mmerr << "mp= " << ext.ext.mp << "\n";
    mmerr << "x = " << new_x << "\n";
    ERROR ("unexpected situation");
  }
  (const_cast<Field*> (&ext)) -> x= new_x;
}

TMPL Ball
eval (const Field& ext, const Element& p1) {
  increase_precision (ext);
  Ball sum= 0;
  for (int i=deg(p1); i>=0; i--) {
    sum *= ext.x;
    sum += as<Ball> (p1[i]);
  }
  return sum;
}

TMPL Ball
eval (const Field& ext1, const Field& ext2, const vector<C>& v) {
  increase_precision (ext1);
  increase_precision (ext2);
  Ball sum1= 0;
  for (int i1=deg(ext1.ext.mp)-1; i1>=0; i1--) {
    Ball sum2= 0;
    for (int i2=deg(ext2.ext.mp)-1; i2>=0; i2--) {
      sum2 *= ext2.x;
      sum2 += as<Ball> (v[i1*deg(ext2.ext.mp) + i2]);
    }
    sum1 *= ext1.x;
    sum1 += sum2;
  }
  return sum1;
}

TMPL Ball
as_ball (const algebraic<C,Field >& a) {
  return eval (field (a), value (a));
}

TMPL Center_type(Ball )
as_floating (const algebraic<C,Field >& a) {
  return center (as_ball (a));
}

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

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

TMPL inline Element
mul (const Field& ext, const Element& p1, const Element& p2) {
  return mul (ext.ext, p1, p2);
}

TMPL inline Element
invert (const Field& ext, const Element& p1) {
  return invert (ext.ext, p1);
}

TMPL inline Element
div (const Field& ext, const C& c1, const Element& p2) {
  return div (ext.ext, c1, p2);
}

TMPL inline Element
div (const Field& ext, const Element& p1, const Element& p2) {
  return div (ext.ext, p1, p2);
}

TMPL inline bool
is_zero (const Field& ext, const Element& p1) {
  if (deg (p1) <= 0) return is_zero (ext.ext, p1);
  Ball y= eval (ext, p1);
  if (is_non_zero (y)) return false;
  Polynomial ann= annihilator (ext, p1);
  if (ann[0] == 0 && is_non_zero (eval (derive (ann), y))) return true;
  nat old_precision= mmx_bit_precision;
  mmx_bit_precision *= 2;
  bool r= is_zero (ext, p1);
  mmx_bit_precision= old_precision;
  return r;
}

TMPL inline int
sign (const Field& ext, const Element& p1) {
  if (deg (p1) <= 0) return sign (ext.ext, p1);
  Ball y= eval (ext, p1);
  if (is_non_zero (y)) return sign (y);
  Polynomial ann= annihilator (ext, p1);
  if (ann[0] == 0 && is_non_zero (eval (derive (ann), y))) return 0;
  nat old_precision= mmx_bit_precision;
  mmx_bit_precision *= 2;
  int r= sign (ext, p1);
  mmx_bit_precision= old_precision;
  return r;
}

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

TMPL Field
join (const Field& ext1, const Field& 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
  if (deg (ext1.ext.mp) == 1) {
    z1= Polynomial (-ext1.ext.mp[0]/ext1.ext.mp[1]);
    z2= Polynomial (C(1), (nat) 1);
    return ext2;
  }
  else if (deg (ext2.ext.mp) == 1) {
    z1= Polynomial (C(1), (nat) 1);
    z2= Polynomial (-ext2.ext.mp[0]/ext2.ext.mp[1]);
    return ext1;
  }
  else if (hard_eq (ext1, ext2)) {
    z1= Polynomial (C(1), (nat) 1);
    z2= Polynomial (C(1), (nat) 1);
    return ext1;
  }
  else {
    nat n= deg (ext1.ext.mp) * deg (ext2.ext.mp);
    matrix<C> m= transpose (pow_matrix (ext1.ext, ext2.ext));
    matrix<C> u= invert (range (m, 0, 0, n, n));
    vector<C> v= - (u * column (m, n));
    v << C(1);
    Polynomial mp= Polynomial (v);
    Polynomial sf= square_free (mp);
    Extension ext= Extension (sf);
    v= fill<C> (C(0), n);
    v[deg(ext2.ext.mp)]= C(1);
    z1= Element (u * v);
    if (deg (sf) < deg (mp)) z1= rem (z1, sf);
    v= fill<C> (C(0), n);
    v[1]= C(1);
    z2= Element (u * v);
    if (deg (sf) < deg (mp)) z2= rem (z2, sf);
    Ball x;
    nat old_precision= mmx_bit_precision;
    while (true) {
      x= eval (ext1, ext2, column (m, 1));
      if (shrink (ext.mp, x)) break;
      mmx_bit_precision= mmx_bit_precision << 1;
    }
    mmx_bit_precision= old_precision;
    return Field (ext, x);
  }
}

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

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

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

TMPL Polynomial
annihilator (const Field& ext, const Element& p) {
  return annihilator (ext.ext, p);
}

TMPL Field
normalize (const Field& ext, const Element& p) {
  Polynomial mp= annihilator (ext.ext, p);
  Ball z;
  nat old_precision= mmx_bit_precision;
  while (true) {
    z= eval (ext, p);
    if (shrink (mp, z)) break;
    mmx_bit_precision= mmx_bit_precision << 1;
  }
  mmx_bit_precision= old_precision;
  return Field (Extension (mp), z);
}

/******************************************************************************
* Ramification
******************************************************************************/

TMPL Field
ramify (const Field& ext, nat p) {
  Ball z= (p == 2? sqrt (ext.x): exp (log (ext.x) / ((int) p)));
  return Field (ramify (ext.ext, p), z);
}

/******************************************************************************
* Complex arithmetic
******************************************************************************/

inline void
set_imaginary (algebraic_number& z) {
  typedef ball<complex<floating<> > > B;
  polynomial<rational> mp (vec<rational> (1, 0, 1));
  polynomial<rational> i (vec<rational> (0, 1));
  algebraic_complex_extension ext (mp, imaginary_cst<B> ());
  z= algebraic_number (ext, i);
}

inline algebraic_number
times_i (const algebraic_number& z) {
  return z * imaginary_cst<algebraic_number> ();
}

inline algebraic_number
over_i (const algebraic_number& z) {
  return -z * imaginary_cst<algebraic_number> ();
}

inline algebraic_number
gaussian (const algebraic_real& x, const algebraic_real& y) {
  return algebraic_number (x) + times_i (algebraic_number (y));
}

TMPL inline Field
conj (const Field& ext) {
  return Field (ext.ext, conj (ext.x));
}

inline algebraic_number
conj (const algebraic_number& z) {
  return algebraic_number (conj (field (z)), value (z));
}

inline algebraic_real
Re (const algebraic_number& z) {
  algebraic_number x= normalize ((z + conj (z)) / rational (2));
  algebraic_complex_extension cext= field (x);
  algebraic_real_extension    rext (cext.ext, Re (cext.x));
  return algebraic_real (rext, value (x));
}

inline algebraic_real
Im (const algebraic_number& z) {
  return Re (over_i (z));
}

inline algebraic_real
abs (const algebraic_number& z) {
  algebraic_number x= normalize (sqrt (z * conj (z)));
  algebraic_complex_extension cext= field (x);
  algebraic_real_extension    rext (cext.ext, Re (cext.x));
  return algebraic_real (rext, value (x));
}

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

---