/Users/mourrain/Devel/mmx/algebramix/include/algebramix/series_relaxed.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_relaxed.hpp
* DESCRIPTION: Relaxed univariate power series
* COPYRIGHT  : (C) 2006  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_RELAXED__HPP
#define __MMX__SERIES_RELAXED__HPP
#include <algebramix/series_naive.hpp>
namespace mmx {
#define TMPL template<typename C,typename V>
#define Series series<C,V>
#define Series_rep series_rep<C,V>

/******************************************************************************
* Variant
******************************************************************************/

template<typename V>
struct series_relaxed {};

template<typename F, typename V, typename W>
struct implementation<F,V,series_relaxed<W> >:
    public implementation<F,V,W> {};

/******************************************************************************
* Relaxed multiplication
******************************************************************************/

template<typename U,typename W>
struct implementation<series_multiply,U,series_relaxed<W> >:
  public implementation<series_abstractions,U>
{

TMPL
class mul_series_rep: public Series_rep {
public:
  typedef typename series_polynomial_helper<C,V >::PV PV;
  typedef implementation<polynomial_multiply,PV> Pol;

protected:
  Series f, g;

  // NOTE: The following fields keep track of large segments of zeros
  // in the coefficients. This makes the multiplication algorithm linear
  // in the case when we multiply by a polynomial. This is for instance
  // interesting when computing coefficients of holonomic functions.
  int lnz_f;    // last non-zero coefficient in f
  int lnz_g;    // last non-zero coefficient in g
  nat zeros_f;  // bit 2^p set if f[2^p-1 ... 2^{p+1}-1] = 0
  nat zeros_g;  // bit 2^p set if g[2^p-1 ... 2^{p+1}-1] = 0

public:
  mul_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f (f2), g (g2), lnz_f (-1), lnz_g (-1),
    zeros_f (0), zeros_g (0) {}
  ~mul_series_rep () {}

  syntactic expression (const syntactic& z) const {
    return flatten (f, z) * flatten (g, z); }

  void Set_order (nat l2) {
    // NOTE: far from optimal, but let us not bother...
    if (l2 <= this->l) return;
    mmx_delete<C> (this->a, this->l);
    this->a= mmx_new<C> (l2);
    this->l= l2;
    Pol::clear (this->a, this->l);
    nat k= this->n;
    this->n= 0;
    lnz_f= -1;
    lnz_g= -1;
    /* not so nice in case of zero coefficients...
    if (k>0) {
      nat i= 1;
      while (((i<<1) + 1 < k) && ((i<<1) - 1 < this->l))
        i= (i<<1) + 1;
      Pol::mul (this->a, f->a, g->a, i, i);
      this->n= i;
    }
    */
    while (this->n < k) {
      (void) next ();
      this->n++;
    }
  }

  void Increase_order (nat l) {
    if (this->l < l) {
      nat k=1;
      while (k < (this->l << 1)) k= (k<<1) + 1;
      Set_order (max (k, l));
    }
    increase_order (f, l);
    increase_order (g, l);
  }

  C next () {
    nat k= this->n;
    if (exact_neq (f[k], C(0))) lnz_f= k;
    if (exact_neq (g[k], C(0))) lnz_g= k;
    if (k == 0) {
      if (lnz_f < 0) zeros_f |= 1;
      if (lnz_g < 0) zeros_g |= 1;
      if ((zeros_f&1) == 0 && (zeros_g&1) == 0) this->a[0]= f[0] * g[0];
    }
    else {
      if ((zeros_f&1) == 0) this->a[k] += f[0] * g[k];
      if ((zeros_g&1) == 0) this->a[k] += f[k] * g[0];
      if ((k&1) == 0) {
        nat p= 0;
        k += 2;
        while ((k&1) == 0 && k != 2) {
          p++;
          k= k >> 1;
          nat len= min (this->l, ((k+2)<<p) - 3) - ((k<<p) - 2);
          nat ord= min (((nat) 1)<<p, len);
          if (k == 2) {
            if (lnz_f < ((int) (1<<p)-1)) zeros_f |= (1<<p);
            if (lnz_g < ((int) (1<<p)-1)) zeros_g |= (1<<p);
          }
          C* aux= mmx_new<C> ((ord << 1) - 1);
          if ((zeros_f & (1<<p)) == 0 && lnz_g >= ((int) (((k-1)<<p)-1))) {
            //mmout << "  " << ((1<<p)-1) << " - " << (((1<<p)-1)+ord)
            //<< " x " << ((k-1)<<p)-1 << " - " << ((k-1)<<p)-1+ord << "\n";
            Pol::mul (aux, f->a + (1<<p)-1, g->a + ((k-1)<<p)-1, ord, ord);
            Pol::add (this->a + (k<<p)-2, aux, len);
          }
          if (k != 2)
            if ((zeros_g & (1<<p)) == 0 && lnz_f >= ((int) (((k-1)<<p)-1))) {
              //mmout << "  " << ((k-1)<<p)-1 << " - " << ((k-1)<<p)-1+ord
              //<< " x " << ((1<<p)-1) << " - " << (((1<<p)-1)+ord) << "\n";
              Pol::mul (aux, f->a + ((k-1)<<p)-1, g->a + (1<<p)-1, ord, ord);
              Pol::add (this->a + (k<<p)-2, aux, len);
            }
          mmx_delete<C> (aux, (ord << 1) - 1);
        }
      }
    }
    /*
    mmout << this->n << ", " << type_name<C> () << "\n\t";
    for (nat i= 0; i<=this->n; i++) mmout << f->a[i] << " ";
    Series::set_formula_output (true);
    mmout << "; " << zeros_f << "; " << flatten (f, "z") << "\n";
    Series::set_formula_output (false);
    mmout << "\t";
    for (nat i= 0; i<=this->n; i++) mmout << g->a[i] << " ";
    Series::set_formula_output (true);
    mmout << "; " << zeros_g << "; " << flatten (g, "z") << "\n";
    Series::set_formula_output (false);
    mmout << "\t";
    for (nat i= 0; i<=this->n; i++) mmout << this->a[i] << " ";
    mmout << "\n";
    */
    return this->a [this->n];
  }
};

TMPL static inline Series
ser_mul (const Series& f, const Series& g) {
  typedef mul_series_rep<C,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g))
    return Series (CF(f));
  return (Series_rep*) new Mul_rep (f, g); }

TMPL static inline Series
ser_truncate_mul (const Series& f, const Series& g, nat nf, nat ng) {
  typedef mul_series_rep<C,V> Mul_rep;
  if (is_exact_zero (f) || is_exact_zero (g) || nf == 0 || ng == 0)
    return Series (CF(f));
  return (Series_rep*)
    new Mul_rep (piecewise (f, Series (CF(f)), nf),
                 piecewise (g, Series (CF(g)), ng)); }

}; // implementation<series_multiply,V,series_relaxed<W> >:

#ifdef ALGEBRAMIX_ENABLE_UNSTABLE

/******************************************************************************
* Right composition with polynomials
******************************************************************************/

TMPL Series compose (const Series& f, nat i, nat p, nat q, Series* H);

TMPL
struct partial_pol_compose_series_rep: public Series_rep {
  Series f, h;
  nat i, r;
  partial_pol_compose_series_rep<C>
    (const Series& f2, nat i2, nat p, nat q, Series* H);
  generic expression (const generic& z) const {
    return gen ("partial_polynomial_compose"); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (min (l, r));
    increase_order (h, min (l, r));
  }
  C next ();
};

TMPL
partial_pol_compose_series_rep<C>::partial_pol_compose_series_rep (
  const Series& f2, nat i2, nat p, nat q, Series* H):
    Series_rep (CF(f2)), f (f2), i (i2), r ((p-1) * (q-1) + 1)
{
  assert (p > 1);
  Series h1= compose (f, i, p>>1, q, H);
  Series h2= compose (f, i+ (p>>1), p>>1, q, H);
  int d= p>>1;
  h= h1+ lshiftz (h2 * rshiftz (H[p>>1], d), d);
}

TMPL C
partial_pol_compose_series_rep<C>::next () {
  nat n= this->n;
  if (n == r) h= 0;
  if (n >= r) return 0;
  return h[n];
}

TMPL Series
compose (const Series& f, nat i, nat p, nat q, Series* H) {
  if (p == 1) return Series (f[i]);
  return (Series_rep*) new partial_pol_compose_series_rep<C> (f, i, p, q, H);
}

TMPL
struct pol_compose_series_rep: public Series_rep {
  Series f, *H, h;
  nat q;
  pol_compose_series_rep<C> (const Series& f, const Series& g, nat q);
  ~pol_compose_series_rep<C> () { mmx_classical_delete<Series > (H); }
  syntactic expression (const syntactic& z) const {
    return apply (GEN_COMPOSE, flatten (f, z), flatten (H[1], z)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l);
    increase_order (h, l);
  }
  C next ();
};

TMPL
pol_compose_series_rep<C>::pol_compose_series_rep
  (const Series& f2, const Series& g, nat q2):
    Series_rep (CF(f2)), f (f2), q (q2)
{
  H= mmx_classical_new<Series > (2);
  H[1]= g;
  h= compose (f, 0, 1, q, H);
}

TMPL C
pol_compose_series_rep<C>::next () {
  nat n= this->n;
  nat k= n;
  nat i;
  if (k > 0) while ((k&1) == 0) k= k>>1;
  if (k == 1) {
    if (n > 1) {
      Series* H2= mmx_classical_new<Series > (n+1);
      for (i=1; i<n; i=i<<1) H2[i]= H[i];
      H2[n]= H[n>>1] * H[n>>1];
      mmx_classical_delete<Series > (H);
      H= H2;
    }
    h= compose (f, 0, n<<1, q, H);
    increase_order (h, this->l);
  }
  return h[n];
}

TMPL Series
compose (const Series& f, const Series& g, nat q) {
  if (is_exact_zero (f)) return Series (CF(f));
  ASSERT (g[0] == 0, "constant coefficient must be zero");  
  return (Series_rep*) new pol_compose_series_rep<C> (f, g, q);
}

/******************************************************************************
* General relaxed composition
******************************************************************************/

template<typename C> C
multiply_range (nat first, nat end) {
  if (end == first) return C(1);
  else if (end == first+1) return C (first);
  nat half= (first + end) >> 1;
  return multiply_range<C> (first, half) * multiply_range<C> (half, end);
}

TMPL
class relaxed_compose_series_rep: public Series_rep {
protected:
  Series f, g, h;

public:
  relaxed_compose_series_rep (const Series& f2, const Series& g2):
    Series_rep (CF(f2)), f (f2), g (g2),
    h ((Series_rep*) new compose_series_rep<C,V> (f2, g2)) {}
  ~relaxed_compose_series_rep () {}

  syntactic expression (const syntactic& z) const {
    return apply (GEN_COMPOSE, flatten (f, z), flatten (g, z)); }

  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l);
    increase_order (g, l);
  }

  C next () {
    nat n= this->n;
    if (n==0) return f[0];
    if (n==1) return f[1]*g[1];
    if (n==2) return f[2]*g[1]*g[1]+ f[1]*g[2];
    if (n==3) return (f[3]*g[1]*g[1]+ C(2)*f[2]*g[2])*g[1]+ f[1]*g[3];

    nat p=0;
    while (((p+1) * (1<<(2*p))) <= n) p++;
    if ((n != 4) && ((p * (1<<(2*p-2))) <= (n-1))) return h[n];

    nat nn= (p+1) * (1<<(2*p));
    nat i, q= 2<<p, r= (nn+q-1)/q;
    Series g1 = restrict (g, 0, q);
    Series dg1= derive (g1);
    Series g2 = rshiftz (g - g1, (int) q);

    Series D= f;
    for (i=1; i<r; i++) D= derive (D);
    D= compose (D, g1, q) / multiply_range<C> (1, r); // divides by (r-1)!
    h= D;

    for (i=r-1; i>0; i--) {
      D= (integrate (C(i) * D * dg1)) + f[i-1];
      h= D + lshiftz (h*g2, q);
    }

    increase_order (h, this->l);
    return h[n];
  }
};

TMPL Series
bk_compose (const Series& f, const Series& g) {
  // NOTE: not yet default; needs further efficiency and accuracy twiddling
  if (is_exact_zero (f)) return Series (CF(f));
  ASSERT (g[0] == 0, "constant coefficient must be zero");  
  return (Series_rep*) new relaxed_compose_series_rep<C> (f, g);
}
#endif // UNSTABLE

#undef TMPL
#undef Series
#undef Series_rep
} // namespace mmx
#endif // __MMX__SERIES_RELAXED__HPP