/* pmath_rng.c
 *
 * This set of routines reproduces the Cray RANF family of random
 * number routines bit-for-bit.  Most of the routines here were
 * written by Jim Rathkopf, CP-Division, LLNL, and modified by
 * William C. Biester, B-Division. LLNL.  They were subsequently
 * modified for use in PMATH by Fred N. Fritsch, LC, LLNL.
 *
 * This version contains only the generator proper (equivalent to PM_RANF8
 * in the official PMATH version) and makes the seed and multiplier
 * directly accessible as arrays of doubles.  This was modified for use by
 * the Basis and Python module ranf.c, which contains its own seed and
 * multiplier management routines.  (FNF -- 9/4/96)
 */

/* -----------------------------------------------------------------------

Copyright (c) 1993,1994,1995,1996.
The Regents of the University of California.  All rights reserved.

This work was produced at the University of California, Lawrence
Livermore National Laboratory (UC LLNL) under contract no. W-7405-ENG-48
between the U.S. Department of Energy (DOE) and The Regents of the
University of California (University) for the operation of UC LLNL.
Copyright is reserved to the University for purposes of controlled
dissemination, commercialization through formal licensing, or other
disposition under terms of Contract 48; DOE policies, regulations and
orders; and U.S. statutes.

                         DISCLAIMER
 
This document was prepared as an account of work sponsored by an agency
of the United States Government.  Neither the United States Government
nor the University of California nor any of their employees, makes any
warranty, express or implied, or assumes any liability or responsibility
for the accuracy, completeness, or usefulness of any information,
apparatus, product, or process disclosed, or represents that its use
would not infringe privately-owned rights.  Reference herein to any
specific commercial products, process, or service by trade name,
trademark, manufacturer, or otherwise, does not necessarily constitute
or imply its endorsement, recommendation, or favoring by the United
States Government or the University of California.  The views and
opinions of authors expressed herein do not necessarily state or reflect
those of the United States Government or the University of California,
and shall not be used for advertising or product endorsement purposes.
*/

/* -----------------------------------------------------------------------
 *
 * Service routines defined:
 *    PM_16to24 - Convert 3 16-bit shorts to 2 24-bit doubles
 *    PM_24to16 - Convert 2 24-bit doubles to 3 16-bit shorts
 *    PM_GSeed  - Get the seed as two 24-bit doubles
 *    PM_SSeed  - Set the seed from two 24-bit doubles (unsafe)
 *    PM_GMult  - Get the multiplier as two 24-bit doubles
 *    PM_SMult  - Set the multiplier from two 24-bit doubles (unsafe)
 *
 * User-callable routines defined:
 *    PM_RANF   - The generator itself
 *
 * Currently responsible person:
 *    Fred N. Fritsch
 *    Computer Applications Organization
 *    LCPD, ICF Group
 *    Lawrence Livermore National Laboratory
 *    fnf@llnl.gov
 *
 * Modifications:
 * (See individual SLATEC-format prologues for detailed modification records.)
 *    03-09-93  Integrated Jim's routines into my library system.  (WCB)
 *    10-21-93  Modified to provide RANF8 for PMATH.  (FNF)
 *    10-22-93  Added code for RNFCNT.  (FNF)
 *    11-09-93  Added code for RNMSET.  (FNF)
 *    12-15-93  Modified to process bytes of seeds and multipliers in
 *              in the reverse order, to agree with the way they come
 *              from the Cray.  (FNF)
 *    03-15-94  Changed general pm_binds.h to pm_rnfset.h.  (FNF)
 *    04-25-94  Added SLATEC-format prologues. (FNF)
 *    09-26-95  Simplified code and improved internal documentation.
 *              This included removing unnecessary internal procedures
 *              setseed and seed48. (FNF)
 *    09-27-95  Simplified further by eliminating internal procedure
 *              rnset and no longer used variable Multiplier.  Also,
 *              moved PM_RANF8 so that it appears first. (FNF)
 *    09-28-95  Modified behavior of RNMSET to agree with MATHLIB. (FNF)
 *    10-02-95  Eliminated unions via calls to CV16TO64 and CV64TO16.
 *              Most conditional compilations are now in the rand48_*
 *              procedures. (FNF)
 *    10-05-95  Corrected "overflow" problem to make code work on Cray.
 *              Removed some unnecessary code in the process. (FNF)
 *    10-09-95  Removed conditional compiles.  (No longer needed after
 *              change of 10-02-95.) (FNF)
 *    08-27-96  Removed code for PM_RNSSET, PM_RNMSET, PM_RNSGET, PM_RNFCNT.
 *              Changed name PM_RANF8 to PM_RANF.  Added simple interface
 *              routines to get/set seed and multiplier.  Removed dependency
 *              on CV16TO64 and CV64TO16.  Eliminated Fortran bindings.  (FNF)
 *    09-04-96  Improved descriptions of the get/set routines.  (FNF)
 *    09-05-96  Added definitions of internal procedures PM_16to24 and
 *              PM_24to16 (see ranf.h).  (FNF)
 * -----------------------------------------------------------------------
 */



   /* Get includes, typedefs, and protypes for ranf.c and pmath_rng.c */
#include "ranf.h"

#define BASE24        16777216.0  /* 2^24 */

#define TWO8               256.0  /* 2^8 */
#define TWO16            65536.0  /* 2^16 */
#define TWO48  281474976710656.0  /* 2^48 */

/* ------------------------------------------------------
   The following static global variables retain the state
   of the random number generator between calls.
   ------------------------------------------------------ */
/* The default values to reproduce Cray values */
static double dseed[2] = {16555217.0, 9732691.0} ;
static double dmult[2] = {15184245.0, 2651554.0} ;
static double dadder = 0.0;    /* This is a relic from drand48. */


/* --------------------------------------
      Define in-line mod function.
   -------------------------------------- */
#define MODF(r,m)  (r - (floor(r/m))*m)


/* -----------------------------------------
      Define PMATH RNG internal procedures.
   -----------------------------------------
*/

/* PM_16to24 -- Take a number stored in 3 16-bit unsigned shorts and move it
 *              to 2 doubles each containing 24 bits of data.
 * 
 * This is a local function.
 *
 * Calling sequence:
 *   u16 x16[3];
 *   double x24[2];
 *   PM_16to24 (x16, x24);
 *
 *   x16 :   the 3 16-bit unsigned shorts.
 *   x24 :=  the 2 doubles of 24 bits returned.
 *
 * Return value:
 *   none.
 *
 * Note:
 *   This is the same as PMATH internal procedure rand48_16to24 except
 *   that the order of the entries in x16 have been reversed and its
 *   type declaration has been changed to u16.  It is intended for use
 *   by user-callable routines that set seeds or multipliers for users.
 */
void PM_16to24(u16 x16[3], double x24[2])
{
   double t0, t1, t2, t1u, t1l;

   t0 = (double) x16[0];
   t1 = (double) x16[1];
   t2 = (double) x16[2];

   t1u = floor(t1/TWO8);
   t1l = t1 -t1u*TWO8;
   x24[0] = t0 + t1l*TWO16;
   x24[1] = t1u + t2*TWO8;

#ifdef RAN_DEBUG
fprintf(stderr,"PM_16to24: x16 = %04x %04x %04x\n",x16[2],x16[1],x16[0]);
fprintf(stderr,"PM_16to24: x24 = %.1f %.1f\n",x24[1],x24[0]);
#endif
   return;
}


/* PM_24to16 -- Take a number stored in 2 doubles each containing 24 bits
 *              of data and move it to 3 16-bit unsigned shorts. 
 *
 * Calling sequence:
 *   double x24[2];
 *   u16 x16[3];
 *   PM_24to16 (x24, x16);
 *
 *   x24 :   the 2 doubles of 24 bits.
 *   x16 :=  the 3 16-bit unsigned shorts returned.
 *
 * Return value:
 *   none.
 *
 * Note:
 *   This is the same as PMATH internal procedure rand48_24to16 except
 *   that the order of the entries in x16 have been reversed and its
 *   type declaration has been changed to u16.  It is intended for use
 *   by user-callable routines that get seeds or multipliers for users.
 */
void PM_24to16(double x24[2], u16 x16[3])
{
   double t0u, t0l, t1u, t1l;

   t0u = floor(x24[0]/TWO16);
   t0l = x24[0] - t0u*TWO16;
   x16[0] = (u16)t0l;

   t1u = floor(x24[1]/TWO8);
   x16[2] = (u16)t1u;
   t1l = x24[1] - t1u*TWO8;
   x16[1] = (u16)(t0u + t1l*TWO8);

#ifdef RAN_DEBUG
fprintf(stderr,"PM_24to16: x24 = %.1f %.1f\n",x24[1],x24[0]);
fprintf(stderr,"PM_24to16: x16 = %04x %04x %04x\n",x16[2],x16[1],x16[0]);
#endif
   return;
}

/*   End of internal procedure definitions.
   -----------------------------------------
*/



/* --------------------------------------
      Define service routines.
   --------------------------------------
*/

/* PM_GSeed -- Get the seed as 2 doubles each containing 24 bits of data.
 * 
 * Calling sequence:
 *   double myseed[2];
 *   PM_GSeed(myseed);
 *
 *   myseed :=  the 2 doubles of 24 bits returned.
 *
 * Return value:
 *   none.
 */
void PM_GSeed(double seedout[2])
{
   seedout[0] = dseed[0];
   seedout[1] = dseed[1];

#ifdef RAN_DEBUG
   fprintf(stderr,"PM_GSeed: seed = %.1f %.1f\n", dseed[1], dseed[0]);
#endif
   return;
}

/* PM_SSeed -- Set the seed from 2 doubles each containing 24 bits of data.
 * 
 * Calling sequence:
 *   double myseed[2];
 *   PM_SSeed(myseed);
 *
 *   myseed :  the 2 doubles of 24 bits to be used as the new seed.
 *
 * Return value:
 *   none.
 *
 * Caution:
 *   This routine does not check for valid seeds!  It is intended for
 *   use in a safe interface such as Setranf (in ranf.c).
 *   The elements of the myseed array are the coefficients in the base
 *   2**24 representation of the odd 48-bit integer
 *         seed =  myseed[1]*2**24 + myseed[0],
 *   so myseed must satisfy:
 *     (1) myseed[0] and myseed[1] are integers, stored in doubles;
 *     (2) 0 < myseed[0] < 2**24, 0 <= myseed[1] < 2**24;
 *     (3) myseed[0] is odd.
 */
void PM_SSeed(double seedin[2])
{
   dseed[0] = seedin[0];
   dseed[1] = seedin[1];

#ifdef RAN_DEBUG
   fprintf(stderr,"PM_SSeed: seed = %.1f %.1f\n", dseed[1], dseed[0]);
#endif
   return;
}


/* PM_GMult -- Get the multiplier as 2 doubles each containing 24 bits 
 *             of data.
 * 
 * Calling sequence:
 *   double mymult[2];
 *   PM_GMult(mymult);
 *
 *   mymult :=  the 2 doubles of 24 bits returned.
 *
 * Return value:
 *   none.
 */
void PM_GMult(double multout[2])
{
   multout[0] = dmult[0];
   multout[1] = dmult[1];

#ifdef RAN_DEBUG
   fprintf(stderr,"PM_GMult: mult = %.1f %.1f\n", dmult[1], dmult[0]);
#endif
   return;
}

/* PM_SMult -- Set the multiplier from 2 doubles each containing 24 bits 
 *             of data.
 * 
 * Calling sequence:
 *   double mymult[2];
 *   PM_SMult(mymult);
 *
 *   mymult :  the 2 doubles of 24 bits to be used as the new multiplier.
 *
 * Return value:
 *   none.
 *
 * Caution:
 *   This routine does not check for valid multipliers!  It is intended for
 *   use in a safe interface such as Setmult (in ranf.c).
 *   The elements of the mymult array are the coefficients in the base
 *   2**24 representation of the odd 48-bit integer
 *         mult =  mymult[1]*2**24 + mymult[0].
 *   Since we require  1 < mult < 2**46, mymult must satisfy:
 *     (1) mymult[0] and mymult[1] are integers, stored in doubles;
 *     (2) 1 < mymult[0] < 2**24, 0 <= mymult[1] < 2**22;
 *     (3) mymult[0] is odd.
 */
void PM_SMult(double multin[2])
{
   dmult[0] = multin[0];
   dmult[1] = multin[1];

#ifdef RAN_DEBUG
   fprintf(stderr,"PM_SMult: mult = %.1f %.1f\n", dmult[1], dmult[0]);
#endif
   return;
}



/* --------------------------------------
      Define the generator itself.
   --------------------------------------
*/

/*DECK PM_RANF
C***BEGIN PROLOGUE  PM_RANF
C***PURPOSE  Uniform random-number generator.
C            The pseudorandom numbers generated by C function PM_RANF
C            are uniformly distributed in the open interval (0,1).
C***LIBRARY   PMATH
C***CATEGORY  L6A21
C***TYPE      REAL*8 (PM_RANF-8)
C***KEYWORDS  RANDOM NUMBER GENERATION, UNIFORM DISTRIBUTION
C***AUTHOR  Rathkopf, Jim, (LLNL/CP-Division)
C           Fritsch, Fred N., (LLNL/LC/MSS)
C***DESCRIPTION
C     (Portable version of Cray MATHLIB routine RANF.)
C     (Equivalent to Fortran-callable RANF8.)
C *Usage:
C        double r;
C        r = PM_RANF();
C
C *Function Return Values:
C     r        A random number between 0 and 1.
C
C *Description:
C     PM_RANF generates pseudorandom numbers lying strictly between 0
C     and 1. Each call to PM_RANF produces a different value, until the
C     sequence cycles after 2**46 calls.
C
C     PM_RANF is a linear congruential pseudorandom-number generator.
C     The default starting seed is
C                SEED = 4510112377116321(oct) = 948253fc9cd1(hex).
C     The multiplier is 1207264271730565(oct) = 2875a2e7b175(hex).
C
C *See Also:
C     This is the same generator used by the Fortran generator family
C     SRANF/DRANF/RANF8.
C     The starting seed for PM_RANF may be set via PM_SSeed.
C     The current PM_RANF seed may be obtained from PM_GSeed.
C     The current PM_RANF multiplier may be obtained from PM_GMult.
C     The PM_RANF multiplier may be set via PM_SMult (changing the
C     multiplier is not recommended).
C
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   930308  DATE WRITTEN
C           (Date from Biester's math_rnf.c.)
C   931021  Changed name from B_RANF to PM_RANF8.  (FNF)
C   940425  Added SLATEC-format prologue.  (FNF)
C   950922  Re-ordered code to eliminate old EMULRAND routine. (FNF)
C   951005  Eliminated unnecessary upper part computations.
C           Added fmod invocations in t2_48 computation to avoid
C           generating numbers too large for Cray floating point. (FNF)
C   951020  Replaced fmod with macro MODF to avoid library calls. (FNF)
C   951025  Removed MODF from first term of t2_48. (FNF)
C   960827  Changed name from PM_RANF8 to PM_RANF, eliminated reference
C           counter, and removed Fortran binding.  (FNF)
C***END PROLOGUE  PM_RANF
C
C*Internal notes:
C
C    This routine generates pseudo-random numbers through a 48-bit
C    linear congruential algorithm.  This emulates the drand48 library
C    of random number generators available on many, but not all,
C    UNIX machines.
C
C Algorithm:
C    x(n+1) = (a*x(n) + c)  mod m
C
C   where the defaults for the standard UNIX rand48 are:
C
C                   double name    hexdecimal          decimal
C    x: seed       -dseed[0],[1]  1234abcd330e      20017429951246
C    a: multiplier -dmult[0],[1]     5deece66d         25214903917
C    c: adder      -dadder                   b                  11
C    m: 2**48                    1000000000000     281474976710656
C
C     24-bit defaults (decimal) (lower bits listed first)
C    x: dseed[0],[1] = 13447950.0, 1193131.0
C    a: dmult[0],[1] = 15525485.0, 1502.0
C    c: dadder       = 11.0
C
C   The Cray defaults used in this code are:
C
C                   double name    hexdecimal           octal
C    x: seed       -dseed[0],[1]  948253fc9cd1    4510112377116321
C    a: multiplier -dmult[0],[1]  2875a2e7b175    1207264271730565
C    c: adder      -dadder                   0                   0
C    m: 2**48                    1000000000000   10000000000000000
C
C     24-bit defaults (decimal) (lower bits listed first)
C    x: dseed[0],[1] = 16555217.0, 9732691.0
C    a: dmult[0],[1] = 15184245.0, 2651554.0
C    c: dadder       = 0.0
C
C Return value:
C   double random number such that 0.0< d <1.0
C
C**End
 */
double PM_RANF ()
{
   double t1_48, t2_48, t1_24[2], t2_24;
   double d;

   /* perform 48-bit arithmetic using two part data */
   t1_48 = dmult[0]*dseed[0] + dadder;
   t2_48 = dmult[1]*dseed[0] + MODF(dmult[0]*dseed[1],BASE24);
       /* First term safe, since dmult[1] < 2**22 for default mult */

   t1_24[1] = floor(t1_48/BASE24);       /* upper part */
   t1_24[0] = t1_48 - t1_24[1]*BASE24;   /* lower part */

   t2_24 = MODF(t2_48, BASE24);      /* lower part */

   t2_48 = t2_24 + t1_24[1];

   t2_24 = MODF(t2_48, BASE24);      /* discard anything over 2**24 */

   d = (dseed[0] + dseed[1]*BASE24)/TWO48;

   dseed[0] = t1_24[0];
   dseed[1] = t2_24;

   return (d);
}