#include "stdafx.h"
#include "hnum_cblas.h"
/* -- LAPACK auxiliary routine (version 2.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
February 29, 1992
Purpose
=======
DLARAN returns a random real number from a uniform (0,1)
distribution.
Arguments
=========
ISEED (input/output) INT array, dimension (4)
On entry, the seed of the random number generator; the array
elements must be between 0 and 4095, and ISEED(4) must be
odd.
On exit, the seed is updated.
Further Details
===============
This routine uses a multiplicative congruential method with modulus
2**48 and multiplier 33952834046453 (see G.S.Fishman,
'Multiplicative congruential random number generators with modulus
2**b: an exhaustive analysis for b = 32 and a partial analysis for
b = 48', Math. Comp. 189, pp 331-344, 1990).
48-bit integers are stored in 4 integer array elements with 12 bits
per element. Hence the routine is portable across machines with
integers of 32 bits or more.
=====================================================================
*/
// =====================================================================
// Espen Harlinn:
// This implementation is a slightly modified version version
// of the code shipped with the SuperLU library, as it makes
// no sense to decrement the iseed pointer, and then add 1 to
// the array offsets to adjust for the decrement.
// =====================================================================
namespace harlinn
{
namespace numerics
{
namespace SuperLU
{
double dlaran_(int *iseed)
{
// multiply the seed by the multiplier modulo 2**48
int it4 = iseed[3] * 2549;
int it3 = it4 / 4096;
it4 -= it3 << 12;
it3 = it3 + iseed[2] * 2549 + iseed[3] * 2508;
int it2 = it3 / 4096;
it3 -= it2 << 12;
it2 = it2 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322;
int it1 = it2 / 4096;
it2 -= it1 << 12;
it1 = it1 + iseed[0] * 2549 + iseed[1] * 2508 + iseed[2] * 322 + iseed[3] * 494;
it1 %= 4096;
// return updated seed
iseed[0] = it1;
iseed[1] = it2;
iseed[2] = it3;
iseed[3] = it4;
// convert 48-bit integer to a real number in the interval (0,1)
return ((double) it1 + ((double) it2 + ((double) it3 + (double) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4;
}
};
};
};