#include <cassert>
#include <iostream>
#include <vector>
#include <climits>
#include <cmath>
using namespace std;

#include <aclock.h>
#include <bernsteintest.h>
#include <binoppropertiestest.h>
#include <commandline.h>
#include <d2homogeneoustest.h>
#include <disktest.h>
#include <fixedpointtest.h>
#include <functest.h>
#include <gausselimtest.h>
#include <halfspaceD2test.h>
#include <halfspaceD3test.h>
#include <histogramtest.h>
#include <linechoppedtest.h>
#include <linetest.h>
#include <mathlib.h>
#include <mathlibtest.h>
#include <planetest.h>
#include <point.h>
#include <pointtest.h>
#include <powerseriestest.h>
#include <print.h>
#include <randomtest.h>
#include <simplexD1linkedtest.h>
#include <vectest.h>
#include <zero.h>

template<>
double zero<double>::val = 1E-15;

typedef point2<double> pt2;


/*
void test01()
{
  RR::SetPrecision(200);
  RR::SetOutputPrecision(200);

  RR x;
  cin >> x; 

  cout << "x=" << x << endl;
}

void test02()
{
  RR x;
  RR y;

  x = 0.2;

  arcsin<RR> w;

  RR::SetPrecision(200);
  RR::SetOutputPrecision(30);

  unsigned int imax=20;
  for (unsigned int i=2; i<imax; ++i)
    cout << w(y,x,i) << endl;
}
*/

/*
template<class T>
void print(vector<T> const & v)
{
  for (unsigned int i=0; i<v.size(); ++i)
    cout << v[i] << " ";
  cout << endl;
}
*/


void test03()
{

  vector<double> a1(3);
  a1[0] = 1;
  a1[1] = 0;
  a1[2] = 0;

  vector<double> a2(3), a3(3);

  a2[0] = 0;
  a2[1] = 1;
  a2[2] = 0;

  //surfaceplane< vector<double> >()(a2,a3,a1);

  print(a1);
  print(a2);
  cout << endl;
  //crossprod< vector<double> >()(a3,a1,a2);
  crossproduct::eval(a3,a1,a2);

  print(a3);

  cout << "And the reverse process, from a normal" << endl;
  cout << "generate a surface plane." << endl;

  vector<double> u(3), v(3), t(3);
  surfaceplane< vector<double> >()(u,v,t,a3); 
  cout << "u=";
  print(u);
  cout << "v=";
  print(v);
  cout << "w=";
  print(a3);

}

class f1
{
public:

  void operator () (double& y, doublec & x)
  {
    y = (x-3.0)*(x+2)*x;
  }
};

void test04()
{

  convexity<double> c;
  bool res;

  f1 f;

  unsigned int n;
  double a;
  double b;

  cout << "Testing convextity test." << endl;
  cout << "Enter n a b where n is the number of points" << endl;
  cout << "  sampled on (a,b). n=0 terminates." << endl;

  for ( ; ; )
  {

    cin >> n;
    if (n==0)
      return;

    cin >> a;
    cin >> b;

    c(res,f,n,a,b);
    cout << "res=" << res << endl;
  }

}

/*
void test05()
{

  orthoproj
  < 
    double, 
    point2<double>, 
    innerprod2<double>,
    2
  > p;

  p.v[0] = point2<double>(1.0,2.0); 
  p.v[1] = point2<double>(-2.0,1.0);

  point2<double> c(6.0,4.0);

  p(c);

  cout << "Verify 2D projection." << endl;

  cout << "c=" << c << endl;
  cout << p.v[0] * p.a[0] << endl;
  cout << p.v[1] * p.a[1] << endl;
  cout << p.v[0] * p.a[0] + p.v[1] * p.a[1] << endl;
  cout << endl;

  cout << "Verify 3D projection." << endl;

  point3<double> e(4.0, 2.6, -3.0);
  orthoproj<double,point3<double>,innerprod3<double>,3> w;
  w.v[0] = point3<double>(2.5,3.0,-1.0);
  point3<double> t;
  surfaceplane<point3<double> >() (w.v[1],w.v[2],t,w.v[0]);
  w(e);

  cout << "e=" << e << endl;
  cout << w.v[0] * w.a[0] << endl;
  cout << w.v[1] * w.a[1] << endl;
  cout << w.v[2] * w.a[2] << endl;
  cout << w.v[0] * w.a[0] + w.v[1] * w.a[1] + w.v[2] * w.a[2] << endl;
  cout << endl;

}
*/

void test06()
{

  cout << "Triangle Area Calculation" << endl;
  
  cout << "x0 y0 x1 y1 x2 y2" << endl;
  double x0,y0,x1,y1,x2,y2;
  cin >> x0;
  cin >> y0;
  cin >> x1;
  cin >> y1;
  cin >> x2;
  cin >> y2;
  double f;
  trianglearea(f,x0,y0,x1,y1,x2,y2);
  cout << "area=" << f << endl;

  
  cout << "area=" << 0.5*determinant::calcD2(x1-x0,x2-x0,y1-y0,y2-y0) << endl;

}


void test07()
{
  pt2 x;

  d2matsolve(x, pt2(3.0,2.0), pt2(1.0,1.0), pt2(7.0,3.0) );
  cout << "x=" << x << endl;

}







void test10()
{
  cout << "Testing Area of Triangle Formula in 3D" << endl;

  point3<double> p0(0.0,0.0,0.0);
  point3<double> p1(3.0,0.0,0.0);
  point3<double> p2(3.0,4.0,0.0);

  double a1;
  trianglearea(a1,p0.x,p0.y,p1.x,p1.y,p2.x,p2.y);
  cout << SHOW(a1) << endl;

  double a2;

  trianglearea(a2,p0,p1,p2);
  cout << SHOW(a2) << endl;


}

void test11()
{
  cout << "Testing solvequadratic" << endl;
  cout << "Solve (x-2)*(x+3)=0" << endl;
  // x^2 + x - 6 = 0
  pt2 x;
  solvequadratic(x.x,x.y,1.0,1.0,-6.0);
  cout << "x=" << x << endl;
}

void test12()
{
  cout << "Test intersecting lines equation" << endl;
 
  pt2 a(1.0,1.0);
  pt2 b(5.0,1.0);

  pt2 c(2.0,-1.0);
  pt2 d(2.0, 7.0);

  double t0;
  double t1;

  bool res;

  res = lineintersection(t0,t1,a,b,c,d);
  cout << SHOW(res) << endl;
  cout << "(" << a << "," << b << "): " << a+(b-a)*t0 << endl;
  cout << "(" << c << "," << d << "): " << c+(d-c)*t1 << endl;

}

void test13()
{
  cout << "Test rotating a 2D point about the origin" << endl;

  point2<double> x(1.0,0.0);

  cout << SHOW(x) << endl;

  doublec pi=3.141592654;
  transrotate2D r(pi*0.5);

  r.eval(x);
  cout << SHOW(x) << endl;
}

void test014()
{
  cout << "Testing 3D matrix multiplication." << endl;

  doublec m1[] = 
  { 
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0,
    7.0, 8.0, 9.0
  };

  point3<double> const x0(1.0,1.0,1.0);
  point3<double> x(x0);

  cout << SHOW(x) << endl;

  cout << "multiply with double * matrix" << endl;

  point3<double> y(x);
  matrixmult(x,m1,y);
  cout << SHOW(x) << endl;

  point3<double> r0(m1[0],m1[1],m1[2]);
  point3<double> r1(m1[3],m1[4],m1[5]);
  point3<double> r2(m1[6],m1[7],m1[8]);

  cout << "multiply with rows as 3D points for the matrix" << endl;
  matrixmult(x,r0,r1,r2,x0);
  cout << SHOW(x) << endl;
  
}



void test016()
{
  cout << "Testing line segment intersection" << endl;

  point2<double> p1(0.0,0.0);
  point2<double> p2(1.0,0.0);
  point2<double> q1(0.5,0.0);
  point2<double> q2(0.5,1.0);

  cout << SHOW(p1) << " " << SHOW(p2) << endl;
  cout << SHOW(q1) << " " << SHOW(q2) << endl;

  cout << "Without division algorithm: " << endl;

  bool const res = lineSegmentIntersection(p1,p2,q1,q2,1E-5);
  cout << SHOW(res) << endl;
  cout << endl;

  cout << "With One division" << endl;
  double tp,tq;
  bool const res2 = lineIntersection(tp,tq,p1,p2,q1,q2,1E-5);
  cout << SHOW(res2) << endl;
  cout << SHOW(tp) << "  (" << p1 + (p2-p1)*tp << ")" << endl;
  cout << SHOW(tq) << "  (" << q1 + (q2-q1)*tq << ")" << endl;
  cout << endl;
}

void test017data( vector< point2<double> > & v, uintc n)
{
  doublec w = 1.0 / LONG_MAX;

  for (uint i=0; i<n*4; ++i)
    v.push_back( point2<double>(w*rand(),w*rand()) );
}

void test017()
{
  cout << "Testing line segment intersection algorithms" << endl;
  cout << "Algorithm A has no divisions, Algorithm B has one" <<endl;
  cout << "Generating test data : random points in a unit cube" << endl;

  uintc n = 100000;
  cout << n << " lines are being intersected" << endl;
  
  vector< point2<double> > v;
  test017data(v,n);

  vector<bool> b0(n);

  uint k;

  aclock clk;

  clk.measure();
  for (uint i=0; i<n; ++i)
  {
    k = i*4;
    b0[i] = lineSegmentIntersection(v[k],v[k+1],v[k+2],v[k+3],1E-5);
  } 
  clk.measure();
  cout << SHOW(clk.diff_ms()) << endl;

  cout << "Algorithm B has one division" << endl;
  vector<bool> b1(n);
  double t1,t2;

  clk.measure();
  for (uint i=0; i<n; ++i)
  {
    k = i*4;
    b1[i] = lineSegmentIntersection(t1,t2,v[k],v[k+1],v[k+2],v[k+3],1E-5);
  } 
  clk.measure();
  cout << SHOW(clk.diff_ms()) << endl;

  uint count=0;
  for (uint i=0; i<n; ++i)
  {
    if (b0[i]!=b1[i])
      ++count;
  }

  cout << "Checking the two algorithms for differences in intersection results" << endl;
  cout << "  differences = " << count << endl;

  cout << "My results : for optimized compiler option there was no speed difference" << endl;
  cout << "  between algorithms A and B. ie the division was as fast as the integer version." << endl;
  cout << "  Since the division algorithm B also calculates the intersection point it is superior." << endl;
  cout << "  This result must be because of the processor." << endl;

}

void test018()
{
  cout << "Comparing timings for line intersection with algorithm B and C" << endl;

  uintc n = 100000;
  cout << n << " lines are being intersected" << endl;

  vector< point2<double> > v;
  test017data(v,n);

  vector<bool> b0(n);

  uint k;

  aclock clk;
  double t1,t2;

  cout << "Algorithm C uses a 2D matrix solver" << endl;

  clk.measure();
  for (uint i=0; i<n; ++i)
  {
    k = i*4;
    b0[i] = lineintersection(t1,t2,v[k],v[k+1],v[k+2],v[k+3]);
  } 
  clk.measure();
  cout << SHOW(clk.diff_ms()) << endl;

  cout << "Algorithm B has one division" << endl;
  vector<bool> b1(n);

  clk.measure();
  for (uint i=0; i<n; ++i)
  {
    k = i*4;
    b1[i] = lineIntersection(t1,t2,v[k],v[k+1],v[k+2],v[k+3],1E-5);
  } 
  clk.measure();
  cout << SHOW(clk.diff_ms()) << endl;

  uint count=0;
  for (uint i=0; i<n; ++i)
  {
    if (b0[i]!=b1[i])
      ++count;
  }

  cout << "Checking the two algorithms for differences in intersection results" << endl;
  cout << "  differences = " << count << endl;
}



void test021()
{
  cout << "Testing the volume of a tetraheron calc." << endl;

  doublec a2[] = 
  {
    0.0, 0.0, 0.0,
    1.0, 0.0, 0.0,
    1.0, 1.0, 0.0,
    0.0, 1.0, 0.0,
    0.0, 0.0, 1.0,
    1.0, 0.0, 1.0,
    1.0, 1.0, 1.0,
    0.0, 1.0, 1.0
  };

  uint n=8;

  typedef point3<double> pt3;
  vector< pt3 > v2;
  for (uint i=0; i<n; ++i)
    v2.push_back( pt3(a2[i*3],a2[i*3+1],a2[i*3+2]) ); 


  double volume;
  tetrahedronvolume(volume,v2[0],v2[1],v2[5],v2[6]);
  cout << SHOW(volume) << endl;

  double s=0.0;
  s += abs(volume);
  tetrahedronvolume(volume,v2[0],v2[6],v2[2],v2[1]);
  cout << SHOW(volume) << endl;
  s += abs(volume);
  tetrahedronvolume(volume,v2[2],v2[6],v2[0],v2[3]);
  cout << SHOW(volume) << endl;
  s += abs(volume);
  cout << SHOW(s) << endl;

  cout << "The sum should be 0.5 or half the volume of a unit cube" << endl;

}

void test022()
{
  cout << "Combination" << endl;
  cout << "Enter n k" << endl;
  uint n,k;
  cin >> n >> k;
  double res;
  mathcombination(res,n,k);
  cout << SHOW(res) << endl;
}







int main(int argc, char** argv)
{
  commandline cmd(argc,argv);
  uint prog(0);
  cmd.mapvar(prog,"prog");

  switch (prog)
  {
    case 0:
      cout << "Test programs for mathlib workspace" << endl;
      cout << "$./main prog=1..2      - Gaussian elimination tests." << endl;
      cout << "$./main prog=40..41    - Power series tests" << endl;

      cout << "$./main prog=17        - Testing line segment intersection algorithms" << endl;
      cout << "$./main prog=16        - Testing two line segment intersection" << endl;
      cout << "                         (without division)" << endl;
      cout << "$./main prog=15        - Testing d2homogeneous class." << endl;
      cout << "$./main prog=14        - Testing 3D matrix multiplication." << endl;
      cout << "$./main prog=21        - Testing tetrahedron volume calc." << endl;
      cout << "$./main prog=22..23    - Testing vec" << endl;
      cout << "$./main prog=25..30    - Testing func.h" << endl;
      cout << "$./main prog=31..35    - Testing random.h" << endl;
      cout << "$./main prog=55        - Testing intervalintersection." << endl;
      cout << "$./main prog=70..74    - Testing permfunc and binopproperties classes." << endl;

      cout << "$./main prog=131       - Line test." << endl;
      cout << "$./main prog=140       - Histogram test." << endl;
      cout << "$./main prog=145       - " << simplexD1linkedtest::doc[1] << endl;
      cout << "$./main prog=180-181    - Testing funcstate.h" << endl;
      break;

    case 1: gausselimtest::test01(); break;
    case 2: gausselimtest::test02(); break;

    case 18: test018(); break;
    case 17: test017(); break;
    case 16: test016(); break;
    case 15: d2homogeneoustest01(); break;
    case 14: test014(); break;

    case 21: test021(); break;
    case 22: vectest::test01(); break;
    case 23: vectest::test02(); break;

    case 25: functest::test01(); break;
    case 26: functest::test02(); break;
    case 27: functest::test03(); break;
    case 28: functest::test04(); break;
    case 29: functest::test05(); break;
    case 30: functest::test06(); break;
    case 31: return randomtest::test01(); 
    case 32: return randomtest::test02(); 
    case 33: return randomtest::test03();
    case 34: return randomtest::test04();
    case 35: return randomtest::test05();
    case 36: randomtest::test06(); break;
    case 37: randomtest::test07(); break;

    case 40: powerseriestest::test01(); break;
    case 41: powerseriestest::test02(); break;

    case 45: test022(); break;
    case 46: bernsteintest::test01(); break;

    case 50: mathlibtest::test01(argc,argv); break;
    case 51: mathlibtest::test02(argc,argv); break;
    case 52: mathlibtest::test03(); break;

    case 55: mathlibtest::test04(argc,argv); break;
    case 56: mathlibtest::test05(); break;
    case 57: mathlibtest::test06(argc,argv); break;
    case 58: mathlibtest::test07(argc,argv); break;

    case 60: planetest::test01(argc,argv); break;
    case 61: planetest::test02(argc,argv); break;
    case 62: planetest::test03(argc,argv); break;

    case 65: pointtest::test01(); break;
    case 66: pointtest::test02(); break;
    case 67: pointtest::test03(); break;

    case 70: binoppropertiestest::test01(); break;
    case 71: binoppropertiestest::test02(); break;
    case 72: binoppropertiestest::test03(); break;
    case 73: binoppropertiestest::test04(); break;
    case 74: binoppropertiestest::test05(); break;

    case 80: halfspaceD2test::test01(); break;
    case 81: halfspaceD2test::test02(); break;
    case 82: return halfspaceD2test::unittest01(); 
    case 85: halfspaceD3test::test01(); break;
    case 86: halfspaceD3test::test02(); break;
    case 87: halfspaceD3test::test03(argc,argv); break;

    case 90: disktest::test00(); break;

    case 100: mathlibtest::test010(); break;
    case 101: mathlibtest::test011(); break;
    case 105: mathlibtest::test015(); break;

    case 110: fixedpointtest::test02(); break;
    case 111: fixedpointtest::test03(); break;

    case 131: linetest::test01(); break;
    case 132: linetest::test02(); break;
    case 133: return linetest::unittest01();

    case 135: linechoppedtest::unittest01(); break;
    case 136: linechoppedtest::unittest02(); break;

    case 140: histogramtest::test01(argc,argv); break;
    
    case 145: simplexD1linkedtest::test01(); break;

    case 180: functest::test07(); break;
    case 181: functest::test08(); break;

    default: cout << "error:  No case handled." << endl; return 1;
  }

  return 0;
}