/******************************************************************************
* MODULE     : series_test.cpp
* DESCRIPTION: Test series
* COPYRIGHT  : (C) 2004  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.
******************************************************************************/

#include <basix/double.hpp>
#include <algebramix/series.hpp>
#include <algebramix/series_relaxed.hpp>
#include <algebramix/series_vector.hpp>
#include <algebramix/series_matrix.hpp>
#include <algebramix/series_elementary.hpp>
#include <algebramix/series_integer.hpp>
#include <algebramix/series_rational.hpp>

using namespace mmx;

template<typename C,typename V> void
arithmetic_simple_test () {
  typedef series<C,V> Series;

  Series z (1, 1);
  Series f (5, 3);
  Series g= 1 / (1 - z);
  Series h= seq (C (3), C (6), C (2));
  
  mmout << "z\t= " << z << "\n";
  mmout << "f\t= " << f << "\n";
  mmout << "g\t= " << g << "\n";
  mmout << "h\t= " << h << "\n";
  mmout << "-f\t= " << -f << "\n";
  mmout << "f+g\t= " << f+g << "\n";
  mmout << "f-g\t= " << f-g << "\n";
  mmout << "g-f\t= " << g-f << "\n";
  mmout << "f*g\t= " << f*g << "\n";
  mmout << "2*f\t= " << 2*f << "\n";
  mmout << "f*2\t= " << f*2 << "\n";
  mmout << "f/g\t= " << f/g << "\n";
  mmout << "f/3\t= " << f/3 << "\n";
  mmout << "f+=g\t= " << (f+=g) << "\n";
  mmout << "f-=g\t= " << (f-=g) << "\n";
  mmout << "f*=g\t= " << (f*=g) << "\n";
  mmout << "f/=g\t= " << (f/=g) << "\n";
  mmout << "f*=11\t= " << (f*=11) << "\n";
  mmout << "f/=11\t= " << (f/=11) << "\n";
  mmout << "f <* 2\t= " << lshiftz (f, 2) << "\n";
  mmout << "f >* 2\t= " << lshiftz (f, -2) << "\n";
  mmout << "g <* 2\t= " << lshiftz (g, 2) << "\n";
  mmout << "g >* 2\t= " << lshiftz (g, -2) << "\n";
  mmout << "g o z^3\t= " << dilate (g, 3) << "\n";
  mmout << "g'\t= " << derive (g) << "\n";
  mmout << "`g\t= " << integrate (g) << "\n";
  mmout << "f=g\t= " << (f==g) << "\n";
  mmout << "f!=g\t= " << (f!=g) << "\n";
  mmout << "f<g\t= " << (f<g) << "\n";
  mmout << "f<=g\t= " << (f<=g) << "\n";
  mmout << "f>g\t= " << (f>g) << "\n";
  mmout << "f>=g\t= " << (f>=g) << "\n\n";
}

template<typename C,typename V> void
elementary_simple_test () {
  typedef series<C,V> Series;
  Series z (1, 1);
  Series f (5, 3);

  mmout << "sqrt (z sin z)\t= " << sqrt (z * sin (z)) << "\n";
  mmout << "1/(1+f)\t\t= " << 1/(1+f) << "\n";
  mmout << "exp(f)\t\t= " << exp (f) << "\n";
  mmout << "log(1+f)\t= " << log (1+f) << "\n";
  mmout << "log_init(1+f,3)\t= " << log_init (1+f, C(3)) << "\n";
  mmout << "exp(z*exp(f))\t= " << exp (z * exp (f)) << "\n";
  mmout << "sqrt(1+f)\t= " << sqrt (1+f) << "\n";
  mmout << "(1+f)^7\t\t= " << pow (1+f, C(7)) << "\n";
  mmout << "cos_sin(f)\t= " << cos_sin (f) << "\n";
  mmout << "sin(f)\t\t= " << sin (f) << "\n";
  mmout << "cos(f)\t\t= " << cos (f) << "\n";
  mmout << "tan(f)\t\t= " << tan (f) << "\n\n";
}

template<typename C,typename V> void
division_simple_test () {
  typedef series<C,V> Series;
  Series zero (0);
  Series one (1);
  Series z (1, 1);
  Series a= one / (one - z);
  Series b= a * a;
  
  mmout << "a\t\t= " << a << "\n";
  mmout << "b\t\t= " << b << "\n";
  mmout << "b quo z\t\t= " << quo (b, z) << "\n";
  mmout << "b quo (2 + z)\t= " << quo (b, 2 + z) << "\n";
  mmout << "b quo 0\t\t= " << quo (b, zero) << "\n";
  mmout << "b rem z\t\t= " << rem (b, z) << "\n";
  mmout << "b rem (2 + z)\t= " << rem (b, 2 + z) << "\n";
  mmout << "b rem 0\t\t= " << rem (b, zero) << "\n";
}

template<typename C,typename V> void
solve_simple_test () {
  typedef series<C,V> Series;
  typedef vector_fixed<vector_naive,fixed_value<nat,2> > VV;
  Series z (1, 1);
  Series f (5, 3);
  vector<C,VV> vc (C(0), 2); vc[0]= C(1);
  vector<Series,VV> vv (C(0), 2); vv[0]= z; vv[1]= f;
  vector<Series,VV> v= solve_lde_init (vv, vc);
  mmout << "v\t= " << v << "\n";
  mmout << "ser v\t= " << as_series (v) << "\n";

  matrix<Series > m (Series (0), 2, 2);
  m (0, 0)= Series (1);
  m (1, 1)= Series (z);
  mmout << "m\t= " << m << "\n";
  mmout << "ser m\t= " << as_series (m) << "\n";
  mmout << "slde m\t= " << solve_lde (m) << "\n";
}

template<typename C, typename V>
class polya_series_rep : public recursive_series_rep<C,V> {
public:
  polya_series_rep ():
    recursive_series_rep<C,V> (format<C> ()) {}
  syntactic expression (const syntactic& z) const {
    return apply ("polya", z); }
  series<C,V> initialize () {
    series<C,V> f= this->me ();
    this->initial (0)= C(1);
    return lshiftz (f*f*f + 2*dilate (f, 3), 1) / 3; }
};

template<typename C,typename V> inline series<C,V>
polya_series () {
  series_rep<C,V>* rep= new polya_series_rep<C,V> ();
  return recursive (series<C,V> (rep));
}

template<typename C,typename V> void
polya_simple_test () {
  typedef series<C,V> Series;
  Series::set_output_order (100);
  mmout << "polya\t= " << polya_series<C,V> () << "\n";
}

int
main () {
  try {
    arithmetic_simple_test<rational,series_naive> ();
    elementary_simple_test<double,series_naive> ();
    division_simple_test<integer,series_naive> ();
    solve_simple_test<rational,series_naive> ();
    polya_simple_test<double,series_naive> ();
  }
  catch (exception& e) {
    mmerr << e << "\n";
    return -1;
  }
}