/*
  Hammersley generation of points in a D-parallelepid

	Generate N points into a hyperparallelepid (D dimensions) more
  "regularly" than the random method. 
	Two options:
	- pure classical Hammersley method
	- a slightly improved one

	Last update: 2008-11-22

 Copyright (C) Maurice Clerc 2008 <Maurice.Clerc@WriteMe.com>
  
 * Free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * See the GNU General Public License for more details:
	http://www.gnu.org/licenses/
 with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#define DMax 5 // Max dimension of the search space
#define Nmax 101  // Max number of points

#include <stdio.h>
#include <stdlib.h>

double alea(double a, double b);
int aleaInteger(int min, int max);

int main()
{
  double a;
  int d;
  int D; // Dimension of the search space
  int deb,dd;
  int k, kp;
  int N;
  int option;
  int pp;
  double Phi[Nmax][DMax];
  double xMin[DMax], xMax[DMax];

  int P[DMax]={2,3,5,7,11}; // A sequence of prime numbers
														// You may modify it
														// Note: a possible variant is to 
														// randomly permute it before the generation
  double z;

  FILE *f_hamm_hal;
  f_hamm_hal=fopen("f_hamm_hal.txt","w");

  D=2; // Dimension
  N=100; // Number of points
  option=1; // 0 => pure Hammersley (deterministic)
			// 1 => improved (no point on (0, ...0)), and 
			// 			the first dimension is chosen at random

	// Define the D-parallelepid (not necessarily a hypercube)
  for (d=0;d<D;d++)
  {
		xMin[d]=0;
		xMax[d]=100;
  }

  switch (option)
  {
		default: // Pure Hammersley
		for (k=0;k<N;k++) // For each point
		{
			Phi[k][0]=xMin[0]+(xMax[0]-xMin[0])*(double)k/N;
			for(d=1;d<D;d++) // for each dimension
			{				// compute the position
			pp=P[d-1];
			kp=k;
			z=0;
			while (kp>0)
			{
				a=kp % P[d-1];     // kp modulo P[d-1]
				z=z+(double)a/pp;
				kp=kp/P[d-1]; // Integer division
				pp=pp*P[d-1];
			}
			Phi[k][d]=xMin[d]+z*(xMax[d]-xMin[d]);
			}
		}
		break;
		case 1: // Improved
		deb=aleaInteger(0,D-1);

		for (k=0;k<N;k++) // For each point
		{
			Phi[k][deb]=xMin[deb]+(xMax[0]-xMin[0])*(double)k/(N-1);
	  
			for(d=1;d<D;d++) // for each dimension
			{				// compute the position
			dd=deb+d; if(dd>=D) dd=0;
			pp=P[dd];
			kp=k+1;
			z=0;
			while (kp>0)
			{
				a=kp % P[dd];     // kp modulo P[d]
				z=z+(double)a/pp;
				kp=kp/P[dd]; // Integer division
				pp=pp*P[dd];
			}
			Phi[k][dd]=xMin[dd]+z*(xMax[dd]-xMin[dd]);
			}
		}
		break;
  }
	
  // Display the positions
  for (k=0;k<N;k++)
  {
		printf("\n");
		for(d=0;d<D;d++)
		{
			printf("%3.3f ",Phi[k][d]);
		}
  }

  // Save the positions
  for (k=0;k<N;k++)
  {
		for(d=0;d<D;d++)
		{
			fprintf(f_hamm_hal,"%3.3f ",Phi[k][d]);
		}
		fprintf(f_hamm_hal,"\n");
  }

	return (0);
}
//================================================
double alea(double a, double b) 
/* Random real  number between a and b */
{
	double 	r;

	// Normally, RAND_MAX = 32767 = 2^15-1
	r=(double)rand()/RAND_MAX;
	// Not very good. May be replaced by KISS
	//r=(double)rand_kiss()/RAND_MAX_KISS;

	r= a+r*(b-a);
	return r;
}
//================================================
int aleaInteger(int min, int max)
/* Random integer number between min and max */
{
	double 	r;

	if(min>=max) return min;

	r=alea(min,max);
	return (int)(r+0.5);
}