/******************************************************************************
* MODULE     : ball_test.cc
* DESCRIPTION: Test real 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.hpp>
//#include <numerix/interval.hpp>
using namespace mmx;

#define R floating<>
#define Ball ball<R >
//#define Ball interval<R >

void
test_constants () {
  mmout << "*** Remarkable constants\n\n";
  mmout << "smallest\t= " << Smallest (Ball) << "\n";
  mmout << "largest\t= " << Largest (Ball) << "\n";
  mmout << "fuzzy\t= " << Fuzz (Ball) << "\n";
  mmout << "\n";
}

void
test_add_simple (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Additive arithmetic at precision " << prec << "\n\n";
  Ball x= 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= Ball (2);
  Ball y= make_interval<Ball> (R (-3), R (2));
  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 << "\n";
}

void
test_mul_progression (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Multiplicative progression at precision " << prec << "\n\n";
  int i;
  Ball x= 3;
  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= 3;
  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= 0.3;
  for (i=0; i<50; i++) {
    mmout << i << "\t" << x << "\n";
    x *= x;
  }
  mmout << "\n";
}

void
test_muller (int prec) {
  mmx_bit_precision= prec;
  mmout << "*** Muller's sequence at precision " << prec << "\n\n";
  int i;
  Ball a= 99;
  Ball b= -600;
  Ball d= -7;
  Ball x= 4;
  for (i=0; i<50; i++) {
    mmout << i << "\t" << x << "\n";
    x= (a*x + b) / (x + d);
  }
  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= 3e20;
  Ball z= x - x/y;
  mmout << "log (2)\t\t= " << Log2 (Ball) << "\n";
  mmout << "pi\t\t= " << Pi (Ball) << "\n";
  mmout << "euler\t\t= " << Euler (Ball) << "\n";
  mmout << "NaN\t\t= " << Nan (Ball) << "\n";
  mmout << "Infty\t\t= " << Infinity (Ball) << "\n";
  mmout << "Fuzz\t\t= " << Fuzz (Ball) << "\n";
  mmout << "x\t\t= " << x << "\n";
  mmout << "y\t\t= " << y << "\n";
  mmout << "z\t\t= " << z << "\n";
  mmout << "sqrt (x)\t= " << sqrt (x) << "\n";
  mmout << "sqrt (y)\t= " << sqrt (y) << "\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 cos (z)\t= " << acos (z) << "\n";
  mmout << "arc tan (x)\t= " << atan (x) << "\n";
  mmout << "arc tan (y)\t= " << atan (y) << "\n";
  mmout << "arc tan (x)\t= " << atan (x) << "\n";
  mmout << "atan2 (x, y)\t= " << atan2 (x, y) << "\n";
  mmout << "atan2 (y, x)\t= " << atan2 (y, x) << "\n";
  mmout << "atan2 (0, -x)\t= " << atan2 (Ball (0), -x) << "\n";
  mmout << "min (x, y)\t= " << min (x, y) << "\n";
  mmout << "max (x, y)\t= " << max (x, y) << "\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_muller (prec);
  test_elementary_funcs (prec);
}

int
main () {
  test_constants ();
  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 << "1+(-infty)\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)^2\t= " <<
    (-Infinity (Ball)) * (-Infinity (Ball)) << "\n";
  mmout << "infty/infty\t= " <<
    Infinity (Ball) / Infinity (Ball) << "\n";
  mmout << "(-infty)/infty\t= " <<
    (-Infinity (Ball)) / Infinity (Ball) << "\n";
  mmout << "(-infty)/2\t= " << (-Infinity (Ball)) / Ball (2) << "\n";
  mmout << "|infty|\t\t= " << abs (Infinity (Ball)) << "\n";
  mmout << "|-infty|\t= " << abs (-Infinity (Ball)) << "\n";
  mmout << "center(infty)\t= " << center (Infinity (Ball)) << "\n";
  mmout << "center(-infty)\t= " << center (-Infinity (Ball)) << "\n";
  mmout << "radius(infty)\t= " << radius (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";
  mmout << "(-infty) << 3\t= " << ((-Infinity (Ball)) << 3) << "\n";
  mmout << "(-infty) >> 3\t= " << ((-Infinity (Ball)) >> 3) << "\n";
  return 0;
}