synaps/upol/square_free.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS library.
 *   Author(s): F.Carreras, B.Mourrain, O.Ruatta, C.Tanasescu
 *   $Id: SquareFree.h,v 1.1 2005/07/11 08:12:26 mourrain Exp $
 ********************************************************************/
#ifndef SYNAPS_UPOL_SQUAREFREE_H_
#define SYNAPS_UPOL_SQUAREFREE_H_
//--------------------------------------------------------------------
#include <list>
#include <vector>
#include <synaps/init.h>
#include <synaps/arithm/pow.h>
#include <synaps/arithm/sign.h>
#include <synaps/arithm/gmp.h>
#include <synaps/base/Seq.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/gcd.h>


__BEGIN_NAMESPACE_SYNAPS

/*************************************************************************/
template<class C>
UPolDse<C> square_free(const UPolDse<C> & p)
{
  if(p==0)
    return UPolDse<C>(0);
  else
    {
      typedef UPolDse<C> Pol;
      Pol h = p, dp = diff(h);
      Pol g = gcd(h,dp);
      h/=g;
      return h;
    }
}
/*************************************************************************/
template<class C>
UPolDse<UPolDse<C> > square_free(const UPolDse<UPolDse<C> > & p)
{
  typedef UPolDse<UPolDse<C> > PolPol;
  typedef UPolDse<C> Pol;
  PolPol h = p;
  PolPol dp = diff(p);
  PolPol g = gcd(h,dp);
  if (degree(g) >= 0){
  C tmp = content(g[0]);
  for(int i=1;i <= degree(g);++i){
      tmp = gcd(tmp,content(g[i]));
  }
  for(int i=0;i<=degree(g);i++)
    g[i]/=tmp;
  }
  h/=g;
  return h;
}

/*************************************************************************/
template < typename RT >
Seq< UPolDse<RT> >
square_free_factorisation( const UPolDse<RT>& f)
{

  typedef UPolDse<RT>    Poly;

  Seq< Poly > res;

  res.resize( UPOLDAR::degree(f) );

  Poly df = UPOLDAR::diff(f);

  Poly c = gcd(f, df);
  
  if ( UPOLDAR::degree( c) <= 0) {
      res[0] = f;
      res.resize(1);
      return res;
  }
  
  int i = 1;
  int n = 0;
  Poly w = f / c;
  Poly y = df / c;
  Poly z = y - UPOLDAR::diff(w);
  Poly g;
  
  while ( true ) {
      g = gcd( w, z);
      if ( UPOLDAR::degree( g) > 0) {
          res[i] = g;
      }
      ++i;
      w /= g;
      y = z / g;
      z = y - UPOLDAR::diff( w);
      if ( UPOLDAR::degree(z) <= 0) break;
  }
  res[i] = w;
  res.resize(i+1);
  return res;
}


template < typename RT >
Seq< UPolDse<RT> >
square_free_factorization( const UPolDse<RT>& f, SturmSeq<RT>& sth)
{

  typedef UPolDse<RT>    Poly;
  Seq< Poly > res;

  res.resize( UPOLDAR::degree(f) );

  Poly df = UPOLDAR::diff(f);

  Poly c;
  sth = SturmSeq<RT>( f / content(f), df, SUBRES());

  typename SturmSeq<RT>::reverse_iterator it = sth.rep().rbegin();
  
  for (int i = 0; i < sth.size(); ++i) {
    if ( degree( sth[i]) >= 0 )  sth[i] /= content(sth[i]);
  }

  for(; (it!= sth.rep().rend()) && (*it == 0) ; ++it);
  RT temp = content( *it);
  c = (*it) / temp;
  
  if ( UPOLDAR::degree(c) <= 0) {
      res[0] = f;
      res.resize(1);
      return res;
  }
  
  int i = 1;
  int n = 0;

  Poly w = f / c;
  Poly y = df / c;
  Poly z = y - UPOLDAR::diff(w);
  Poly g;
  
  while ( true ) {
      g = gcd( w, z);
      if ( UPOLDAR::degree( g) > 0) {
          res[i] = g;
      }
      ++i;
      w /= g;
      y = z / g;
      z = y - UPOLDAR::diff( w);
      if ( UPOLDAR::degree(z) <= 0) break;
  }
  res[i] = w;
  res.resize(i+1);
  return res;

}

__END_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
#endif // SYNAPS_UPOL_SQUAREFREE_H_