include/continewz/homotopy_floating.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : homotopy_floating.hpp
* DESCRIPTION: Certified homotopies
* COPYRIGHT  : (C) 2009  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_HOMOTOPY_FLOATING_HPP
#define __MMX_HOMOTOPY_FLOATING_HPP
#include <continewz/homotopy_euler.hpp>

namespace mmx {
#define Function slp_tangent<C>

template<typename V= homotopy_euler, typename W= homotopy_euler>
struct homotopy_floating {};

template<typename V> struct homotopy_variant_helper<floating<V> > {
  typedef homotopy_floating<> HV; };
template<typename V> struct homotopy_variant_helper<complex<floating<V> > > {
  typedef homotopy_floating<> HV; };

/******************************************************************************
* Certified homotopies
******************************************************************************/

template<typename C, typename V, typename W>
struct homotopy_stepper_rep<C,homotopy_floating<V,W> >:
  public homotopy_stepper_rep<C,homotopy_abstract>
{
public:
  typedef Fast_type(C) FC;
  typedef homotopy_floating<V,W> S;
  typedef homotopy_stepper<FC,V> Fast_stepper;
  typedef homotopy_stepper<C ,W> Slow_stepper;
  typedef homotopy_stepper<C ,S> Similar_stepper;

public:
  nat          prec;
  nat          half_prec;
  void*        sub_stepper;
  Slow_stepper end_stepper;

public:
  homotopy_stepper_rep (const Function& f):
    homotopy_stepper_rep<C,homotopy_abstract> (f),
    prec (mmx_bit_precision),
    half_prec (((prec + 64) / 128) * 64),
    end_stepper (f)
  {
    if (prec <= 64) sub_stepper= (void*) new Fast_stepper (fast (f));
    else {
      nat old_bit_precision= mmx_bit_precision;
      mmx_bit_precision= half_prec;
      sub_stepper= (void*) new Similar_stepper (copy (f));
      mmx_bit_precision= old_bit_precision;
    }
  }

  ~homotopy_stepper_rep () {
    if (prec <= 64) delete ((Fast_stepper*) sub_stepper);
    else {
      nat old_bit_precision= mmx_bit_precision;
      mmx_bit_precision= half_prec;
      delete ((Similar_stepper*) sub_stepper);
      mmx_bit_precision= old_bit_precision;
    }
  }

  vector<C>
  hom (const vector<C>& init, nat k, const C& t1) {
    return homp (init, k, vec<C> (init[k], t1));
  }

  vector<C>
  homp (const vector<C>& init, nat k, const vector<C>& p) {
    //mmerr << "From x = " << init << " via x[" << k << "] = " << p
    //<< " [floating stepper]" << "\n";
    nat old_bit_precision= mmx_bit_precision;
    mmx_bit_precision= prec;
    vector<C> v;
    if (prec <= 64) {
      //mmerr << "Double precision\n";
      Fast_stepper stepper= *((Fast_stepper*) sub_stepper);
      v= slow<vector<C> > (homx (stepper, fast (init), k, fast (p)));
    }
    else {
      //mmerr << "New precision= " << half_prec << "\n";
      Similar_stepper stepper= *((Similar_stepper*) sub_stepper);
      v= unary_map<copy_op> (homx (stepper, init, k, p));
    }
    v= homx (end_stepper, v, k, p);
    mmx_bit_precision= old_bit_precision;
    return v;
  }
};

#undef Function
} // namespace mmx
#endif // __MMX_HOMOTOPY_FLOATING_HPP

---