From: Martin Str|mberg <ams AT ludd DOT luth DOT se>
Message-Id: <199811221725.SAA04744@father.ludd.luth.se>
Subject: More rand()
To: djgpp-workers AT delorie DOT com (DJGPP-WORKERS)
Date: Sun, 22 Nov 1998 18:25:46 +0100 (MET)
X-Mailer: ELM [version 2.4ME+ PL15 (25)]
MIME-Version: 1.0
Content-Type: text/plain; charset=iso-8859-1
Content-Transfer-Encoding: 8bit
Reply-To: djgpp-workers AT delorie DOT com

Here's a first try. I needed an erand48(), so I threw in the whole bunch
while I was going. But the only one I really have tried erand48().

Any suggestions?


Right,

							MartinS

----- begin rand48.h -----
/* Copyright (C) 1998 Martin Strömberg. */

/* This function returns a double uniformly distributed in [0, 1). It
   uses an internal state common to all the three functions drand48,
   lrand48 and mrand48. */
extern double
drand48(void);

/* This function returns a double uniformly distributed in [0, 1). It
   uses the state supplied. */
extern double 
erand48(
	unsigned short state[3]
	);

/* This function returns an unsigned interger uniformly distributed in
   [0, 2^31]. It uses an internal state common to all the three
   functions drand48, lrand48 and mrand48. */
extern unsigned long
lrand48(void);

/* This function returns an unsigned interger uniformly distributed in
   [0, 2^31]. It uses the state supplied. */
extern unsigned long
nrand48(
	unsigned short state[3]
	);

/* This function returns a signed interger uniformly distributed in
   [0, 2^31]. It uses an internal state common to all the three
   functions drand48, lrand48 and mrand48. */
extern long
mrand48(void);

/* This function returns a signed interger uniformly distributed in
   [0, 2^31].  It uses the state supplied. */
extern long
jrand48(
	unsigned short state[3]
	);

/* This function sets the high order 32-bits of the internal state
   used by the funtions drand48, lrand48 and mrand48. The low order
   16-bits are set to the arbitrary value 0x330e. Addidtionaly this
   function resets the multiplier and additiver (possbly changed by a
   call to lcong48) to the default. */
extern void
srand48(
	long seedval
	);

/* This function sets the internal state to the supplied
   state. Addidtionaly this function resets the multiplier and
   additiver (possbly changed by a call to lcong48) to the default. */
extern unsigned short *
seed48(
       unsigned short seed16v[3]
       );

/* This function sets the internal state to the state supllied in
   param[0], param[1] and param[2]. It sets the multiplier to
   (param[5]*2^32 + param[4]*2^16 + param[3]) and the additiver to
   param[6]. */
extern void
lcong48(
	unsigned short param[7]
	);

----- end rand48.h -----

----- begin rand48.c -----
/* Copyright (C) 1998 Martin Strömberg. */

#include "rand48.h"

#define RAND48_MULT0 (0xe66d)
#define RAND48_MULT1 (0xdeec)
#define RAND48_MULT2 (0x0005)
#define RAND48_ADD (0x000b)

static unsigned short state[3] = {1, 0, 0};
static unsigned short multiplier0 = RAND48_MULT0;
static unsigned short multiplier1 = RAND48_MULT1;
static unsigned short multiplier2 = RAND48_MULT2;
static unsigned short additiver = RAND48_ADD;

static void next( 
		 unsigned short state[]
		 )
{
  unsigned short new_state[3];
  unsigned long tmp;

  tmp = state[0] * multiplier0 + additiver;
  new_state[0] = (unsigned short)(tmp  & 0xffff);

  tmp = (tmp >> 8*sizeof(unsigned short)) 
      + state[0] * multiplier1 
      + state[1] * multiplier0; 
  new_state[1] = (unsigned short)(tmp & 0xffff);

  tmp = (tmp >> 8*sizeof(unsigned short))
      + state[0] * multiplier2
      + state[1] * multiplier1
      + state[2] * multiplier0;
  new_state[2] = (unsigned short)(tmp & 0xffff);

  memcpy(state, new_state, 3*sizeof(unsigned short));
}

double drand48(void)
{
  return(erand48(state));
}

double erand48(
	       unsigned short state[3]
	       )
{
  int i; /* Counter. */
  double pot = 0.5; /* A potential of 0.5. */
  double result = 0.0;

  next(state);

  for(i = 0; i < 8*sizeof(unsigned short); i++)
  {
    if( (state[2] << i) & 0x8000 )
    {
      result += pot;
    }
    pot /= 2.0;
    if( (state[1] << i) & 0x8000 )
    {
      result += pot;
    }
    pot /= 2.0;
    if( (state[0] << i) & 0x8000 )
    {
      result += pot;
    }
    pot /= 2.0;
  }

  return(result);

}


unsigned long lrand48(void)
{
  return(nrand48(state));
}


unsigned long nrand48(
		      unsigned short state[3]
		      )
{

  next(state);
  return( ((unsigned long)state[2]) * 0x10000 + ((unsigned long)state[1]) );

}


long mrand48(void)
{
  return(jrand48(state));
}


long jrand48(
	     unsigned short state[3]
	     )
{

  next(state);
  if( (state[2] & 0x8000) )
  {
    return( -1.0 * ((long)(state[2] & 0x7fff)) * 0x10000 + ((unsigned long)state[1]) );
  }
  else
  {
    return( ((long)(state[2] & 0x7fff)) * 0x10000 + ((unsigned long)state[1]) );
  }

}

void srand48(
	     long seedval
	     )
{

  /* Restore default multipliers and additiver. */
  multiplier0 = RAND48_MULT0;
  multiplier1 = RAND48_MULT1;
  multiplier2 = RAND48_MULT2;
  additiver = RAND48_ADD;

  /* Setup the new state. */
  state[0] = 0x330e;
  state[1] = (seedval & 0xffff);
  state[2] = ( (seedval >> 16) & 0xffff);

}

unsigned short *seed48(
		       unsigned short state_seed[3]
		       )
{
  static unsigned short old_state[3];

  /* Restore default multipliers and additiver. */
  multiplier0 = RAND48_MULT0;
  multiplier1 = RAND48_MULT1;
  multiplier2 = RAND48_MULT2;
  additiver = RAND48_ADD;

  /* Remember old state. */
  memcpy(old_state, state, 3*sizeof(unsigned short));

  /* Setup the new state. */
  memcpy(state, state_seed, 3*sizeof(unsigned short));

  return(old_state);
}

void lcong48(
	     unsigned short param[7]
	     )
{
  
  /* Set the state. */
  state[0] = param[0];
  state[1] = param[1];
  state[2] = param[2];

  /* Set the multipilers. */
  multiplier0 = param[3];
  multiplier1 = param[4];
  multiplier2 = param[5];

  /* Set the additiver. */
  additiver = param[6];

}

----- end rabd48.c -----