synaps/upol/bound.h

Go to the documentation of this file.

/*******************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): Elias P. TSIGARIDAS <et@di.uoa.gr>
 *              J.P. Pavone, GALAAD, INRIA
 ********************************************************************/
#ifndef SYNAPS_UPOL_BOUND_H
#define SYNAPS_UPOL_BOUND_H

#include <synaps/init.h>
#include <synaps/base/Debug_util.h>
#include <synaps/base/type.h>
#include <synaps/arithm/gmp.h>
#include <synaps/arithm/double.h>
#include <vector>

__BEGIN_NAMESPACE_SYNAPS

struct NISP{};

struct NISN{};

struct Cauchy {};


template < typename POLY,
           typename M,
           typename T >
T bound(const POLY& p, const M& mth, type::As<T>)
{
  return T(0);
}








template < typename T >
struct abs_max : public std::unary_function<T, void>
{
  abs_max() : max(T(-1)) {}
  void operator() (const T& x)
  {
    using std::abs;
    T temp = abs(x);
    if (temp > max) max = temp;
  }
  T max;
};


template <typename FT, typename POLY>
FT bound(const POLY& p, Cauchy)
{
  using std::abs;
  DEBUG::print("Cauchy");

  typedef typename POLY::value_type     RT;

  abs_max<RT> M = std::for_each(p.begin(), p.end()-1, abs_max<RT>());
  RT an =  abs(UPOLDAR::lcoeff(p));
  return let::convert(RT(M.max + an), an, type::As<FT>());
}


template < typename FT, typename POLY>
FT
max_bound(const POLY& p, Cauchy)
{
  return bound< FT >(p, Cauchy());
}

template < typename FT, typename POLY>
FT
min_bound(const POLY& p, Cauchy)
{
  POLY f(p);
  UPOLDAR::reverse(f, UPOLDAR::degree(p));
  FT tmp =  bound<FT>(f, Cauchy());
  return let::convert(denominator(tmp), numerator(tmp), type::As<FT>());
}

template <typename FT, typename POLY>
FT bound(const POLY& p, NISP )
{
  typedef typename POLY::value_type      RT;
  //  typedef rationalof<RT>                 NTR;
  //  typedef typename NTR::T                FT;
  //typedef double FT;
  FT max(0);
  FT sum(0);

  unsigned i;
  unsigned n = p.size();
  
  // check the degree 
  i = n; while ( p[--i] == 0 ); n = i+1;
  if ( p[i] > 0 )
    {
      i = 0;
      while ( p[i] > 0 && i < n ) i++;
      if ( i == n ) return 0;
      
      for ( unsigned k = i+1; k < n; k ++ )
        if ( p[k] > 0 ) sum += p[k];

      for ( unsigned k = i; k < n; k ++ )
        {
          if ( p[k] < 0 ) 
            { FT tmp = -p[k]/sum; if ( tmp > max ) max = tmp; }
          else 
            sum -= p[k];
        };
    }
  else
    {
      i = 0;
      while ( p[i] < 0 && i < n ) i++;
      if ( i == n ) return 0;
      for ( unsigned k = i+1; k < n; k ++ )
        if ( p[k] < 0 ) sum += p[k];
      for ( unsigned k = i; k < n; k ++ )
        if ( p[k] > 0 ) { FT tmp = -p[k]/sum; if ( tmp > max ) max = tmp; }
        else sum -= p[k];
    };
  max += 1;
  return max;
}

template <typename FT, typename POLY>
FT bound(const POLY& p, NISN )
{
  typedef typename POLY::value_type      RT;
  
  POLY tmp = p;
  for ( unsigned i = 1; i < p.size(); i += 2 ) 
    tmp[i] = -tmp[i];
  return  -bound<FT>(tmp,NISP());
};

__END_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
#endif // SYNAPS_UPOL_BOUND_H