synaps/upol/gcd.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS library.
 *   Author(s): F.Carreras, B.Mourrain, O.Ruatta, C.Tanasescu
 *   $Id: Gcd.h,v 1.1 2005/07/11 08:12:26 mourrain Exp $
 ********************************************************************/
#ifndef SYNAPS_UPOL_GCD_H
#define SYNAPS_UPOL_GCD_H
//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/upol.h>
#include <synaps/upol/SturmSeq.h>
#include <synaps/base/Debug_util.h>
//--------------------------------------------------------------------

__BEGIN_NAMESPACE_SYNAPS

template<class C> void contentcoeff( C & a, const C & c ) { a = c; };
template<class C, class POL>
void contentcoeff( C & tmp, const POL & p)
{
  //  typedef typename POL::value_type coeff_type;
  //  typedef typename POL::size_type size_type;
  tmp = C(0);
  int d = degree(p);
  for(int i = 0; i <= d; i++)
    {
      C b;
      contentcoeff(b,p[i]);
      if ( b != C(0) )
        {
          tmp = gcd(tmp,b);
        };
    };
}
template<class C, class UPOL> 
  C contentcoeff( const UPOL & p ) { C tmp; contentcoeff(tmp,p); return tmp; };
template<class POL>
typename POL::value_type content(const POL & p)
{
  typedef typename POL::value_type coeff_type;
  typedef typename POL::size_type size_type;
  size_type d = degree(p);
  coeff_type tmp = p[0];
  for(size_type i = 1; i <= d; i++)
    tmp = gcd(tmp,p[i]);
  return  tmp;
}
template<class R>
double content(const UPolDse<double, R > & p)
{
  double r=1.;
  if (degree(p)>=0)
        r = p[degree(p)];
  return r;
}


template<class C, class R>
UPolDse<C, R> gcd(const UPolDse<C, R> & P, const UPolDse<C, R> & Q)
{
  typedef SturmSeq<C,UPolDse<C,R>,std::vector<UPolDse<C,R> > > stha_t;
  if ((degree(P) < 0))
        return Q;
  else if (degree(Q) < 0)
        return P;
  else if ((degree(P) == 0) || (degree(Q) ==0))
        return UPolDse<C,R>(1);
  else {
    stha_t s(P,Q,HABICHT());
    typename stha_t::reverse_iterator it = s.rep().rbegin();
    for(; (it!= s.rep().rend()) && (*it == 0) ; ++it);
    UPolDse<C,R> RES=*it;
    C tmp=content(RES);
    RES/=tmp;
    return RES;
  }
}


template<class C, class R>
UPolDse<C,R> lcm(const UPolDse<C,R> & p, const UPolDse<C,R> & q)
{
  return (p/gcd(p,q))*q;
}


template<class C> inline
bool gcd_coeffs_stop(const C& p, const type::false_t & ) { return p == C(1); };
template<class C> inline
bool gcd_coeffs_stop(const C& p, const type::true_t &  ) { return degree(p) == 0; };
template<class C> inline
bool gcd_coeffs_stop(const C& p ) { return gcd_coeffs_stop(p, typename ispolynomial<C>::T() ); };

template< typename C, class R >
C gcd_coeffs( const UPolDse<C,R> & pol )
{
  typedef UPolDse<C,R> polynom;
  typedef typename polynom::const_iterator iterator;
  typedef typename polynom::value_type coeff;
  coeff p;
  iterator b = pol.begin();
  ++b;
  p = gcd( *(pol.begin()), *b );
  ++b;
  while( b != pol.end() && ! gcd_coeffs_stop(p)  ) { 
    p = gcd(p,*b);
    ++b; 
  };
  
  if ( b != pol.end() ) return C(1); 
  
  return p;
};
  
  


__END_NAMESPACE_SYNAPS

#endif // SYNAPS_UPOL_GCD_H