synaps/usolve/Aberth.h

Go to the documentation of this file.

/**********************************************************************
 *    This file is part of the source code of the SYNAPS library. 
 *    Author(s): B. Mourrain, GALAAD, INRIA                    
 *    $Id: Aberth.H,v 1.5 2005/09/28 06:30:38 mourrain Exp $
 *********************************************************************/
#ifndef SYNAPS_SOLVE_ABERTH_H
#define SYNAPS_SOLVE_ABERTH_H
#ifdef SYNAPS_WITH_GMP
//----------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/base/Seq.h>
#include <synaps/arithm/gmp.h>
#include <synaps/arithm/CC.h>
#include <mpsolve/mps_poly.h>
//====================================================================
__BEGIN_NAMESPACE_SYNAPS
namespace arithm {
template <> struct realof_<CC> {typedef RR T;};
}
//----------------------------------------------------------------------
namespace aberth_solver
{
#ifdef SYNAPS_WITH_GMP
  inline void assign(mpf_t r, double x)           { mpf_set_d (r,x);}
  inline void assign(mpf_t r, signed long int x)  { mpf_set_si(r,x);}
  inline void assign(mpf_t r, unsigned long int x){ mpf_set_ui(r,x);}
  inline void assign(mpf_t r, mpz_t x)            { mpf_set_z(r,x);}
  inline void assign(mpz_t r, double x)           { mpz_set_d(r,x);}
  inline void assign(mpf_t r, const mpf_t & x)    { mpf_set(r,x);}
  inline void assign(mpc_t r, mpc_t x)            { mpc_set(r,x);}
  inline void assign(mpz_t r, const ZZ& x)  { mpz_set(r, x.get_mpz_t());}
  inline void assign(mpf_t r, const ZZ& x)  { mpf_set_z(r,x.get_mpz_t());}
  inline void assign(mpf_t r, const RR& x)  { mpf_set(r,x.get_mpf_t());}
  inline void assign(mpf_t r, const QQ& x)
  { 
    //    std::cout<<x<<" ";
    RR u=RR(x.get_num())/RR(x.get_den()); mpf_set(r,u.get_mpf_t());
    //    std::cout<<u<<std::endl;
}

  inline void assign(double& d,const mpf_t& r)
  {
    d=mpf_get_d(r);
  }

  template<class T>
  bool assign(std::complex<T> & z, const mpc_t & x, const T& eps)
  {
    T r,i;
    assign(r,x->r);
    assign(i,x->i);
    z = std::complex<T>(r,i);
    return true;
  }

  template<class T>
  inline bool assign(CC& z, const mpc_t & x, const T & eps)
  {
    mpc_set(&z.rep(),x);
    return true;
  }
  inline bool assign(RR& r, const mpc_t & x, const RR & eps)
  {
    mpf_set(r.get_mpf_t(),x->r);
    RR i; mpf_set(i.get_mpf_t(),x->i);
    //DBG_PRINT(x->r<<" "<<x->i);
    if(i<eps && i>-eps) 
      return true;
    else 
      return false;
  }

  inline bool assign(double& d, const mpc_t & x, const RR & eps)
  {
    RR r;
    mpf_set(r.get_mpf_t(),x->r);
    RR i; mpf_set(i.get_mpf_t(),x->i);
    //DBG_PRINT(x->r<<" "<<x->i);
    d=to_double(r);
    if(i<eps && i>-eps) 
      return true;
    else 
      return false;
  }
  // inline void assign(mpz_t r, const Scl<MPZ>& x)  { mpz_set(r,&x.rep());}
  // inline void assign(mpf_t r, const Scl<MPZ>& x)  { mpf_set_z(r,&x.rep());}
  // inline void assign(mpf_t r, const Scl<MPF>& x)  { mpf_set(r,&x.rep());}
#endif
}
//--------------------------------------------------------------------
template<class S, class P>
void run_mpsolve(Seq<S>& r, const P& t, int pi, int po, char* G= "iannc")
{
   mpspoly_t p;

  /* set flags */
  DOLOG = false;
  DOWARN = true;
  DOSORT = true;

  /* set streams */
  outstr = stdout;
  logstr = stderr;

  /* set default values */
  prec_in = pi;         /* if != -1 ignore precision from file */
  prec_out = po;        /* 2*DBL_DIG; default output precision */
  goal = new char[5];
  strncpy(goal, G, 5);  /* default goal */
  //DBG_PRINT(goal<<"\n");
  /* Read polynomial */
  //  read_poly(instr, p);
  p->data_type[0]='d';
  p->data_type[1]='r';
  if(pi==0)
    p->data_type[2]='i';
  else
    p->data_type[2]='f';
  p->prec_in = prec_in;
  mp_set_prec(p->prec_in);
  p->deg = degree(t);
  p->n = p->deg;
  int num_coeff = degree(t) + 1;
  allocate_poly(p);

  if(p->data_type[2]=='i')
    for (int i = 0; i < num_coeff; i++)
      aberth_solver::assign(p->mip_r[i],t[i]);
  else
    for (int i = 0; i < num_coeff; i++)
      aberth_solver::assign(p->mfpr[i],t[i]);
  //mpz_out_str(stdout,10,p->mip_r[i]);

  validate_poly(p, num_coeff);

  /* Set polynomial */
  set_poly(p);

  /* allocate global variables */
  allocate_data();
  /*#G  atexit(free_data);      call free_data at exit */
  mpsolve();

  /* copy roots */
  copy_roots();
  //DBG_PRINT(n);

  S s;
  bool valid;

  typename realof<S>::T eps=1;
  for(int i=0;i<po/2;++i) eps/=10;

  if(n != p->deg)
    r.push_back(S(0));
  for (int i = 0; i < n; i++)
    {
      valid = aberth_solver::assign(s,mroot[i],eps);
      if(valid) 
        {
          r.push_back(s);
        }
    }
}
//----------------------------------------------------------------------
template<class C> inline int input_prec(const C & x) { return -1; }
template<class C> inline int output_prec(const C & x) { return -1; }
template<> inline int input_prec(const double & x) { return 16; }
template<> inline int output_prec(const double & x) { return 16; }

//====================================================================

template<class C>
struct SlvAberth {
  typedef C sol_t;
  int pin, pout;
  SlvAberth():pin(0), pout(DBL_DIG){}
  SlvAberth(int o):pout(o), pin(0){}
  SlvAberth(int o, int i):pin(i), pout(o){}
};

//====================================================================
template<class C>
struct IslAberth {
  typedef C sol_t;
  int pin, pout;
  IslAberth():pin(0), pout(DBL_DIG){}
  IslAberth(int o):pout(o), pin(0){}
  IslAberth(int o, int i):pin(i), pout(o){}
};
//----------------------------------------------------------------------
template<class C, class P> 
void runmpsolve(Seq<C>& r, const P& t, int pi, int po, const SlvAberth<C> & mth)
{
  run_mpsolve(r,t,pi,po,"aannc");
}

template<class C, class P> 
void runmpsolve(Seq<C>& r, const P& t, int pi, int po, const IslAberth<C> & mth)
{
  run_mpsolve(r,t,pi,po,"iannc");
}


#ifdef SYNAPS_WITH_GMP
template<class P> 
int runmpsolve(Seq<RR>& r, const P& t, int pi, int po, const SlvAberth<RR> & mth)
{
  run_mpsolve(r,t,pi,po,"aRrnc");
}

template<class P> 
int runmpsolve(Seq<RR>& r, const P& t, int pi, int po, const IslAberth<RR> & mth)
{
  run_mpsolve(r,t,pi,po,"iRrnc");
}
#endif
//----------------------------------------------------------------------
template<class S, class P>
Seq<S> solve(const P &p, const SlvAberth<S> & A)
{
  Seq<S> r;
  run_mpsolve(r,
              p,
              std::max(A.pin,  input_prec(p[0])), 
              std::max(A.pout, output_prec(p[0])),
              "iannc");
  return r;
}
//--------------------------------------------------------------------
template<>
Seq<double> solve(const UPolDse<double> &p, const SlvAberth<double> & A)
{
  UPolDse<QQ> q(p.size(),AsSize());
  for(unsigned i=0;i<p.size();i++) 
    let::assign(q[i],p[i]);
  return solve(q,A);
}
//--------------------------------------------------------------------
template<class P,class S>
Seq<S> solve(const P &p, const IslAberth<S> & A)
{
  Seq<S> r;
  runmpsolve(r,
             p,
             std::max(A.pin,  input_prec(p[0])), 
             std::max(A.pout, output_prec(p[0])),
             A);
  return r;
}
//--------------------------------------------------------------------
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif // SYNAPS_WITH_GMP
#endif // SYNAPS_SOLVE_ABERTH_H