/******************************************************************************
* MODULE     : ball_complex_test.cpp
* DESCRIPTION: Test complex balls
* 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.
******************************************************************************/

#include <numerix/ball_complex.hpp>
//#include <numerix/complex_interval.hpp>
using namespace mmx;

#define R floating<>
#define C complex<R >
#define Real_ball ball<R >
#define Ball ball<C >
#define Radius unary_return_type_helper<radius_op,Ball >::RET
//#define Real_ball interval<R >
//#define Ball complex<Real_ball >

void
test_add_simple (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Additive arithmetic at precision " << prec << "\n\n";
  Ball x= gaussian (Real_ball (1), Real_ball (2));
  Ball y= 7;
  mmout << "x\t= " << x << "\n";
  mmout << "y\t= " << y << "\n";
  mmout << "-x\t= " << (-x) << "\n";
  mmout << "x+x\t= " << (x+x) << "\n";
  mmout << "x+y\t= " << (x+y) << "\n";
  mmout << "x-x\t= " << (x-x) << "\n";
  mmout << "x-y\t= " << (x-y) << "\n";
  mmout << "(x-x)+y\t= " << ((x-x)+y) << "\n";
  mmout << "\n";
}

void
test_mul_simple (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Multiplicative arithmetic at precision " << prec << "\n\n";
  Ball x= Ball (2) + Ball (1);
  mmout << "x\t= " << x << "\n";
  mmout << "x * x\t= " << x*x << "\n";
  mmout << "x * -x\t= " << x*(-x) << "\n";
  mmout << "-x * x\t= " << (-x)*x << "\n";
  mmout << "-x * -x\t= " << (-x)*(-x) << "\n";
  mmout << "x / x\t= " << (x/x) << "\n";
  mmout << "x / -x\t= " << (x/(-x)) << "\n";
  mmout << "-x / x\t= " << ((-x)/x) << "\n";
  mmout << "-x / -x\t= " << ((-x)/(-x)) << "\n";
  mmout << "x*x / x\t= " << ((x*x)/x) << "\n";
  mmout << "\n";
}

void
test_other_ops (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Other operations " << prec << "\n\n";
  Ball x= gaussian (Real_ball (2), Real_ball (3));
  Ball y= make_ball<Ball > (C (-1), Radius (5) >> 1);
  mmout << "x\t= " << x << "\n";
  mmout << "y\t= " << y << "\n";
  mmout << "x << 5\t= " << (x << 5) << "\n";
  mmout << "x >> 3\t= " << (x >> 3) << "\n";
  mmout << "y << 5\t= " << (y << 5) << "\n";
  mmout << "abs(x)\t= " << (abs (x)) << "\n";
  mmout << "abs(-x)\t= " << (abs (-x)) << "\n";
  mmout << "abs(y)\t= " << (abs (y)) << "\n";
  mmout << "arg(x)\t= " << (arg (x)) << "\n";
  mmout << "arg(-x)\t= " << (arg (-x)) << "\n";
  mmout << "arg(y)\t= " << (arg (y)) << "\n";
  mmout << "\n";
}

void
test_mul_progression (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Multiplicative progression at precision " << prec << "\n\n";
  int i;
  Ball x= gaussian (Real_ball (1), Real_ball (2));
  Ball y= x;
  for (i=0; i<50; i++) {
    mmout << i << "\t" << y << "\n";
    y *= x;
  }
  mmout << "\n";
}

void
test_mul_overflow (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Multiplicative overflows at precision " << prec << "\n\n";
  int i;
  Ball x= gaussian (Real_ball (1), Real_ball (2));
  for (i=0; i<50; i++) {
    mmout << i << "\t" << x << "\n";
    x *= x;
  }
  mmout << "\n";
}

void
test_mul_underflow (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Multiplicative underflows at precision " << prec << "\n\n";
  int i;
  Ball x= gaussian (Real_ball (0.3));
  // Ball x= gaussian (Real_ball (0.3), Real_ball (0.1));
  // FIXME: bug in MPFR when replacing by uncommented line?
  for (i=0; i<50; i++) {
    mmout << i << "\t" << x << "\n";
    x *= x;
  }
  mmout << "\n";
}

void
test_elementary_funcs (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Elementary functions at precision " << prec << "\n\n";
  Ball x= 1;
  Ball y= gaussian (Real_ball (3e6));
  // Ball y= gaussian (Real_ball (3e20));
  // FIXME: bug in MPFR when replacing by uncommented line?
  Ball z= x - x/y;
  mmout << "x\t\t= " << x << "\n";
  mmout << "y\t\t= " << y << "\n";
  mmout << "z\t\t= " << z << "\n";
  mmout << "pi\t\t= " << Pi(Ball) << "\n";
  mmout << "i\t\t= " << Imaginary(Ball) << "\n";
  mmout << "NaN\t\t= " << Nan(Ball) << "\n";
  mmout << "Infty\t\t= " << Infinity(Ball) << "\n";
  mmout << "Fuzz\t\t= " << Fuzz(Ball) << "\n";
  mmout << "sqrt (x)\t= " << sqrt (x) << "\n";
  mmout << "sqrt (y)\t= " << sqrt (y) << "\n";
  mmout << "exp (x)\t\t= " << exp (x) << "\n";
  mmout << "exp (y)\t\t= " << exp (y) << "\n";
  mmout << "exp (-x)\t= " << exp (-x) << "\n";
  mmout << "exp (-y)\t= " << exp (-y) << "\n";
  mmout << "log (x)\t\t= " << log (x) << "\n";
  mmout << "log (y)\t\t= " << log (y) << "\n";
  mmout << "cos (x)\t\t= " << cos (x) << "\n";
  mmout << "cos (y)\t\t= " << cos (y) << "\n";
  mmout << "sin (x)\t\t= " << sin (x) << "\n";
  mmout << "sin (y)\t\t= " << sin (y) << "\n";
  mmout << "tan (x)\t\t= " << tan (x) << "\n";
  mmout << "tan (y)\t\t= " << tan (y) << "\n";
  mmout << "arc sin (1/2)\t= " << asin (x>>1) << "\n";
  mmout << "arc sin (z)\t= " << asin (z) << "\n";
  mmout << "arc cos (1/2)\t= " << acos (x>>1) << "\n";
  mmout << "arc tan (x)\t= " << atan (x) << "\n";
  mmout << "arc tan (y)\t= " << atan (y) << "\n";
  mmout << "arc tan (z)\t= " << atan (z) << "\n";
  mmout << "\n";
}

void
test_all (int prec) {
  test_add_simple (prec);
  test_mul_simple (prec);
  test_other_ops (prec);
  test_mul_progression (prec);
  test_mul_overflow (prec);
  test_mul_underflow (prec);
  test_elementary_funcs (prec);
}

int
main () {
  test_all (30);
  test_all (100);
  test_all (200);

  mmout << "infty\t\t= " << Infinity (Ball) << "\n";
  mmout << "-infty\t\t= " << -Infinity (Ball) << "\n";
  mmout << "infty+1\t\t= " << Infinity (Ball) + Ball (1) << "\n";
  mmout << "1-infty\t\t= " << Ball (1) - Infinity (Ball) << "\n";
  mmout << "2 infty\t\t= " << Ball (2) * Infinity (Ball) << "\n";
  mmout << "1/infty\t\t= " << Ball (1) / Infinity (Ball) << "\n";
  mmout << "infty^2\t\t= " << Infinity (Ball) * Infinity (Ball) << "\n";
  mmout << "infty/infty\t= " << Infinity (Ball) / Infinity (Ball) << "\n";
  mmout << "|infty|\t\t= " << abs (Infinity (Ball)) << "\n";
  mmout << "center(infty)\t= " << center (Infinity (Ball)) << "\n";
  mmout << "radius(infty)\t= " << radius (Infinity (Ball)) << "\n";
  mmout << "infty << 3\t= " << (Infinity (Ball) << 3) << "\n";
  mmout << "infty >> 3\t= " << (Infinity (Ball) >> 3) << "\n";
  return 0;
}