// Alogrithm 659, Collected Algorithm from ACM
// This is the C version of halton sequences 
// By Yaohang Li
// Dept. of Comp. Sci.
// Florida State Univ.
// 7/2000

#include <stdio.h>
#include <math.h>

int s;
double e;
double prime[]={2.0,3.0,5.0,7.0,11.0,13.0,17.0,19.0,23.0,29.0,31.0,37.0,41.0,43.0,47.0,53.0,59.0,61.0,67.0,71.0,73.0,79.0,83.0,89.0,97.0,101.0,103.0,107.0,109.0,113.0,127.0,131.0,137.0,139.0,149.0,151.0,157.0,163.0,167.0,173.0};

int inhalt(int dimen, int atmost, double tiny, double *quasi)
{
  double delta;
  int i;

// check dimen

  s=dimen;
  if (s<1||s>40)
    return(-1);

// compute and check tolerance

  e=0.9*(1.0/(atmost*prime[s-1])-10.0*tiny);
  delta=100*tiny*(double)(atmost+1)*log10((double)atmost);
  if (delta>=0.09*(e-10.0*tiny))
    return(-2);

// now compute first vector
  for (i=0;i<s;i++)
    {
      prime[i]=1/prime[i];
      quasi[i]=prime[i];
    }

  return(0);
}

int gohalt(double *quasi)
{
  double t,f,g,h;
  int i;

// generate quasi one compoment at a time using radix prime[k] for 
// component k

  for (i=0;i<s;i++)
    {
      t=prime[i];
      f=1.0-quasi[i];
      g=1.0;
      h=t;
      while ((f-h)<e)
// this checks whether q+h>1-e
        {
          g=h;
          h*=t;
        }
      quasi[i]=g+h-f;
    }

  return(0);
}

void main()
{
  double quasi[40];
  int dimen, n;
  int i,j;
  double tiny;

  dimen=2;
  n=100000;
  tiny=0;
  if (inhalt(dimen,n,tiny,quasi)<0)
    {
      printf("Error in dimension or tiny\n");
      return(-1);
    }

  for (i=0;i<n;i++)
    {
      gohalt(quasi);
      for (j=0;j<dimen;j++)
        printf("%f ",quasi[j]);
      printf("\n");
    }
}