#define MATHEPOOL

#include <list>
#include <fstream>
#include <string>
#include <CGAL/Cartesian.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Point_2.h>
#include <CGAL/Random.h>

using namespace CGAL;
using namespace std;

typedef Gmpq number;
typedef Cartesian<number> K;
typedef Point_2<K> Point;

void construct_grid (list<Point> & pts, int n)
{
   for (int i = 1; i <= n; ++i)
   {
      number x (2*i-1,2*n);
      for (int j = 1; j <= n; ++j)
      {
	 number y (2*j-1,2*n);
	 pts.push_back (Point(x,y));
      }
   }
}

void construct_random (list<Point> & pts, int n)
{
   Random randgen;
   for (int i = 0; i < n; ++i)
      pts.push_back (Point(randgen.get_double(),randgen.get_double()));
}

void construct_gridrand (list<Point> & pts, int n)
{
   Random randgen;
   for (int i = 0; i < n; ++i)
   {
      number x (i,n);
      for (int j = 0; j < n; ++j)
      {
	 number y (j,n);
	 pts.push_back (Point(x+randgen.get_double()/n,y+randgen.get_double()/n));
      }
   }
}

void construct_lattice (list<Point> & pts, int n, double a)
{
   for (int i = 0; i < n; ++i)
   {
      number x (i,n);
      number y (a*i-floor(a*i));
      pts.push_back (Point(x,y));
   }
}

void construct_corput (list<Point> & pts, int n)
{
   for (int i = 0; i < n; ++i)
   {
      number x (i,n);
      double ri = 0;
      double div = 1;
      div /= 2;
      for (int ic = i; ic > 0; ic >>= 1, div /= 2)
	 if (ic%2 == 1)
	    ri += div;
      number y (ri);
      pts.push_back (Point(x,y));
   }
}

int print_usage ()
{
   cout << "usage: pointlist <filename> <type> options\n";
   cout << "    where <type> is { grid | random | gridrand | lattice | corput }\n";
   cout << "    options for <type> =\n";
   cout << "      = grid :     n > 0  -- square length (yields n*n points)\n";
   cout << "      = random :   n > 0  -- total number of points\n";
   cout << "      = gridrand : n > 0  -- same as for grid\n";
   cout << "      = lattice :  n > 0  -- total number of points\n";
   cout << "                   optional: a > 0  -- irrational number ...good luck :)\n";
   cout << "      = corput :   n > 0  -- total number of points\n";
   return -1;
}

int main (int argc, char *args[])
{
   enum { undef, grid, random, gridrand, lattice, corput } ptstype;

   if (argc < 4)
      return print_usage();
   
   // parse filename
   if (strlen(args[1])==0)
      return print_usage();

   // parse type
   if (strcmp(args[2],"grid")==0)
      ptstype = grid;
   else if (strcmp(args[2],"random")==0)
      ptstype = random;
   else if (strcmp(args[2],"gridrand")==0)
      ptstype = gridrand;
   else if (strcmp(args[2],"lattice")==0)
      ptstype = lattice;
   else if (strcmp(args[2],"corput")==0)
      ptstype = corput;
   else
      return print_usage();

   // parse size
   int n = strtol (args[3], NULL, 10);
   if (n < 1)
      return print_usage();

   // parse optional `a' parameter for lattice
#ifdef MATHEPOOL
   double a = (std::sqrt((double)5)+1)/2;
#else
   double a = (sqrt((long double)5)+1)/2;
#endif
   if ((ptstype == lattice) && (argc > 4))
      a = strtod (args[4], NULL);

   // construct points
   list<Point> ptlist;
   switch (ptstype)
   {
      case grid:
	 construct_grid (ptlist, n);
	 break;
      case random:
	 construct_random (ptlist, n);
	 break;
      case gridrand:
	 construct_gridrand (ptlist, n);
	 break;
      case lattice:
	 construct_lattice (ptlist, n, a);
	 break;
      case corput:
	 construct_corput (ptlist, n);
	 break;
      default:
	 cout << "whoops, this shouldn't happen...\n";
	 break;
   }

   // write
   if (ptlist.size() > 0)
   {
      ofstream outfile (args[1]);
      if (outfile)
	 for (list<Point>::const_iterator p = ptlist.begin(); p != ptlist.end(); ++p)
	    outfile << *p << "\n";
   }
       
   return 0;
}