include/continewz/analytic_meta.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : analytic_meta.hpp
* DESCRIPTION: Routines M: Analytic C -> Analytic C, such that M(f) and f
*              are fundamentally equal, but only differ in the way
*              the underlying algorithms for analytic continuation work.
* COPYRIGHT  : (C) 2010  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_ANALYTIC_META_HPP
#define __MMX_ANALYTIC_META_HPP
#include <continewz/analytic.hpp>

namespace mmx {
#define TMPL template<typename C, typename V>
#define Analytic analytic<C,V>
#define Analytic_rep analytic_rep<C,V>
#define Assume assume_bound<typename unvectorize<R >::val>
#define Scalar typename unvectorize<C>::val
#define Radius Abs_type(Scalar)
#define Series series<C>
#define Series_rep series_rep<C>
#define Polynomial polynomial<C>
#define R Abs_type(C)

/******************************************************************************
* Fast evaluation in limited precision
******************************************************************************/

TMPL
class fast_eval_analytic_rep: public Analytic_rep {
protected:
  typedef Fast_type(C) FC;
  typedef typename unvectorize<FC>::val FScalar;

  Analytic        f;
  polynomial<FC>* fp;

public:
  inline fast_eval_analytic_rep (const Analytic& f2):
    Analytic_rep (CF(f2)), f (f2), fp (NULL) {}
  inline ~fast_eval_analytic_rep () {
    if (fp != NULL) mmx_delete_one (fp); }
  void Clear_cache (nat which) const {
    Analytic_rep::Clear_cache (which);
    clear_cache (f, which); }
  syntactic expression (const syntactic& z) const {
    return apply ("fast_eval", flatten (f, z)); }
  Series Expand () const {
    return expand (f); }
  Radius Radius_bound (nat order) const {
    return radius_bound (f, order); }
  R Tail_bound (const Radius& r, nat order, Assume& a) const {
    return tail_bound (f, r, order, a); }
  Analytic Move (const Scalar& z) const {
    return fast_eval (move (f, z)); }
  C Eval (const Scalar& z) const {
    nat prec = precision (as<FScalar> (promote (0, z)));
    nat order= (nat) (mmx_order_ratio * ((double) prec));
    if (fp == NULL) {
      fast_eval_analytic_rep<C,V>* me=
        const_cast<fast_eval_analytic_rep<C,V>*> (this);
      polynomial<C> p= truncate (f, order);
      me->fp= mmx_new_one<polynomial<FC> > (fast (p));
    }
    FScalar fz= fast (z);
    FC sum= (*fp)[order-1];
    for (int i= order-2; i>=0; i--)
      sum= sum * fz + (*fp)[i];
    return slow<C> (sum); }
  Analytic Derive () const {
    return fast_eval (derive (f)); }
};

TMPL Analytic
fast_eval (const Analytic& f) {
  // radial (f) returns the same analytic functions as f,
  // but with special optimizations for fast numeric evaluation
  return (Analytic_rep*) new fast_eval_analytic_rep<C,V> (f);
}

/******************************************************************************
* Cache continuation along last path
******************************************************************************/

TMPL
class analytic_next {
MMX_ALLOCATORS
public:
  Scalar   next_point;
  Radius   next_radius;
  Analytic next_function;
};

TMPL
class cache_last_analytic_rep: public Analytic_rep {
protected:
  Analytic f;
  analytic_next<C,V>* next;

protected:
  void Set_next (const Scalar& z1) const {
    //mmout << "    next -> " << z1 << "\n";
    typedef cache_last_analytic_rep<C,V> Me;
    Me* me= const_cast<Me* > (this);
    if (next == NULL) me->next= mmx_new_one<analytic_next<C,V> > ();
    me->next->next_point   = z1;
    me->next->next_function= cache_last (move (f, z1));
    me->next->next_radius  = radius_bound (next->next_function);
  }

public:
  inline cache_last_analytic_rep (const Analytic& f2):
    Analytic_rep (CF(f2)), f (f2), next (NULL) {}
  inline ~cache_last_analytic_rep () {
    if (next != NULL) mmx_delete_one (next); }
  void Clear_cache (nat which) const {
    Analytic_rep::Clear_cache (which);
    clear_cache (f, which); }
  syntactic expression (const syntactic& z) const {
    return apply ("cache_last", flatten (f, z)); }
  Series Expand () const {
    return expand (f); }
  Radius Radius_bound (nat order) const {
    return radius_bound (f, order); }
  R Tail_bound (const Radius& r, nat order, Assume& a) const {
    return tail_bound (f, r, order, a); }
  Analytic Move (const Scalar& z) const {
    return cache_last (move (f, z)); }
  C Eval (const Scalar& z) const {
    return eval (f, z); }
  Analytic Move (const list<Scalar >& l2) const {
    //mmout << "l2= " << l2 << "\n";
    if (is_nil (l2)) return this->me ();
    Radius r1= shrink (radius_bound (f));
    list<Scalar > l= simplify (l2, r1);
    //mmout << "l = " << l << "\n";
    Scalar z1= car (l);
    if (next == NULL) Set_next (z1);
    if (hard_eq (z1, next->next_point))
      return move (next->next_function, cdr (l));
    Scalar z2= next->next_point;
    Radius r2= shrink (next->next_radius);
    if (is_nil (cdr (l)) || abs (z2 - z1) >= r2) {
      //mmout << "z1= " << z1 << "\n";
      //mmout << "z2= " << z2 << "\n";
      //mmout << "r1= " << r1 << "\n";
      //mmout << "r2= " << r2 << "\n";
      Set_next (z1);
      return move (next->next_function, cdr (l));
    }
    else return move (next->next_function, cons (z1 - z2, cdr (l))); }
  C Eval (const list<Scalar >& l2) const {
    if (is_nil (l2)) return this->me () [0];
    Radius r1= shrink (radius_bound (f));
    list<Scalar > l= simplify (l2, r1);
    Scalar z1= car (l);
    if (is_nil (cdr (l))) return eval (f, z1);
    if (next == NULL) Set_next (z1);
    if (hard_eq (z1, next->next_point))
      return eval (next->next_function, cdr (l));
    Scalar z2= next->next_point;
    Radius r2= shrink (next->next_radius);
    if (is_nil (cdr (l)) || abs (z2 - z1) >= r2) {
      Set_next (z1);
      return eval (next->next_function, cdr (l));
    }
    else return eval (next->next_function, cons (z1 - z2, cdr (l))); }
  Analytic Derive () const {
    return cache_last (derive (f)); }
};

TMPL Analytic
cache_last (const Analytic& f) {
  // cache_last (f) returns the same analytic functions as f,
  // but with special optimizations for accelerating continuation
  // along several slightly different paths
  return (Analytic_rep*) new cache_last_analytic_rep<C,V> (f);
}

TMPL Analytic
radial (const Analytic& f) {
  return cache_last (fast_eval (f));
}

/******************************************************************************
* Rough analytic continuation based on heuristic algorithms
* for the computation of the convergence radius and tail bounds
******************************************************************************/

TMPL
class rough_analytic_rep: public Analytic_rep {
  Analytic f;

public:
  inline rough_analytic_rep (const Analytic& f2):
    Analytic_rep (CF(f2)), f (f2) {}
  void Clear_cache (nat which) const {
    Analytic_rep::Clear_cache (which);
    clear_cache (f, which); }
  syntactic expression (const syntactic& z) const {
    return apply ("rough", flatten (f, z)); }
  Series Expand () const {
    return expand (f); }
  Radius Radius_bound (nat order) const {
    return heuristic_radius (truncate (f, order), order); }
  R Tail_bound (const Radius& r, nat order, Assume& a) const {
    return tail_bound (f, r, order, a); }
  Analytic Move (const Scalar& z) const {
    return rough (move (f, z)); }
  Analytic Derive () const {
    return rough (derive (f)); }
};

TMPL Analytic
rough (const Analytic& f) {
  return (Analytic_rep*) new rough_analytic_rep<C,V> (f);
}

#undef TMPL
#undef Analytic
#undef Analytic_rep
#undef Assume
#undef Scalar
#undef Radius
#undef Series
#undef Series_rep
#undef Polynomial
#undef R
} // namespace mmx
#endif // __MMX_ANALYTIC_META_HPP

---