synaps/usolve/Exclusion.h

Go to the documentation of this file.

/********************************************************************
 *      This file is part of the source code of SYNAPS kernel.
 *      Author(s): B. Mourrain, M. Racois, GALAAD, INRIA    
 *      History: First draft, BM. 28/05/2000.
 *        Cleaning adjacent box; M. Racois 24/07/2000.
 *      $Id: Exclusion.h,v 1.1 2005/07/11 10:03:51 jppavone Exp $
 ********************************************************************/
#ifndef SYNAPS_SOLVE_EXCLUSION_H
#define SYNAPS_SOLVE_EXCLUSION_H
/********************************************************************/
#include <synaps/init.h>
#include <synaps/base/type.h>
#include <synaps/base/Seq.h>
#include <synaps/upol/UPOLDAR.m>

#include <cmath>


__BEGIN_NAMESPACE_SYNAPS

using std::abs;

struct Exclusion 
{
  double eps;

  Exclusion(): eps(1.e-6) {}
  Exclusion(double e): eps(e) {}
  Exclusion(const Exclusion & E): eps(E.eps) {}
  template<class R>
  typename R::value_type operator()(const R & p, 
                                   const typename R::value_type & x0,
                                   const typename R::value_type & r) const;
};

//--------------------------------------------------------------------
template<class R>
typename R::value_type Exclusion::operator()(const R & p, 
                               const typename R::value_type & x0,
                               const typename R::value_type & r)  const
{
  R tmp=p;
  using UPOLDAR::shift; shift(tmp,x0);
  tmp[0]= abs(tmp[0]);
  for(int i=1; i <=tmp.get_degree(); i++) tmp[i]= -abs(tmp[i]);
  return tmp(r);
}

//--------------------------------------------------------------------
template<class C> struct box 
{
  C center, radius;
  box(): center(0), radius(0){}
  box(C c, C r): center(c), radius(r){}
  box(const box<C> & b): center(b.center), radius(b.radius){}
  box<C>& operator=(const box<C> & b) 
    {center =b.center; radius=b.radius; return *this;}
};

//--------------------------------------------------------------------
template<class POL, class C>
  Seq<C> solve(const POL& p,const Exclusion & EX, const C& a,const C& b)
{
        std::vector<box<C> >  lst;
  Seq<C> sol;
  lst.push_back(box<C>((a+b)/2, abs(b-a)/2));
  typename std::vector<box<C> >::iterator ib=lst.begin();
  while(lst.size()!=0)
    {
      //      cout <<"B("<<ib->center<<","<<ib->radius<<")"<<endl;
      if(EX(p,ib->center,ib->radius)<0)
        {
          if(ib->radius<EX.eps) 
          { 
            sol.push_back(ib->center);
            lst.erase(ib);ib=lst.end();ib--;
          } else {
            C nr=(ib->radius/2);
            box<C> nb(ib->center+nr,nr);
            ib->center-=nr;
            ib->radius/=2;
            ib=lst.insert(ib,nb);
            ib++;
          }
        }
      else
        {          
          lst.erase(ib);
          ib=lst.end();ib--;
        }
    }
  //  cout<<sol<<endl;
  
   Seq<C> liste;
   for(typename std::vector<C>::iterator iter=sol.begin();iter!=sol.end();iter++){
     C ab=*iter;
     C ac=*iter;
     while ((abs(*iter-*(iter+1))<EX.eps*3) && ((iter+1)!=sol.end()))  
       {iter++;ac=*iter;}
     liste.push_back((ab+ac)/2);
   } 
   return (liste);
}

template<class POL, class C> inline
Seq<C> solve(const POL& p, const box<C> & b, const Exclusion & mth)
{
  return solve(p,b.center-b.radius,b.center+b.radius,mth);
}



__END_NAMESPACE_SYNAPS


/********************************************************************/
#endif // SYNAPS_SOLVE_EXCLUSION_H