/**************************************************************************
 *{@C
 *      Copyright:      1988-2023 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         Utilities
 *      Filename:       UT_Random.c
 *
 *      Author:         Paul Obermeier
 *
 *      Description:    Machine independent pseudo-random number generators.
 *                      The functions defined below are similar to the
 *                      "drand48" and "lrand48" functions available with
 *                      System V Unix. They are re-implemented here so that
 *                      we have identical random number generators available
 *                      on all types of computers.
 *
 *      Additional documentation:
 *                      None.
 *
 *      Exported functions:
 *                      UT_RandomNew
 *                      UT_RandomDelete
 *                      UT_RandomGetInt32
 *                      UT_RandomGetFloat32
 *
 **************************************************************************/

#include <stdio.h>

#include "UT_Compat.h"
#include "UT_Memory.h"

#include "UT_Random.h"
#include "UT_RandomPrivate.h"

static const UInt16
    a[3] = {0xe66d, 0xdeec, 0x0005},
    c = 0xb;

/* A random unsigned integer number is calculated.
   The numbers are uniformly distributed over the interval [ 0, (1<<32) - 1 ).
*/

static UInt32 calcRandom (UT_RandomGen rGen)
{
    UInt32 y[3];

    /* Compute (r * a + c) % (1 << 48), using 48-bit integer arithmetic. */
    y[2] =                   (rGen->r0 * a[2]);
    y[2] = (y[2] & 0xffff) + (rGen->r1 * a[1]);
    y[2] = (y[2] & 0xffff) + (rGen->r2 * a[0]);
    y[1] =                   (rGen->r0 * a[1]);
    y[2] = (y[2] & 0xffff) + (y[1] >> 16);
    y[1] = (y[1] & 0xffff) + (rGen->r1 * a[0]);
    y[2] = (y[2] & 0xffff) + (y[1] >> 16);
    y[0] =                   (rGen->r0 * a[0]);
    y[1] = (y[1] & 0xffff) + (y[0] >> 16);
    y[2] = (y[2] & 0xffff) + (y[1] >> 16);
    y[0] = (y[0] & 0xffff) + c;
    y[1] = (y[1] & 0xffff) + (y[0] >> 16);
    y[2] = (y[2] & 0xffff) + (y[1] >> 16);
    rGen->r0 = y[0];
    rGen->r1 = y[1];
    rGen->r2 = y[2];
    return (rGen->r2 << 16) | rGen->r1;
}

/***************************************************************************
 *[@e
 *      Name:           UT_RandomNew
 *
 *      Usage:          Initialize a new random number generator.
 *
 *      Synopsis:       UT_RandomGen UT_RandomNew (UInt32 seed)
 *
 *      Description:    A new pseudo-random number generator is created and
 *                      initialized with seed value "seed".
 *
 *      Return value:   The identifier for the new random number generator.
 *                      If the generator could not be created, return NULL.
 *
 *      See also:
 *
 ***************************************************************************/

UT_RandomGen UT_RandomNew (UInt32 seed)
{
    RandomGen *pGen;

    pGen = (RandomGen *) UT_MemPerm (RandomGen);
    if (!pGen) {
        return NULL;
    }
    pGen->r0 = 0x330e;
    pGen->r1 = seed & 0xffff;
    pGen->r2 = seed >> 16;
    return pGen;
}

/***************************************************************************
 *[@e
 *      Name:           UT_RandomDelete
 *
 *      Usage:          Delete a random number generator.
 *
 *      Synopsis:       void UT_RandomDelete (UT_RandomGen rGen)
 *
 *      Description:    Delete the random number generator identified by "rGen".
 *
 *      Return value:   None.
 *
 *      See also:
 *
 ***************************************************************************/

void UT_RandomDelete (UT_RandomGen rGen)
{
    if (rGen) {
        UT_MemFree (rGen);
    }
}

/***************************************************************************
 *[@e
 *      Name:           UT_RandomGetInt32
 *
 *      Usage:          Generate a pseudo-random integer number.
 *
 *      Synopsis:       Int32 UT_RandomGetInt32 (UT_RandomGen rGen, Int32 rmin, Int32 rmax)
 *
 *      Description:    A random integer number is generated.
 *                      The generated numbers are uniformly distributed over the
 *                      interval [rmin, rmax).
 *                      Precondition: rmin < rmax
 *
 *      Return value:   The integer random number.
 *
 *      See also:
 *
 ***************************************************************************/

Int32 UT_RandomGetInt32 (UT_RandomGen rGen, Int32 rmin, Int32 rmax)
{
    Float64 num;

    num = (Float64)(0x3fffffff & calcRandom(rGen)) / (Float64)0x40000000;
    return (Int32) ((Float64)rmin + ((Float64)(rmax - rmin) * num));
}

UInt32 UT_RandomGetUInt32(UT_RandomGen rGen, UInt32 rmin, UInt32 rmax)
{
  Float64 num;
  num = (Float64)(0x3fffffff & calcRandom(rGen)) / (Float64)0x40000000;
  return (UInt32) ((Float64)rmin + ((Float64)(rmax - rmin) * num));
}

/***************************************************************************
 *[@e
 *      Name:           UT_RandomGetFloat32
 *
 *      Usage:          Generate a pseudeo-random floating point number.
 *
 *      Synopsis:       Float32 UT_RandomGetFloat32     (UT_RandomGen rGen,
 *                                               Float32 rmin, Float32 rmax)
 *
 *      Description:    A random floating point number is generated.
 *                      The generated numbers are uniformly distributed over the 
 *                      interval [rmin, rmax).
 *                      Precondition: rmin < rmax
 *
 *      Return value:   The floating point random number.
 *
 *      See also:
 *
 ***************************************************************************/

Float32 UT_RandomGetFloat32 (UT_RandomGen rGen, Float32 rmin, Float32 rmax)
{
    Float64 num;

    num = (Float64)(0x3fffffff & calcRandom(rGen)) / (Float64)0x40000000;
    return (Float32) (rmin + ((rmax - rmin) * num));
}
