include/algebramix/series_implicit.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_implicit.hpp
* DESCRIPTION: Implicitly defined power series
* 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_SERIES_IMPLICIT_HPP
#define __MMX_SERIES_IMPLICIT_HPP
#include <algebramix/series_vector.hpp>
namespace mmx {
#define TMPL_DEF template<typename C, typename V=typename Series_variant(C)>
#define TMPL template<typename C, typename V>
#define Series series<C,V>
#define Series_rep series_rep<C,V>
#define VC vector<C>
#define VSeries series<vector<C>,V>
#define VSeries_rep series_rep<vector<C>,V>
#define UC unknown<C,V>
#define UC_rep unknown_rep<C,V>
#define USeries series<UC >
#define USeries_rep series_rep<UC >
#define Solver_rep solver_series_rep<C,V>
TMPL class solver_series_rep;

/******************************************************************************
* Unknown coefficients
******************************************************************************/

TMPL_DEF
class unknown_rep REP_STRUCT {
public:
  Solver_rep* f;   // underlying series
  C   b;           // base
  C*  s;           // slope
  nat l;           // number of allocated elements
  nat i1;          // starting coefficient
  nat i2;          // ending coefficient

  inline void normalize () {
    nat k=i1;
    while (i1 < i2 && is_exact_zero (s[i2-k-1])) i2--;
    while (i1 < i2 && is_exact_zero (s[i1-k])) i1++;
    if (i1 != k) {
      //mmerr << "    normalize " << i1 << " -- " << i2 << "\n";
      nat d= i1-k;
      C* s2= mmx_new<C> (i2-i1);
      for (nat i=0; i<i2-i1; i++) s2[i]= s[i+d];
      mmx_delete<C> (s, l);
      s= s2;
      l= i2-i1;
    }
  }
    
  inline unknown_rep ():
    f (NULL), b (), s (mmx_new<C> (0)), l (0), i1 (0), i2 (0) {}
  inline unknown_rep (const C& b2):
    f (NULL), b (b2), s (mmx_new<C> (0)), l (0), i1 (0), i2 (0) {}
  unknown_rep (Solver_rep* f2, const C& b2, C* ss, nat s2, nat e2):
    f (f2), b (b2), s (ss), l (e2 - s2), i1 (s2), i2 (e2) { normalize (); }
  inline virtual ~unknown_rep () {
    mmx_delete<C> (s, l); }
};

TMPL_DEF
class unknown {
INDIRECT_PROTO_2 (unknown, unknown_rep, C, V)
public:
  inline unknown ():
    rep (new UC_rep ()) {}
  template<typename T> inline unknown (const T& b):
    rep (new UC_rep (as<C> (b))) {}
  inline unknown (Solver_rep* f, const C& b, C* s, nat i1, nat i2):
    rep (new UC_rep (f, b, s, i1, i2)) {}
};
INDIRECT_IMPL_2 (unknown, unknown_rep, typename C, C, typename V, V)

TMPL inline bool
is_known (const UC& c) {
  return c->i2 == c->i1;
}

TMPL inline C
known (const UC& c) {
  ASSERT (c->i2 == c->i1, "cast failed");
  return c->b;
}

/******************************************************************************
* Prototype for implicit series
******************************************************************************/

TMPL_DEF
class solver_series_rep: public VSeries_rep {
public:
  nat m;
  vector<USeries > eqs;
  nat cur_n;
  vector<UC> sys;

public:
  inline solver_series_rep (nat m, const vector<vector<C> >& init);
  virtual void Increase_order (nat l);
  virtual syntactic name_component (nat i);
  virtual vector<USeries > initialize () = 0;
  vector<USeries > me ();
  VC next ();
};

/******************************************************************************
* Mandatory routines for unknown coefficients
******************************************************************************/

TMPL bool
is_exact_zero (const UC c) {
  return c->i1 == c->i2 && is_exact_zero (c->b);
}

template<typename Op, typename C, typename V> nat
unary_hash (const UC& c) {
  register nat i, h= 78460;
  if (c->i1 == c->i2) return Op::op (c->b) ^ h;
  h += (c->i1 << 3) ^ Op::op (c->b);
  for (i=0; i<c->i2 - c->i1; i++)
    h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (c->s[i]);
  return h;
}

template<typename Op,typename C,typename V> inline bool
binary_test (const UC& c1, const UC& c2) {
  if (Op::not_op (c1->b, c2->b)) return false;
  if (c1->i1 == c1->i2 || c2->i1 == c2->i2)
    return c1->i1 == c1->i2 && c2->i1 == c2->i2;
  if (Op::not_op (c1->f, c2->f)) return false;
  if (Op::not_op (c1->i1, c2->i1)) return false;
  if (Op::not_op (c1->i2, c2->i2)) return false;
  for (nat i= c1->i1; i<c1->i2; i++)
    if (Op::not_op (c1->s[i - c1->i1], c2->s[i - c1->i1]))
      return false;
  return true;
}

TMPL syntactic
flatten (const UC& c) {
  syntactic sum= flatten (c->b);
  for (nat i=c->i1; i<c->i2; i++) {
    nat k= i / c->f->m;
    nat j= i % c->f->m;
    sum= sum + flatten (c->s[i-c->i1]) *
               access (c->f->name_component (j), flatten (k));
  }
  return sum;
}

TRUE_IDENTITY_OP_SUGAR(TMPL,UC)
EXACT_IDENTITY_OP_SUGAR(TMPL,UC)

/******************************************************************************
* Basic arithmetic on affine combinations of unknown coefficients
******************************************************************************/

TMPL UC
operator - (const UC& c) {
  //mmerr << "    negate " << c << "\n";
  if (is_exact_zero (c)) return c;
  nat n= c->i2 - c->i1;
  C* s= mmx_new<C> (n);
  for (nat i=0; i<n; i++)
    s[i]= -c->s[i];
  return UC (c->f, -c->b, s, c->i1, c->i2);
}

TMPL UC
operator + (const UC& c1, const C& c2) {
  //mmerr << "    scalar plus " << c1 << ", " << c2 << "\n";
  if (is_exact_zero (c1))
    return UC (c1->f, c2, mmx_new<C> (0), c1->i1, c1->i1);
  if (is_exact_zero (c2)) return c1;
  nat n= c1->i2 - c1->i1;
  C* s= mmx_new<C> (n);
  for (nat i=0; i<n; i++)
    s[i]= c1->s[i];
  return UC (c1->f, c1->b + c2, s, c1->i1, c1->i2);
}

TMPL UC
operator * (const UC& c1, const C& c2) {
  //mmerr << "    scalar times " << c1 << ", " << c2 << "\n";
  if (is_exact_zero (c1)) return c1;
  if (is_exact_zero (c2))
    return UC (c1->f, C (0), mmx_new<C> (0), c1->i1, c1->i1);
  nat n= c1->i2 - c1->i1;
  C* s= mmx_new<C> (n);
  for (nat i=0; i<n; i++)
    s[i]= c1->s[i] * c2;
  return UC (c1->f, c1->b * c2, s, c1->i1, c1->i2);
}

TMPL UC
operator * (const C& c1, const UC& c2) {
  //mmerr << "    scalar times " << c1 << ", " << c2 << "\n";
  if (is_exact_zero (c2)) return c2;
  if (is_exact_zero (c1))
    return UC (c2->f, C (0), mmx_new<C> (0), c2->i1, c2->i1);
  nat n= c2->i2 - c2->i1;
  C* s= mmx_new<C> (n);
  for (nat i=0; i<n; i++)
    s[i]= c1 * c2->s[i];
  return UC (c2->f, c1 * c2->b, s, c2->i1, c2->i2);
}

TMPL UC
operator + (const UC& c1, const UC& c2) {
  //mmerr << "    plus " << c1 << ", " << c2 << "\n";
  if (c1->i1 == c1->i2) return c2 + known (c1);
  if (c2->i1 == c2->i2) return c1 + known (c2);
  ASSERT (c1->f == c2->f, "incompatible unknown coefficients");
  nat i1= min (c1->i1, c2->i1);
  nat i2= max (c1->i2, c2->i2);
  C* s= mmx_new<C> (i2-i1);
  for (nat i= i1; i<i2; i++)
    s[i-i1]=
      (i >= c1->i1 && i < c1->i2? c1->s[i - c1->i1]: C(0)) +
      (i >= c2->i1 && i < c2->i2? c2->s[i - c2->i1]: C(0));
  return UC (c1->f, c1->b + c2->b, s, i1, i2);
}

TMPL UC
operator - (const UC& c1, const UC& c2) {
  //mmerr << "    minus " << c1 << ", " << c2 << "\n";
  if (c1->i1 == c1->i2) return (-c2) + known (c1);
  if (c2->i1 == c2->i2) return c1 + (-known (c2));
  ASSERT (c1->f == c2->f, "incompatible unknown coefficients");
  nat i1= min (c1->i1, c2->i1);
  nat i2= max (c1->i2, c2->i2);
  C* s= mmx_new<C> (i2-i1);
  for (nat i=i1; i<i2; i++)
    s[i-i1]=
      (i >= c1->i1 && i < c1->i2? c1->s[i - c1->i1]: C(0)) -
      (i >= c2->i1 && i < c2->i2? c2->s[i - c2->i1]: C(0));
  return UC (c1->f, c1->b - c2->b, s, i1, i2);
}

TMPL UC
substitute (const UC& c) {
  //mmerr << "    substitute " << c << "\n";
  if (c->i1 == c->i2 || c->i1 >= c->f->n * c->f->m) return c;
  nat i1= min (c->f->n * c->f->m, c->i2);
  nat d= i1 - c->i1;
  nat n= c->i2 - i1;
  C*  s= mmx_new<C> (n);
  C   b= c->b;
  for (nat i=0; i<d; i++) {
    nat k= (i + c->i1) / c->f->m;
    nat j= (i + c->i1) % c->f->m;
    b += c->s[i] * c->f->a[k][j];
  }
  for (nat i=0; i<n; i++)
    s[i]= c->s[i + d];
  return UC (c->f, b, s, i1, c->i2);
}

TMPL UC
operator * (const UC& c1, const UC& c2) {
  //mmerr << "    times " << c1 << ", " << c2 << "\n";
  if (is_exact_zero (c1)) return c1;
  if (is_exact_zero (c2)) return c2;
  if (c1->i1 == c1->i2) return known (c1) * c2;
  if (c2->i1 == c2->i2) return c1 * known (c2);
  UC c1b= substitute (c1);
  UC c2b= substitute (c2);
  if (c1b->i1 == c1b->i2) return known (c1b) * c2b;
  if (c2b->i1 == c2b->i2) return c1b * known (c2b);
  ERROR ("invalid product of unknown coefficients");
}

/******************************************************************************
* Incremental triangulation
******************************************************************************/

TMPL void
reduce (UC& c1, UC& c2) {
  // Triangulation of two "rows" c1 and c2
  // On exit, the simplest "row" is put into c1
  //mmerr << "    reduce " << c1 << ", " << c2 << "\n";
  if (is_exact_zero (c1)) return;
  if (is_exact_zero (c2)) { swap (c1, c2); return; }
  ASSERT (c1->f == c2->f, "incompatible unknown coefficients");
  if (c1->i2 < c2->i2) return;
  if (c2->i2 < c1->i2) { swap (c1, c2); return; }
  if (better_pivot (c1->s[c1->i2 - 1 - c1->i1],
                    c2->s[c2->i2 - 1 - c2->i1])) swap (c1, c2);
  C lambda= c1->s[c1->i2 - 1 - c1->i1] / c2->s[c2->i2 - 1 - c2->i1];
  nat i1= min (c1->i1, c2->i1);
  nat i2= c1->i2;
  C* s= mmx_new<C> (i2-i1-1);
  for (nat i= i1; i<i2-1; i++)
    s[i-i1]=
      (i >= c1->i1? c1->s[i - c1->i1]: C(0)) -
      lambda * (i >= c2->i1? c2->s[i - c2->i1]: C(0));
  c1= UC (c1->f, c1->b - lambda * c2->b, s, i1, i2-1);
}

TMPL void
insert_and_reduce (vector<UC >& sys, UC& c) {
  for (nat i=0; i<N(sys); i++)
    reduce (c, sys[i]);
  if (is_exact_zero (c)) return;
  ASSERT (c->i1 != c->i2, "contradictory equations");
  sys << c;
}

/******************************************************************************
* Unknown and known series
******************************************************************************/

TMPL_DEF
class unknown_series_rep: public USeries_rep {
  Solver_rep* f;
  nat k;
public:
  unknown_series_rep (Solver_rep* f2, nat k2):
    USeries_rep (format<UC > (/*FIXME*/)), f (f2), k (k2)
  {
    this->n= f->n;
    this->Set_order (this->n);
    for (nat i=0; i<f->n; i++) {
      C* s= mmx_new<C> (0);
      this->a[i]= UC (f, f->a[i][k], s, i, i);
    }
  }
  syntactic expression (const syntactic& z) const {
    return syn (f->name_component (k), z); }
  UC next () {
    C* s= mmx_new<C> (1);
    s[0]= C (1);
    return UC (f, C(0), s, this->n * f->m + k, this->n * f->m + k + 1); } 
};

TMPL USeries
unknown_series (Solver_rep* f, nat k) {
  return (USeries_rep*) new unknown_series_rep<C> (f, k);
}

template<typename C, typename V, typename UV>
class known_series_rep: public Series_rep {
public:
  series<UC,UV> f;
  inline known_series_rep (const series<UC,UV>& f2):
    Series_rep (format<C> (/*FIXME*/)), f (f2) {}
  inline syntactic expression (const syntactic& z) const {
    return syn ("known", flatten (f, z)); }
  inline void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  inline C next () { return known (substitute (f[this->n])); }
};

template<typename C, typename V, typename UV> Series
known (const series<UC,UV>& f) {
  return (Series_rep*) new known_series_rep<C,V,UV> (f);
}

/******************************************************************************
* Substitution product
******************************************************************************/

template<typename C, typename V, typename UV>
class subst_mul_series_rep:
   public implementation<series_abstractions,UV>
     ::template binary_series_rep<mul_op,UC,UV> {
public:
  nat f_sh, g_sh;
  Series f_kn, g_kn, inner;
  inline subst_mul_series_rep (const USeries& f, const USeries& g):
    implementation<series_abstractions,UV>
     ::template binary_series_rep<mul_op,UC,UV> (f, g),
    f_sh (1), g_sh (1), f_kn (known (f)), g_kn (known (g)),
    inner (rshiftz (f_kn, (int) f_sh) * rshiftz (g_kn, (int) g_sh)) {}
  UC next () {
    if (this->n >= g_sh &&
        !is_known (substitute (this->f [this->n - g_sh]))) {
      g_sh = g_sh << 1;
      inner= rshiftz (f_kn, (int) f_sh) * rshiftz (g_kn, (int) g_sh);
    }
    if (this->n >= f_sh &&
        !is_known (substitute (this->g [this->n - f_sh]))) {
      f_sh = f_sh << 1;
      inner= rshiftz (f_kn, (int) f_sh) * rshiftz (g_kn, (int) g_sh);
    }
    if (this->n < f_sh + g_sh) {
      UC acc= this->f[0] * this->g[this->n];
      for (nat i=1; i<=this->n; i++)
        acc += this->f[i] * this->g[this->n-i];
      return acc;
    }
    else {
      UC acc= f_kn[0] * this->g[this->n] + this->f[this->n] * g_kn[0];
      for (nat i=1; i<f_sh; i++)
        acc += f_kn[i] * this->g[this->n-i];
      for (nat i=1; i<g_sh; i++)
        acc += this->f[this->n-i] * g_kn[i];
      return acc + inner[this->n - f_sh - g_sh];
    }
  }
};

template<typename C, typename V, typename UV> series<UC,UV>
operator * (const series<UC,UV>& f, const series<UC,UV>& g) {
  if (is_exact_zero (f) || is_exact_zero (g))
    return series<UC,UV> (CF(f));
  return (series_rep<UC,UV>*) new subst_mul_series_rep<C,V,UV> (f, g);
}

/******************************************************************************
* Implicit series
******************************************************************************/

TMPL inline
Solver_rep::solver_series_rep (nat m2, const vector<vector<C> >& init):
  VSeries_rep (format<VC > (/*FIXME*/)), m (m2), cur_n (0)
{
  this->n= N(init);
  this->Set_order (this->n);
  for (nat i=0; i<this->n; i++)
    this->a[i]= init[i];
}

TMPL void
Solver_rep::Increase_order (nat l) {
  VSeries_rep::Increase_order (l);
  for (nat i=0; i<N(eqs); i++)
    increase_order (eqs[i], l);
}

TMPL syntactic
Solver_rep::name_component (nat i) {
  if (m == 1) return syntactic ("f");
  else return access (syntactic ("f"), flatten (i+1));
}

TMPL vector<USeries >
Solver_rep::me () {
  vector<USeries > r= fill<USeries > (m);
  for (nat i=0; i<m; i++)
    r[i]= unknown_series (this, i);
  return r;
}

TMPL VC
Solver_rep::next () {
  //mmerr << "Solving " << this->n << "\n";
  for (nat i=0; i<N(sys); i++)
    sys[i]= substitute (sys[i]);
  VC  ret = fill<C> (m);
  nat done= 0;
  while (true) {
    ASSERT (cur_n < this->n + 100, "too large delay in implicit solve");
    //mmerr << "  Coefficient " << cur_n << "\n";
    for (nat j=0; j<N(eqs); j++) {
      UC c= eqs[j][cur_n];
      //mmerr << "    Equation " << j << "= " << c << "\n";
      insert_and_reduce (sys, c);
    }
    //mmerr << "  System= " << sys << "\n";
    cur_n++;
    for (nat i=0; i<N(sys); i++)
      ASSERT (sys[i]->f == ((Solver_rep*) this) &&
              sys[i]->i1 >= this->n * m && sys[i]->i2 > this->n * m,
              "invalid situation during implicit solving");
    while (N(sys) > 0 && done < m) {
      UC c= sys[N(sys)-1];
      if (c->i2 <= this->n * m + done + 1) {
        nat j1 = c->i1 - this->n * m;
        nat j2 = min (done, c->i2 - this->n * m);
        C   rhs= c->b;
        for (nat j=j1; j<j2; j++)
          rhs += c->s[j-j1] * ret[j];
        if (c->i2 <= this->n * m + done) {
          ASSERT (is_exact_zero (rhs), "contradictory equations"); }
        else {
          ret[done]= -rhs / c->s[done-j1];
          //mmerr << "  Component " << done << "= " << ret[done] << "\n";
          done++;
        }
        sys.secure ();
        inside (sys) -> resize (N(sys) - 1);
      }
      else break;
    }
    if (done == m) return ret;
  }
}

TMPL
class solver_container_series_rep: public VSeries_rep {
  VSeries f;
public:
  solver_container_series_rep (const VSeries& f2):
    VSeries_rep (CF(f2)), f (f2) {
    Solver_rep* rep= (Solver_rep*) f.operator -> ();
    rep->eqs= rep->initialize (); }
  ~solver_container_series_rep () {
    Solver_rep* rep= (Solver_rep*) f.operator -> ();
    rep->eqs= vector<USeries > (); }
  syntactic expression (const syntactic& z) const {
    return flatten (f, z); }
  virtual void Increase_order (nat l) {
    VSeries_rep::Increase_order (l);
    increase_order (f, l); }
  VC next () { return f[this->n]; }
};

TMPL VSeries
solver (const VSeries& f) {
  return (VSeries_rep*) new solver_container_series_rep<C,V> (f);
}

#undef TMPL_DEF
#undef TMPL
#undef Series
#undef Series_rep
#undef VC
#undef VSeries
#undef VSeries_rep
#undef UC
#undef UC_rep
#undef Iseries
#undef Iseries_rep
#undef Solver_rep
} // namespace mmx
#endif // __MMX_SERIES_IMPLICIT_HPP

---