/Users/mourrain/Devel/mmx/shape/include/shape/topological_degree.hpp

Go to the documentation of this file.

#ifndef shape_topological_degree_hpp
#define shape_topological_degree_hpp

#include <realroot/sign_variation.hpp>
#include <realroot/Interval.hpp>
#include <realroot/Seq.hpp>
#include <realroot/polynomial_bernstein.hpp>

#define Polynomial polynomial<C, with<Bernstein> >

namespace mmx {

template<class C>
bool not_const_sign(const Polynomial& f)
// si les coefficients de f sont de signes constants not_const_sign = false, 
// sinon not_const_sign = true
{
  bool b = false;
  int d  = f.size();
  int i,j;
  for (i = 0; (i < d) && (f[i] == 0); i++) ;
  
  if (i < d)
  {
    for (j = i+1; (j < d+1)&&(!b); j++)
    {
      b = (sign(f[i]) * sign(f[j]) < 0);
    }
  }
  return b;
}
//----------------------------------------------------------------------
template<class C>
C eval(const Polynomial& f, const C& x)
// evaluation de f en x
{
  return f[0];
}

//----------------------------------------------------------------------
template<class B, class C>
void separate_roots(Seq<Interval<B> >& L, 
                    const Polynomial& bf, const Polynomial& bg, 
                    double eps)
// isolation de racines reelles de l'equation bf*bg = 0 sur [a,b]
{
  typedef Interval<B> interval_t;
  Seq< interval_t > L0; 
  Seq< Polynomial > Lf, Lg; 
  Polynomial f1,f2,g1,g2;
  interval_t I0;

  std::cout<< bf<< std::endl;
  std::cout<< bg<< std::endl;
  if (not_const_sign(bf) || not_const_sign(bg))
  {
    Lf.push_back(bf);
    Lg.push_back(bg);
    L0 << interval_t(0,1);
  }
  B middle = 0.5;
  while (L0.size())
  {
    int i = L0.size()-1;
    I0 = L0[i];
    L0.resize(i);
    f1 = Lf[i]; f2=f1;
    Lf.resize(i);
    g1 = Lg[i]; g2=g1;
    Lg.resize(i);
    
    if (I0.size() < eps)
    {
      L << I0;
    }
    else
    {
      tensor::split(f1,f2,0,0.5);
      tensor::split(g1,g2,0,0.5);

      if (not_const_sign(f2) || not_const_sign(g2))
      {
        Lf.push_back(f2);
        Lg.push_back(g2);
        L0 << interval_t(middle * I0.lower() + (1-middle) * I0.upper(),
                         I0.upper());
      }
        
      if ((f2[0]==0) || (g2[0]==0))
      {
        Lf.push_back(f2);
        Lg.push_back(g2);
        L0 << interval_t(middle * I0.lower() + (1-middle) * I0.upper(),
                         middle * I0.lower() + (1-middle) * I0.upper());
      }
        
      if (not_const_sign(f1) || not_const_sign(g1))
      {
        Lf.push_back(f1);
        Lg.push_back(g1);
        L0<< interval_t(I0.lower(),
                        middle * I0.lower() + (1-middle) * I0.upper());
      }
    }
  }
}

template<class C>
int topological_degree(const Polynomial& f, const Polynomial& g) 
// calcul le degree topologique deg((f,g),(0,0),box)
{
  typedef Interval<double> interval_t;
  //  double a=0, b=1, c=0, d=1;
  
  std::cout<< "B topo"<< std::endl;
  //  print(std::cout,f); 
  //  std::cout<<"f "<<f;std::cout<< std::endl;
  //  print(std::cout,g); std::cout<< std::endl;
  //  std::cout<<"g "<<g;std::cout<< std::endl;

  Seq< interval_t > sa, sb, sc, sd;

  Polynomial ff;
  tensor::face(ff, f, 0, 0);
  //  std::cout<<"ff "<<ff; std::cout<< std::endl;

  Polynomial fg;   //Polynomial h1(face(g.env,0));
  tensor::face(fg, g, 0, 0);
  //  std::cout<<"fg "<<fg; std::cout<< std::endl;

 //  tensor::face(fb, f, 0, 1);
//   std::cout<<"fb "<<fb; std::cout<< std::endl;
 
//   tensor::face(fc, f, 1, 0);
//   std::cout<<"fc "<<fc; std::cout<< std::endl;

//   tensor::face(fd, f, 1, 1);
//   std::cout<<"fd "<<fd; std::cout<< std::endl;

//   std::cout<< "C topo"<<std::endl;


//   tensor::face(gb, g, 0, 1);
//   std::cout<<gb; std::cout<< std::endl;

//   tensor::face(gc, g, 1, 0);
//   std::cout<<gc; std::cout<< std::endl;

//   tensor::face(gd, g, 1, 1);
//   std::cout<<gd;
  return 0;

  // fa(t) = f(a,t), fb(t) = f(b,t), fc(t) = f(t,c), fd(t) = f(t,d)
  // ga(t) = g(a,t), gb(t) = g(b,t), gc(t) = g(t,c), fd(t) = f(t,d)
  
  // sa la liste d'intervalles qui isolent les racines reelles de fa*ga=0
  // sb la liste d'intervalles qui isolent les racines reelles de fb*gb=0
  // sc la liste d'intervalles qui isolent les racines reelles de fc*gc=0
  // sd la liste d'intervalles qui isolent les racines reelles de fd*gd=0
  
//   double eps = 1e-8; 
//   separate_roots(sa, fa, ga, eps);
//   std::cout << "A "<< sa.size()<< std::endl;
//   separate_roots(sb, fb, gb, eps);
//   std::cout << "B "<< sb.size()<< std::endl;
//   separate_roots(sc, fc, gc, eps);
//   std::cout << "C "<< sc.size()<< std::endl;
//   separate_roots(sd, fd, gd, eps);
//   std::cout << "D "<< sd.size()<< std::endl;
  
//  int ja = sa.size();
//   if (ja < 2)
//   {
//     if (ja == 0) sa.resize(1);
//     sa[0].assign(c, d);
//   }
//   else
//   {
//     sa[0].assign(c, 0.5*(sa[0].upper() + sa[1].lower()));
//     for (int i = 1; i < ja-1; i++)
//     {
//       sa[i].assign(sa[i-1].upper(), 0.5*(sa[i].upper() + sa[i+1].lower()));
//     }
//     sa[ja-1].assign(sa[ja-2].upper(), d);
//   }

//   int jb = sb.size();
//   if (jb < 2)
//   {
//     if (jb == 0) sb.resize(1);
//     sb[0].assign(c, d);
//   }
//   else
//   {
//     sb[0].assign(c, 0.5*(sb[0].upper() + sb[1].lower()));
//     for (int i = 1; i < jb-1; i++)
//     {
//       sb[i].assign(sb[i-1].upper(), 0.5*(sb[i].upper()+sb[i+1].lower()));
//     }
//     sb[jb-1].assign(sb[jb-2].upper(), d);
//   }
  
//   int jc = sc.size();
//   if (jc < 2)
//   {
//     if (jc == 0) sc.resize(1);
//     sc[0].assign(a, b);
//   }
//   else
//   {
//     sc[0].assign(a, 0.5*(sc[0].upper() + sc[1].lower()));
//     for (int i = 1; i < jc-1; i++)
//     {
//       sc[i].assign(sc[i-1].upper(), 0.5*(sc[i].upper() + sc[i+1].lower()));
//     }
//     sc[jc-1].assign(sc[jc-2].upper(), b);
//   }

//   int jd = sd.size();
//   if (jd < 2)
//   {
//     if (jd == 0) sd.resize(1);
//     sd[0].assign(a, b);
//   }
//   else
//   {
//     sd[0].assign(a, 0.5*(sd[0].upper() + sd[1].lower()));
//     for (int i = 1; i < jd-1; i++)
//     {
//       sd[i].assign(sd[i-1].upper(), 0.5*(sd[i].upper()+sd[i+1].lower()));
//     }
//     sd[jd-1].assign(sd[jd-2].upper(), b);
//   }

//   int deg = 0;
  
//   for (int i = 0; i < ja; i++)
//   {
//     C u = sa[i].lower(), v = sa[i].upper();
//     int sfu = sign(eval<C>(fa,(u-c)/(d-c))), 
//       sfv = sign(eval<C>(fa,(v-c)/(d-c))), 
//       sgu = sign(eval<C>(ga,(u-c)/(d-c))), 
//       sgv = sign(eval<C>(ga,(v-c)/(d-c)));
//     if (sfu == sfv)
//       deg += sfu * (sgu - sgv);
//     else
//       deg += sgu * (sfv - sfu);
//   }
  
//   for (int i = 0; i < jb; i++)
//   {
//     C u = sb[i].lower(), v = sb[i].upper();
//     int sfu = sign(eval<C>(fb,(u-c)/(d-c))), 
//       sfv = sign(eval<C>(fb,(v-c)/(d-c))), 
//       sgu = sign(eval<C>(gb,(u-c)/(d-c))), 
//       sgv = sign(eval<C>(gb,(v-c)/(d-c)));
//     if (sfu == sfv)
//       deg += sfu * (sgv - sgu);
//     else
//       deg += sgu * (sfu - sfv);
//   }
  
//   for (int i = 0; i < jc; i++)
//   {
//     C u = sc[i].lower(), v = sc[i].upper();
//     int sfu = sign(eval<C>(fc,(u-a)/(b-a))), 
//       sfv = sign(eval<C>(fc,(v-a)/(b-a))), 
//       sgu = sign(eval<C>(gc,(u-a)/(b-a))), 
//       sgv = sign(eval<C>(gc,(v-a)/(b-a)));
//     if (sfu == sfv)
//       deg += sfu * (sgv - sgu);
//     else
//       deg += sgu * (sfu - sfv);
//   }
  
//   for (int i = 0; i < jd; i++)
//   {
//     C u = sd[i].lower(), v = sd[i].upper();
//     int sfu = sign(eval<C>(fd,(u-a)/(b-a))), 
//       sfv = sign(eval<C>(fd,(v-a)/(b-a))), 
//       sgu = sign(eval<C>(gd,(u-a)/(b-a))), 
//       sgv = sign(eval<C>(gd,(v-a)/(b-a)));
//     if (sfu == sfv)
//       deg += sfu * (sgu - sgv);
//     else
//       deg += sgu * (sfv - sfu);
//   }

//   return (int)(deg/8);
}
//======================================================================
#undef Polynomial
//======================================================================
} // namespace mmx
#endif //shape_topological_degree_hpp