synaps/upol/quotient.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, H. Prieto, GALAAD, INRIA                 *
**********************************************************************
History:
    1998/07/01 first version by HP.
    2000/03/02 updated by BM, use modules.
**********************************************************************/
#ifndef SYNAPS_UPOL_MODUPOL_H
#define SYNAPS_UPOL_MODUPOL_H
//------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/base/type.h>
#include <synaps/upol/UPOLDAR.m>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/horner.h>
#include <synaps/linalg/MatrStr.h>
#include <synaps/linalg/hankel.h>
#ifdef FFT
#include <synaps/upol/mult_fft.h>
#endif
#include <synaps/upol/resultant.h>
//----------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
namespace quotient 
{
  template<class C, class R=UPolDse<C> >
    struct ModUPol
    {
      
      R poly;
      
      typedef typename R::value_type       value_type;
      typedef typename R::size_type        size_type;
      typedef typename R::iterator         iterator;
      typedef typename R::reverse_iterator iterator_r;
      typedef upol::horner<C>              rep_type;
      typedef ModUPol<C,R>                 self_t;
      
      ModUPol():poly() {}
      ModUPol(const R & r): poly(r) {}
      ModUPol(const ModUPol<R> & Mod):  poly(Mod.poly) {}
      
      template<class S>
      void  monomial_to_horner(rep_type& h, const S & pol);
      template<class S>
      void  horner_to_monomial(S& h, const rep_type & h);
      
      // Returns the quotient dimension
      size_type dim() {return UPOLDAR::degree(poly);}
      
      void  assign(rep_type& r, const rep_type & a)   {r=a;}
      void  assign(rep_type& r, const R & a)          {monomial_to_horner(r,a);}
      void  assign(rep_type& r, const value_type & a)
      {r.resize(dim()); r[dim()-1]=a/poly[dim()];}
      void  assign(rep_type& r, int a) {assign(r,value_type(a));}

      void  add(rep_type& r, const rep_type & a, const rep_type& b);
      void  add(rep_type& r, const rep_type & a);
      void  sub(rep_type& r, const rep_type & a, const rep_type& b);
      void  sub(rep_type& r, const rep_type & a);
      void  mul(rep_type& r, const rep_type & a, const rep_type& b);
      void  mul(rep_type& r, const rep_type & a, const C& b);
      void  div(rep_type& r, const rep_type & a, const C& b);
      
      void  square(rep_type& r, const rep_type & a);
      
      void  print(std::ostream & os, const rep_type & a)
      { R p; horner_to_monomial(p,a); UPOLDAR::print(os,p); }
      void  read (std::istream & is, rep_type & a);
    };

    template <class C, class R>
      void ModUPol<C,R>::add(rep_type & r,const rep_type & a, const rep_type& b)
      {
        r.resize(dim());  UPOLDAR::add(r,a,b);
      }
    
      template <class C, class R>
        void ModUPol<C,R>::add(rep_type & r,const rep_type & a)
        {
          r.resize(dim());  UPOLDAR::add(r,a);
        }
      
        template <class C, class R>
          void ModUPol<C,R>::sub(rep_type & r,const rep_type & a, const rep_type& b)
          {
            r.resize(dim());  UPOLDAR::sub(r,a,b);
          }

template <class C, class R>
void ModUPol<C,R>::sub(rep_type & r,const rep_type & a)
{
  r.resize(dim());  UPOLDAR::sub(r,a);
}

template <class C, class R>
void ModUPol<C,R>::square(rep_type & x,const rep_type & f)
{
  assert(&x != &f);

  typedef typename R::value_type  value_type;
  typedef typename R::size_type   size_type;

  size_type degree = dim();
  size_type M = 2*degree;

  // Compute the convolution product poly * f,
  // by a product of the FFT of poly and f
  x.resize(M); VECTOR::init(x,0);

#ifdef FFT
  UPOLDAR::mul_fft(x, poly, f);
#else
  UPOLDAR::mul_index(x, poly, f);
  //  VECTOR::print(std::cout,x);std::cout<<std::endl;
#endif

  // Change the signs of the last degree terms of poly*f
  for (size_type i = M-1; i>=degree; i--) x[i]*= value_type(-1);

  // Now, we compute the convolution of this result and f
  linalg::rep1d<value_type> res(M+degree-1);

#ifdef FFT
  UPOLDAR::mul_fft(res,x,f);
#else
  UPOLDAR::mul_index(res,x,f);
  //  VECTOR::print(std::cout,res);std::cout<<std::endl;
#endif

  // We keep only the middle part.
  for (size_type j = 0, i = degree; i<M; i++,j++)
    x[j] = value_type(-1)*(res[i]);

  for(size_type i = degree; i<x.size();++i) x[i]=0;

  UPOLDAR::checkdegree(x);
  //      delete [] res;
}

//----------------------------------------------------------------------
template <class C,class R>
void ModUPol<C,R>::mul(rep_type & x, const rep_type & a, const rep_type& b)
{
  typedef typename R::size_type   size_type;
  typedef typename R::value_type  value_type;

  size_type degree = dim();
  rep_type res(degree,type::AsSize()), temp(degree,type::AsSize());

  add(temp,a,b);
  square(x,temp);
  sub(temp,a,b);
  square(res,temp);
  sub(x,res);
  VECTOR::div_ext(x,value_type(4));
  //      delete [] res.tab_;
}

//----------------------------------------------------------------------
template <class C, class R>
template<class S>
void ModUPol<C,R>::monomial_to_horner(rep_type &x, const S &  pol)
{
  typedef typename R::value_type  value_type;
  typedef typename R::size_type   size_type;
  typedef typename R::rep_type    rep_type;

  size_type degree = dim();
  rep_type h1(2*degree-1, type::AsSize());
  h1[degree-1]=1/poly[degree];
  for (size_type i = degree; i < 2*degree-1; i++)
    {
      for(size_type j=0; j<degree; j++)
        h1[i] -= h1[i-degree+j]*poly[j];
      h1[i]/=poly[degree];
    }
  //VECTOR::print(std::cout,h1);std::cout<<std::endl;

  /* If x is expressed in the monomial basis and we want it in
     the Horner basis, we compute x = H1 * x. */
  MatrStr<linalg::hankel<value_type> >  P(degree, degree, h1.begin());
  rep_type tmp(pol.begin(),pol.end());
  tmp.resize(degree);
  x.resize(degree); VECTOR::init(x,0);
  using namespace MATRIX; mul_vect(x,P,tmp);
  //  VECTOR::print(std::cout,x);std::cout<<std::endl;
  VECTOR::reverse(x,degree);
}

//----------------------------------------------------------------------
template <class C, class R>
template<class S>
void ModUPol<C,R>::horner_to_monomial(S &x, const rep_type &  pol)
{
  typedef typename R::value_type  value_type;
  typedef typename R::size_type   size_type;
  typedef typename R::rep_type    rep_type;

  size_type degree = dim();
  rep_type h1(2*degree,type::AsSize());
  R r(degree,type::AsSize());

  /* if x is expressed in the Horner basis and we want it
     in the monomial basis, we compute x = P*x. */
  MatrStr<linalg::hankel<value_type> > P =bezout(poly);
  //  cout <<P<<endl;
  S tmp(pol.begin(),pol.end()); tmp.resize(degree);
  VECTOR::reverse(tmp,degree);
  x.resize(degree); VECTOR::init(x,0);
  using namespace MATRIX; mul_vect(x,P,tmp);
 }
//----------------------------------------------------------------------
} //namespace quotient
__END_NAMESPACE_SYNAPS
//----------------------------------------------------------------------
#endif // SYNAPS_UPOL_MODUPOL_H