synaps/usolve/Symbolic.h

#ifndef synaps_solve_Symbolic_H
#define synaps_solve_Symbolic_H
//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/arithm/sign.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/SturmSeq.h>
#include <synaps/usolve/Algebraic.h>
#include <synaps/arithm/symbolic/real.h>
#include <synaps/base/Seq.h>
#include <synaps/usolve/bezier/IslBzBdg.h>
#include <synaps/upol/SturmSeq.h>
//#include <synaps/arithm/double.h>
//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
struct Bounding;
//--------------------------------------------------------------------
template<class C, class SLV= Bezier<C,Bounding> , class X=double>
//--------------------------------------------------------------------
struct Symbolic
{
  X eps;
  Symbolic():eps(1.e-6) {}
  Symbolic(const X& e):eps(e) {}
  Symbolic(const Symbolic & mth): eps(mth.eps) {}
};

namespace symbolic {

  template<class C>
  struct sqrt_rep : public real_rep<C>
  {
  protected:
    real<C> alpha;
  public:
    //  friend class sqrtof<C>;
    sqrt_rep(real<C> x)
      : real_rep<C>(sqrt(x.lower()),sqrt(x.upper())), alpha(x) {};
    
    sqrt_rep(const C & x)
      :
      real_rep<C>(sqrt(x)-real<C>::eps/2,sqrt(x)+real<C>::eps/2),
      alpha(x)        {};
    
    C prec() { return this->upper()-this->lower(); }
    unsigned improve()
    {
      if(prec()> real<C>::eps)
        {
          C m = (this->lower()+this->upper())*0.5;
          int s =(alpha>m*m);
          if(s==1)
            { this->lower()=m; return 1;}
          else if(s==-1)
            { this->upper()=m; return 1;}
          else
            return 0;
        }
      else
        return 0;
    }
    
    void print(std::ostream & os)
    {
      os<<"Sqrt["<<this->lower()<<", "<<this->upper()<<"] of "<<alpha;
    }
  };

  template<class C>
  struct sqrt
  {
  protected:
    sqrt_rep<C>* _rep;
  public:
    sqrt(const real<C> & x): _rep(new sqrt_rep<C>(x)) {}
    sqrt(const C & x): _rep(new sqrt_rep<C>(x)) {}
    inline operator real<C> () { return real<C>(_rep); }
  };
  
  template<class C>
  real<C> Sqrt(real<C> a) {return sqrt<C>(a);}

  template<class C>
  real<C> Sqrt(const C a) {return sqrt<C>(a);}
  
  //--------------------------------------------------------------------
  template<class C,class R> struct root_rep;
  //--------------------------------------------------------------------
  template<class C,class R>
  int sturmcompare(root_rep<C,R>* a, root_rep<C,R>* b)
  {
    if (a->upper() < b->lower())
      return -1;
    else if (b->upper() < a->lower())
      return 1;
    else if (a->poly() == b->poly())
      {
        if (a->su != b->sl)
          return 0;
        else
          if(a->lower() < b->lower())
            return -1;
          else if(a->lower()> b->lower())
            return 1;
          else
            return 0;
      }
    else
      {
        // SturmSeq<R> s(a->poly(), Diff(a->poly())*b->poly(), SUBRES());
        SturmSeq<R> s(a->poly(), b->poly(), SUBRES());
        R l= std::max(a->lower(),b->lower()),
          u= std::min(a->upper(),b->upper());
        int v0 = variation(s,l) - variation(s,u);
        if(v0==0)
          if(a->lower() < b->lower())
            return -1;
          else if(a->lower()> b->lower())
            return 1;
          else
            return 0;
        else
          return v0 * b->su;
      }
  }

  //--------------------------------------------------------------------
  template<class R, class ZT=R>
  struct root_rep: public real_rep<R>, ALGEBRAIC::root_of<ZT>
  {
    int    sl, su;
    
    typedef typename ALGEBRAIC::root_of<ZT>::FT FT;

    root_rep(const UPolDse<ZT> & t, const FT & a, const FT &b)
      : ALGEBRAIC::root_of<ZT>(t,a,b)
    {
      let::assign(this->lower(),a);
      let::assign(this->upper(),b);
      sl=sign(UPOLDAR::eval_horner(this->poly(),this->lower()));
      su=sign(UPOLDAR::eval_horner(this->poly(),this->upper()));
    }
    root_rep(const UPolDse<ZT> & t, const FT & a, const FT &b, unsigned idx)
      : ALGEBRAIC::root_of<ZT>(t,a,b,idx)
    {
      let::assign(this->lower(),a);
      let::assign(this->upper(),b);
      sl=sign(UPOLDAR::eval_horner(this->poly(),this->lower()));
      su=sign(UPOLDAR::eval_horner(this->poly(),this->upper()));
    }
    
    root_rep(char *s, const R& a, const R& b)
      : real_rep<R>(a,b),
        ALGEBRAIC::root_of<ZT>(s)
    {
      sl=sign(UPOLDAR::eval_horner(this->poly(),a));
      su=sign(UPOLDAR::eval_horner(this->poly(),b));
    }
    
    UPolDse<ZT> pol() { return this->poly(); }
    
    unsigned improve()
    {
      if (prec()> real<R>::eps)
        {
          int s = sign(UPOLDAR::eval_horner(this->poly(),R((this->lower()+this->upper())/2)));
          if(s==0) {
            this->lower()=(this->lower()+this->upper())/2; this->upper()=this->lower();
          }
          else if(s*sl<0) {
            this->upper()=(this->lower()+this->upper())/2; su=s;
          }
          else {
            this->lower()=(this->lower()+this->upper())/2; sl=s;
          }
          return 1;
        }
      else
        return 0;
    }
    
    R prec() { return this->upper()-this->lower(); }
    
    int compare(real_rep<R>* b)
    {
#if 1
      if (this->upper() < b->lower()) return -1;
      else if (b->upper() < this->lower()) return 1;
      else if ( (this->prec() < real<R>::eps) &&
                (b->prec() < real<R>::eps) )
        return sturmcompare(this, dynamic_cast<root_rep<R,ZT>*>(b));
      else
        {
          int i1=this->improve(), i2=b->improve();
          if(i1 || i2)
            return (this->compare(b));
        }
#endif
      return sturmcompare(this, dynamic_cast<symbolic::root_rep<R,ZT>*>(b));
    }
    
    void print(std::ostream & os) const
    {
      os<<"Rootof(";
      UPOLDAR::print(os,this->poly());
      os<<", ["<<this->lower()<<", "<<this->upper()<<"])";
    }
  };

  //--------------------------------------------------------------------
  template<class R, class Z=R>
  struct root
  {
    typedef typename ALGEBRAIC::root_of<Z>::FT FT;
  protected:
    root_rep<R,Z>* rep_;
  public:
    root(const UPolDse<Z> & t, const FT& a, const FT& b)
      :rep_(new root_rep<R,Z>(t,a,b)) {}
    root(const UPolDse<Z> & t, const FT& a, const FT& b,unsigned idx)
      :rep_(new root_rep<R,Z>(t,a,b,idx)) {}
    root_rep<R,Z>* operator -> () {return rep_;}
    root_rep<R,Z>* operator -> () const {return rep_;}
    root_rep<R,Z>  rep         () {return *rep_;}

    inline operator real<R> () { return real<R> (rep_); }
  };
} //symbolic
//====================================================================
template<class R> 
std::ostream& operator<<(std::ostream& os, const symbolic::root<R> &a)
{
  a->print(os);
}
//--------------------------------------------------------------------
// Comparison operators:
//--------------------------------------------------------------------
template<class R> inline
bool operator< (symbolic::real<R> a, symbolic::real<R> b) 
{ 
  return (a.compare(b)>0); 
}
template<class R> inline
bool operator> (symbolic::real<R> a, symbolic::real<R> b) 
{  
  return (a.compare(b)<0); 
}
template<class R> inline
bool operator< (const R& a, symbolic::real<R> b) { return (compare(b,a)<0); }
template<class R> inline
bool operator> (const R& a, symbolic::real<R> b) { return (compare(b,a)>0); }
template<class R> inline
bool operator< (symbolic::real<R> a, const R& b) { return (compare(a,b)>0); }
template<class R> inline
bool operator> (symbolic::real<R> a, const R& b) { return (compare(a,b)<0); }

//--------------------------------------------------------------------
template < typename C_T,
           typename UPOL_T >
symbolic::root<C_T,typename UPOL_T::coeff_t> RealRoot(const UPOL_T& p, unsigned idx)
{
  typedef typename UPOL_T::coeff_t RT;
  typedef typename Sturm<RT>::RO_t                   RO_t;
  typedef std::vector<RO_t>                          Cont;
  typedef std::back_insert_iterator< Cont >          OutIter;
  
  Cont sol;
  OutIter tt(sol);
  Sturm<RT> mth;
//  ALGEBRAIC::Solve(p, tt, mth);
  ALGEBRAIC::solve_sturm(p, tt, mth);
  std::stable_sort( sol.begin(), sol.end(), ALGEBRAIC::Refine_compare());

  Seq<RO_t> solutions(sol);
  assert(idx<solutions.size());
  RO_t s=solutions[idx];
  return symbolic::root<C_T,RT>(s.poly(),s.left(),s.right(),idx);
}
//--------------------------------------------------------------------
template<typename T, typename C>
inline bool has_type(symbolic::real_rep<C>* object)
{
    return dynamic_cast<T*>(object) != 0;
}



//--------------------------------------------------------------------
template<class FT, class C, class MTH>
void AppendSol(Seq< symbolic::real<FT> > & res,
               const UPolDse<C> & P, 
               const Symbolic<FT,MTH> & mth,
               const FT& a, const FT & b)//, const SturmSeq<C> & s)
{
  assert(a<=b);
  //  std::cout<<"Closed root "<<P<<" "<<a<<" "<<b<<std::endl;
  //  typedef typename NumberTraits<C>::FT FT;
  typedef FT sol_t;
  Seq<Interval<sol_t> > sol = solve(P,IslBzBdg<FT>(1e-2),a,b);
  //Seq<sol_t> sol; AppendSol(sol,P,Bezier<sol_t,SyNum>(1e-2),a,b);
  //  AppendSol(sol,P,Bezier<sol_t,Isolation>(1e-2),a,b);
  int v0=0;
  if(v0==1)
    res.push_back(symbolic::root<FT,C>(P,a,b));
  else if(v0>1)
    {
      std::cout<<"Isolation still needed "<<a<<" "<<b<<std::endl;
//       FT m=(a+b)/2;
//       AppendSol(res,P,mth,a,m);
//       AppendSol(res,P,mth,m,b);
    }
}
//--------------------------------------------------------------------
template<class FT, class C, class MTH>
void AppendSol(Seq< symbolic::real<FT> > & res,
               const UPolDse<C> & P, const Symbolic<FT,MTH> & mth)
{
  using let::assign;
  typedef typename NumberTraits<C>::FT  field_t;
  typedef FT sol_t;
  SturmSeq<C> s(P, diff(P), HABICHT());
  if(degree(s[s.size()-1])!=0)
    {
      std::cout<<"Not square free"<<std::endl;
      //      AppendSol(res, P/s[s.size()-1],mth);
    }
  else 
    {
      Seq<Interval<double> > sol = solve(P, IslBzBdg<double>(1e-3));
      //  Seq<double> sol = Isole(P, Bezier<double,Bounding>(1e-3));
      //Seq<sol_t> sol = Isole(P, MTH(1e-6));
      //  std::cout<<sol<<std::endl;
      
      FT a,b;
      for(unsigned i=0; i<sol.size(); i++)
        {
          assign(a,sol[i].inf()); assign(b,sol[i].sup());
          //      AppendSol(res,P,Bezier<field_t,Isolation>(),a,b);
          //      std::cout<<a<<" "<<b<<" "<<v0<<std::endl;
          if(a!= b)
            {
              int v0 = variation(s,a) - variation(s,b);
              if(v0==1)
                res.push_back(symbolic::root<FT,C>(P, a, b));
              else if(v0>1)
                AppendSol(res,P,mth,a,b);
        }
          else if(P(a)==0)
            res.push_back(symbolic::root<FT,C>(P, a, a));
          
        }
    }
}

//--------------------------------------------------------------------
template<class C, class FT>
Seq< symbolic::real<FT> >
solve(const UPolDse<C> & P, const Symbolic<FT> & mth)
{ 
  Seq< symbolic::real<FT> > res;
  //  UPolDse<C> Pr = square_free(P);
  UPolDse<C> Pr = P;
  AppendSol(res, Pr, mth);
  return res;
}
//--------------------------------------------------------------------
__END_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
#endif // synaps_solve_Symbolic_H