synaps/upol/UPOLDAR.m

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.
*      Author(s): B. Mourrain, GALAAD, INRIA
*
*      Modification: et@di.uoa.gr
*                    Add lcoeff, tcoeff functions
*                    Add diff function that returns the derivative
**********************************************************************/
#ifndef synaps_upol_upoldar_H
#define synaps_upol_upoldar_H
//--------------------------------------------------------------------
#include <iostream>
#include <synaps/init.h>
#include <synaps/arithm/let.h>
#include <synaps/arithm/sign.h>
#include <synaps/linalg/VECTOR.m>
//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//====================================================================
namespace let {
  template<class A,class B> void assign(A&, const B&);
}
namespace VECTOR {}
//====================================================================

namespace UPOLDAR 
{

  using VECTOR::add;
  using VECTOR::sub;
  using VECTOR::mul_ext;
  using VECTOR::div_ext;

  //-------------------------------------------------------------------
  template <class R>
  int degree(const R & p)
  {
    // std::cout <<"(Using default degree) "<<p<<std::endl;
    int i = p.size()-1;
    while(i>-1 && p[i]==0) i--;
    return i;
  }
  //-------------------------------------------------------------------
  template< class R >
  inline
  typename R::value_type
  lcoeff(const R& p)
  {
    return p;
  }
  //-------------------------------------------------------------------
  template< class R >
  inline
  typename R::value_type
  tcoeff(const R& p)
  {
    return p;
  }
  //-------------------------------------------------------------------
//   template <class output, class R>
//   inline output & print(output & os, const R & p)
//   {
//     bool plus=false;
//     if ( degree(p)<0) os << '(' << 0 << ')';
//     else {
//     for(int i= degree(p); i!= -1; i--)
//       {
//         if(p[i]!=0) {
//           if(plus) os<<'+'; else plus=true;
//           print_coeff(os,p[i]);
//           if(i==0);

//           else  {
//             os<<"*x0"; if(i !=1) os<<'^'<<i;
//           }
//         }
//       }
//       if(!plus) os << '(' << 0 << ')';
//     }
//     return os;
//   }
  //-------------------------------------------------------------------
  template <class C>
  std::ostream & print_as_coeff(std::ostream & os, const C & c)
  {
    if(with_plus_sign(c)) 
      os <<"+";
    os<<c;
    return os;
  }
  
  template <class R>
  inline std::ostream & print(std::ostream & os, const R & p,char* x = "x")
  {
    typedef typename R::value_type coeff_t;
    bool plus=false;
    if ( degree(p)<0) os << '(' << 0 << ')';
    else {
    for(int i= degree(p); i!= -1; i--)
      {
        if(p[i]!= 0) {
          plus=true;
          if(i!=degree(p))
            print_as_coeff(os,p[i]);
          else
            os <<p;
          if(i!=0)
            {
              os<<"*"<<x; 
              if(i !=1) os<<'^'<<i;
            }
        }
      }
      if(!plus) os << '(' << 0 << ')';
    }
    return os;
  }
  //--------------------------------------------------------------------
  template <class R, class C>
  inline void set_monomial(R & x, const C & c, unsigned n)
  {

    x.resize(n+1);
    VECTOR::init(x,(typename R::value_type)0);
    x=c;
    checkdegree(x);

  }
  //--------------------------------------------------------------------
  template <class R, class S>
  inline void add_cst( R & r, const R & a, const S & c )
  {
    
  };
  
  template <class R, class S>
  inline void add_cst(R & r, const S & c)
  {
    r+=c;
  }
  //-------------------------------------------------------------------
  template <class R,class A, class B>
  inline void mul_index(R & r, const A & a, const B & b)
  {
    int da=degree(a),
      db=degree(b),i,j;
    for (i=0;i<=db;i++) {
      typename R::value_type bi = b[i];
      if  (bi==0) continue;
      for (j=0;j<=da;j++) {r[i+j] += a[j]*bi;}
    }
  }
  //-------------------------------------------------------------------
  template <class R> inline
  void mul_index_it(R & r, const R & a, const R & b) 
  {
    typename R::const_iterator ia, ib;
    typename R::iterator       ir, it;
    ir=r.begin();
    for (ib=b.begin();ib !=b.end()&& ir != r.end();++ib) {
      typename R::value_type c = *ib;
      it=ir;
      if  (c==0) {++ir;continue;}
      for (ia=a.begin();ia!=a.end();++ia,++it) {
        (*it)+=(*ia)*c;
      }
      ++ir;
    }
  }
  //-------------------------------------------------------------------
  template <class R> inline
  void mul(R & r, const R & a, const R & b)
  {
    VECTOR::init(r,0);
    mul_index(r,a,b);
  }
  //-------------------------------------------------------------------
template <class R> inline
void mul(R & a, const R & b)
{
  if(degree(a)>=0 && degree(b)>=0)
    {
      R ta(a);
      a.resize(degree(a)+degree(b)+1);
      VECTOR::init(a,0);
      mul_index(a,ta,b);
    }
  else
    VECTOR::init(a,0);
}
//--------------------------------------------------------------------
template <class R> inline
void mul_index(R & a, const R & b) 
{
    R ta(a);
    a.resize(degree(a)+degree(b)+1);
    VECTOR::init(a,0);
    mul(a,ta,b);
}
//--------------------------------------------------------------------
template <class R>  inline
void mul_karatsuba(R& r, const R & a, const R & b)
{
     mul_index(r,a,b);
}
//----------------------------------------------------------------------
template <class C, class R>
inline C eval_horner(const R & p, const C & c)
{
  using namespace let;
  C r,cf;
  if(degree(p)>0) {
        typedef typename R::const_reverse_iterator const_reverse_iterator;
        const_reverse_iterator it=p.rbegin();
        assign(r,*it); it++;
        for(; it !=p.rend(); ++it) {r*=c; assign(cf,*it); r+=cf;}
        return r;
  } else if(degree(p)>-1)
  {
        assign(r,p[0]); return r;
  }
  else
        return C(0);
}
//----------------------------------------------------------------------
template <class C, class R>
inline C eval_horner_idx(const R & p, const C & c)
{
  if(degree(p)>0) {
    int d=degree(p);
    C r=p;
    for(int i=d-1; i>-1; --i) {r*=c; r+= p;}
    return r;
  } else if(degree(p)>-1)
    return p[0];
  else
    return C(0);
}

//--------------------------------------------------------------------
template<typename POL, typename X> inline
int sign_at(const POL & p, const X& x)
{
  X v= eval_horner(p,x);
  return sign(v);
}
//--------------------------------------------------------------------
template <class C, class R>
inline C eval(const R & p, const C & c)
{
  return eval_horner(p,c);
}

template<class R>
void div_rem(R & q, R & a, const R & b)
{
  typedef typename R::value_type CF;
  typedef typename R::iterator iterator;
  CF lb = b;
  q=0;
  while (degree(a)>=degree(b)) {
    CF lt(a[degree(a)]/lb);
    int s = degree(a)-degree(b);
    for ( int i = 0; i < degree(b)+1; i ++ )
      a[i+s] -= lt*b[i];
    if ( degree(q) < s ) q.resize(s+1);
    q = lt;
  }

}

//----------------------------------------------------------------------
template<class R> inline void checkdegree(R & p) { 
  int d = p.size()-1; 
  if ( d < 0 ) return;
  while ( p[d] == 0 ) d--;
  p.resize(d+1);
}
//----------------------------------------------------------------------
template <class R> inline  void diff(R & r, const R & p);
//----------------------------------------------------------------------
template <class R> inline  R diff(const R & p);
//----------------------------------------------------------------------
template<class R> void reciprocal(R & w,const R & p);
//----------------------------------------------------------------------
template<class R> void reverse(R & p, typename R::size_type n);
//----------------------------------------------------------------------
template<class R>
typename R::value_type derive(const R p,typename R::value_type x);
//----------------------------------------------------------------------
template<class R> void reduce(R & p, const typename R::size_type & e);
//----------------------------------------------------------------------
template<class R>
  void shift(R & t, const typename R::value_type & x0);
template<class R>
  void shift(R & t,const R & p,const typename R::value_type & x0);
//----------------------------------------------------------------------
template<class R, class C> 
  void scale(R & t,const R & p,const C & l);
//----------------------------------------------------------------------
template<class R> inline
void diff(R & r, const R & p) // version + simple
{
  int n,i;
  r.resize(n=degree(p));
  for( i = 1; i < n+1; ++i) r[i-1] =p[i]*i;
}
  // ancienne version
  // template <class R> inline
  // void diff(R & r, const R & p)
  // {
  //   if(degree(p)>0)//  assert(degree(p)>0);
  //    {
  //      for(unsigned i=1;i<p.size();++i)
  //        r[i-1] =p[i]*i;
  //      r.resize(degree(p));
  //    } else r.resize(0);
  // }

//---------------------------------------------------------------------
template <class R> inline
R diff(const R & p)
{
  R r(p); 
  diff(r,p);
  return r;
}
// ancienne version
// template <class R> inline 
// R diff(const R & p)
// {
//   R r(p);
//   if(degree(p)>0) {
//      for(unsigned i = 1; i < p.size(); ++i)
//        r[i-1] = p[i] * i;
//      r.resize(degree(p));
//   } else r.resize(0);
//   return r;
// }
//---------------------------------------------------------------------
template <class R>
void reciprocal(R & w,const R & p)
{
  typedef typename R::size_type size_type;
  typedef typename R::value_type C;

  const size_type deg = UPOLDAR::degree(p)+1;
  R w0(deg),v(deg),xp(p);
  size_type K = (size_type) (log(p.degree()+1)/log(2)),j=1;
  C p0 = xp;

  xp /= p0;
  w0 =xp/xp;
  v = w0;
  v *= C(-1.0);
  v+=C(2);
  w=v/xp;
  w0=w;
  while (j <= K) {
    UPOLDAR::reduce(w0,pow(2,j+1));
    w=xp*w0;
    v =w;
    v*=C(-1.0);
    v+=C(2);
    UPOLDAR::reduce(v,pow(2,j+1));
    w=w0*v;
    w0=w;
    j++;
  }
   UPOLDAR::reduce(w,deg);
   w/=p0;
}
//-------------------------------------------------------------------------
template <class T> 
inline void reduce(T & p, const typename T::size_type & e)
{
  // il me semble que tout ca c'est la meme chose que p.resize(e) 
  T temp(e);
  for (typename T::size_type i = 0; i < e; i++) temp[i]=p[i];
  p.resize(e);
  //  p.degree()=e-1;
  for (typename T::size_type i = 0; i < e; i++) p[i]=temp[i];
}
//-------------------------------------------------------------------------
template <class T> // version + simple
inline void reverse(T & p, int n)
{ for ( int i = 0; i < n/2; i++ ) std::swap(p[i],p[n-i]) ; }

// template <class T> // ancienne version
// inline void reverse(T & p,typename T::size_type n)
// {
//   for (typename T::size_type i = 0, j=n-1; i<j; i++, j--)
//     {
//       swap(p[i],p[j]);
//     }
// }

//----------------------------------------------------------------------------
template<class O, class R, class I> inline
void eval( O& p, O& dp, const R& Pol, const I& x )
{
  int n = p.size()-1;
  p  = Pol;
  dp = O(0);
  for ( int j = n-1; j >= 0; dp = dp*x + p, p =p*x + Pol, j-- ); 
};

template <class R> inline // wrap sur eval
typename R::value_type derive( const R& Pol, const typename R::value_type& x )
{ 
  typename R::value_type p,res;
  eval(p,res,Pol,x); 
  return res;
};

// template <class R> // ancienne version
// inline typename R::value_type
// derive(const R p, typename R::value_type x)
// {
//   typename R::size_type degree = p.degree();
//   typename R::value_type  deriv(0);
//   for (typename R::size_type j=1; j<=p.degree();j++)
//     deriv+=p[j]*j*pow(x,j-1);
//   return (typename R::value_type) deriv;
// }

//--------------------------------------------------------------------
template< class R> // version + simple
void shift( R& p, const typename R::value_type& c )
{
  int j,k,s;
  for ( s = p.size(), j = 0; j <= s-2; j++ )
    for( k = s-2; k >= j; p += c*p, k-- );
};
//--------------------------------------------------------------------
template<class R> 
void shift(R & r,const R & p,const typename R::value_type & x0)
{
  r=p; shift(r,x0);
}

// ancienne version
// {
//   R dp(p);
//   r[0]=eval_horner(dp,x0);
//   for(unsigned i=1; i< p.size();i++){
//     diff(dp,dp);
//     VECTOR::div(dp,i);
//     r[i]=eval_horner(dp,x0);
//   }
// }
//--------------------------------------------------------------------

template<class R, class C>
void scale(R & r, const R & p, const C & l)
{
  r.resize(p.size());
  r=p;
  typename R::value_type s=l;
  for(unsigned i=1; i< p.size();i++){
    r[i]=p[i]*s;
    s *=l;
  }
}

template<class R>
void coeff_modulo( R & r, const typename R::value_type & x )
{
  for ( typename R::iterator i = r.begin(); i != r.end(); *i %= x, ++i );
};

} //UPOLDAR
/********************************************************************/
__END_NAMESPACE_SYNAPS
/*********************************************************************/
#endif // SYNAPS_UPOL_UPOLDAR_H