halton.c File Reference

Detailed Description

Implementation of Halton and modified Halton sequences.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, R. Soliday, H. Shang, Y. Wang

Definition in file halton.c.

#include "mdb.h"

Go to the source code of this file.

Functions
int32_t	startHaltonSequence (int32_t *radix, double value)
	Initialize and start a new Halton sequence.
int32_t	restartHaltonSequence (long ID, double value)
	Restart an existing Halton sequence from a new initial value.
double	nextHaltonSequencePoint (long ID)
	Get the next point in a Halton sequence.
int32_t	power (int32_t a, int32_t b, int32_t m)
int32_t	primes ()
int	inhalt (int dimen, int atmost, double tiny, double *quasi)
int32_t	startModHaltonSequence (int32_t *radix, double tiny)
	Start a modified Halton sequence.
int32_t	restartModHaltonSequence (long ID, double tiny)
	Restart an existing modified Halton sequence.
int32_t	generateModHaltSequence (double *quasi, double *dq, double *wq, long ID)
double	nextModHaltonSequencePoint (long ID)
	Retrieve the next point from the modified Halton sequence.

Function Documentation

generateModHaltSequence()

int32_t generateModHaltSequence	(	double *	quasi,
		double *	dq,
		double *	wq,
		long	ID )

Definition at line 550 of file halton.c.

                                                                                {
  int32_t i = (int32_t)ID, j, k, ytemp[40], xtemp[40], ltemp, mtemp;
  double t, f, g, h;
  /* generate quasi one compoment at a time using radix prime[k] for 
     component k */
 
  t = iprime[i];
  f = 1.0 - quasi[i];
  g = 1.0;
  h = t;
  while ((f - h) < eError)
  /* this checks whether q+h>1-eError */
  {
    g = h;
    h *= t;
  }
  quasi[i] = g + h - f;
 
  k = 0;
  mtemp = nextPoint[i];
  ltemp = prime[i];
 
  while (mtemp != 0) {
    ytemp[k] = mtemp % ltemp;
    mtemp = mtemp / ltemp;
    k++;
  }
  /*generating Optimal primitive root  */
  for (j = 0; j < k; j++) {
    /* xtemp[j] = (ytemp[j]*power(primroots[i/10][i%10], nextn%ltemp, ltemp))%ltemp; */
    xtemp[j] = (warnockOpt[i] * power(primroots[i / 10][i % 10], ytemp[j], ltemp)) % ltemp;
    xtemp[j] -= ytemp[j];
  }
  dq[i] = 0;
  t = iprime[i];
  for (j = 0; j < k; j++) {
    dq[i] += xtemp[j] * t;
    t *= iprime[i];
  }
  dq[i] += quasi[i];
  /* generating Warnock Optimal sequences */
  for (j = 0; j < k; j++) {
    /*  if( i%2==0)xtemp[j]= (ytemp[j]*warnockOpt[i])%ltemp;
          else xtemp[j]= (ytemp[j]*warnockOpt[i]*warnockOpt[i])%ltemp;
      */
    xtemp[j] = (ytemp[j] * power(warnockOpt[i], i + 1, ltemp)) % ltemp;
    /* if(i>=3)
         xtemp[j]=((prime[i-1]-1)*power(warnockOpt[i],ytemp[j],ltemp))%ltemp;
     else  xtemp[j]= (power(warnockOpt[i],ytemp[j],ltemp))%ltemp;
      */
    xtemp[j] -= ytemp[j];
  }
 
  wq[i] = 0;
  t = iprime[i];
  for (j = 0; j < k; j++) {
    wq[i] += xtemp[j] * t;
    t *= iprime[i];
  }
  wq[i] += quasi[i];
 
  nextPoint[i]++;
 
  return (0);
}

inhalt()

int inhalt	(	int	dimen,
		int	atmost,
		double	tiny,
		double *	quasi )

Definition at line 468 of file halton.c.

                                                              {
  double delta;
  int i /*,m*/;
 
  /* check dimen */
  sDim = dimen;
  if (sDim < 1 || sDim > 1000)
    return (-1);
 
  /* compute and check tolerance*/
 
  eError = 0.9 * (1.0 / (atmost * prime[sDim - 1]) - 10.0 * tiny);
  delta = 100 * tiny * (double)(atmost + 1) * log10((double)atmost);
  if (delta >= 0.09 * (eError - 10.0 * tiny))
    return (-2);
 
  /* now compute first vector */
 
  /*m=1;*/
  for (i = 0; i < sDim; i++) {
    iprime[i] = 1.0 / prime[i];
    quasi[i] = iprime[i];
    /*m=i*prime[i];*/
    nextPoint[i] = 2;
  }
 
  /*printf("largest prime=%d, %d \n",prime[sDim-1], m);*/
 
  return 0;
}

nextHaltonSequencePoint()

double nextHaltonSequencePoint

(

long

ID

)

Get the next point in a Halton sequence.

Computes the next value in the specified Halton sequence.

Parameters

ID	The sequence ID.

Returns

The next point in the sequence, or -1 if the ID is invalid.

Definition at line 107 of file halton.c.

                                        {
  double r, f, value;
 
  ID -= 1;
 
  if (ID > sequencesInUse || ID < 0)
    return -1;
 
  f = 1 - lastPointValue[ID];
  r = 1. / R[ID];
  while (f <= r)
    r = r / R[ID];
  value = lastPointValue[ID] + (R[ID] + 1) * r - 1;
  lastPointValue[ID] = value;
  return value;
}

nextModHaltonSequencePoint()

double nextModHaltonSequencePoint

(

long

ID

)

Retrieve the next point from the modified Halton sequence.

Generates and returns the next value for a given sequence ID in the modified Halton sequence.

Parameters

ID	The sequence ID.

Returns

The next point in the modified Halton sequence.

Definition at line 624 of file halton.c.

                                           {
  static double *dq = NULL, *wq = NULL;
 
  ID -= 1;
  if (dq == NULL)
    dq = malloc(sizeof(*dq) * sDim);
  if (wq == NULL)
    wq = malloc(sizeof(*wq) * sDim);
  if (!modSequenceInUse) {
    fprintf(stderr, "ModHalton sequence not started.\n");
    exit(1);
  }
  if (ID < 0 || ID > sDim - 1) {
    fprintf(stderr, "Invalid ID (%ld) provided\n", ID);
    exit(1);
  }
  generateModHaltSequence(quasi, dq, wq, ID);
 
  return wq[ID];
}

power()

int32_t power	(	int32_t	a,
		int32_t	b,
		int32_t	m )

Definition at line 444 of file halton.c.

                                               {
  int32_t i, c = 1;
  for (i = 0; i < b; i++)
    c = (c * a) % m;
  return c;
}

primes()

int32_t primes

(

)

Definition at line 451 of file halton.c.

                 {
  int32_t i, j, a[MAX_D + 1];
  for (a[1] = 0, i = 2; i <= MAX_D; i++)
    a[i] = 1;
  for (i = 2; i <= MAX_D / 2; i++)
    for (j = 2; j <= MAX_D / i; j++)
      a[i * j] = 0;
  for (i = 1, j = 0; i <= MAX_D; i++) {
    if (a[i]) {
      prime[j] = i;
      iprime[j] = (double)i;
      j++;
    }
  }
  return j;
}

restartHaltonSequence()

int32_t restartHaltonSequence	(	long	ID,
		double	value )

Restart an existing Halton sequence from a new initial value.

Resets the specified Halton sequence so that subsequent calls will continue from the given initial value.

Parameters

ID	The sequence ID to restart.
value	The new starting value.

Returns

1 on success, or -1 if the ID is invalid.

Definition at line 88 of file halton.c.

                                                     {
  ID -= 1;
 
  if (ID > sequencesInUse || ID < 0)
    return -1;
 
  lastPointValue[ID] = value;
 
  return 1;
}

restartModHaltonSequence()

int32_t restartModHaltonSequence	(	long	ID,
		double	tiny )

Restart an existing modified Halton sequence.

Resets internal parameters so that the modified Halton sequence restarts.

Parameters

ID	The sequence ID.
tiny	A small tolerance value.

Returns

1 on success, or -1 on failure.

Definition at line 538 of file halton.c.

                                                       {
  int32_t dimen = 12, total_points = 100000;
 
  tiny = 0;
  if (inhalt(dimen, total_points, tiny, quasi) < 0) {
    fprintf(stderr, "Unable to start modHalton sequence.\n");
    return -1;
  }
 
  return 1;
}

startHaltonSequence()

int32_t startHaltonSequence	(	int32_t *	radix,
		double	value )

Initialize and start a new Halton sequence.

Initializes a new Halton sequence with the given radix and starting value. If the provided radix is not a prime or is non-positive, a suitable prime is chosen.

Parameters

radix	Pointer to an integer specifying the desired prime radix. If non-positive, a prime is chosen automatically.
value	The starting value for the sequence.

Returns

The sequence ID (1-based) on success, or 0 on failure.

Definition at line 38 of file halton.c.

                                                          {
  int32_t ID;
  if ((sequencesInUse == 0 &&
       (!(lastPointValue = malloc(sizeof(*lastPointValue))) ||
        !(R = malloc(sizeof(*R))))) ||
      (!(lastPointValue = realloc(lastPointValue, (sequencesInUse + 1) * sizeof(*lastPointValue))) ||
       !(R = realloc(R, (sequencesInUse + 1) * sizeof(*R)))))
    return 0;
  if (*radix > 0) {
    if (is_prime(*radix) != 1)
      return 0;
    R[sequencesInUse] = *radix;
    ID = sequencesInUse;
  } else {
    /* generate a new, unique prime number for use as the radix */
    long i, passed;
    if ((ID = sequencesInUse) < N_SEQ_PREDEFINED)
      /* try one of the favored values */
      *radix = Rvalues[ID];
    else
      *radix = 2;
    passed = 0;
    while (!passed) {
      passed = 1;
      for (i = 0; i < sequencesInUse; i++) {
        if (R[i] == *radix) {
          passed = 0;
          (*radix)++;
          while (is_prime(*radix) != 1)
            (*radix)++;
        }
      }
    }
    R[ID] = *radix;
  }
  lastPointValue[ID] = value;
  sequencesInUse++;
  return ID + 1;
}

startModHaltonSequence()

int32_t startModHaltonSequence	(	int32_t *	radix,
		double	tiny )

Start a modified Halton sequence.

Initializes a modified Halton sequence with predefined prime bases and internal parameters.

Parameters

radix	Pointer to an integer to store the chosen prime radix.
tiny	A small tolerance value.

Returns

The sequence ID (1-based) on success, or -1 on failure.

Definition at line 508 of file halton.c.

                                                            {
  int32_t modID, dimen = 12, total_points = 100000;
  tiny = 0;
  /*  check dimen*/
  if (!modSequenceInUse) {
    /*generate primes
      primes(); */
    /* commented out primes generation and put it as global constants to skip generate it every time when start the sequence*/
    if (!quasi)
      quasi = malloc(sizeof(*quasi) * sDim);
    if (inhalt(dimen, total_points, tiny, quasi) < 0) {
      fprintf(stderr, "Unable to start modHalton sequence.\n");
      return -1;
    }
  }
  modID = modSequenceInUse;
  *radix = prime[modID];
  modSequenceInUse++;
  return modID + 1;
}