mirror of
https://github.com/php/pecl-math-stats.git
synced 2026-03-23 23:22:18 +01:00
376 lines
12 KiB
C
376 lines
12 KiB
C
#include "randlib.h"
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],
|
|
Xlg1[32],Xlg2[32],Xa1vw,Xa2vw;
|
|
long Xqanti[32];
|
|
void advnst(long k)
|
|
/*
|
|
**********************************************************************
|
|
void advnst(long k)
|
|
ADV-a-N-ce ST-ate
|
|
Advances the state of the current generator by 2^K values and
|
|
resets the initial seed to that value.
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Advance_State from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
k -> The generator is advanced by2^K values
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
|
|
long g,i,ib1,ib2;
|
|
long qrgnin;
|
|
/*
|
|
Abort unless random number generator initialized
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(qrgnin) goto S10;
|
|
fputs(" ADVNST called before random generator initialized - ABORT\n",
|
|
stderr);
|
|
exit(1);
|
|
S10:
|
|
gscgn(0L,&g);
|
|
ib1 = Xa1;
|
|
ib2 = Xa2;
|
|
for(i=1; i<=k; i++) {
|
|
ib1 = mltmod(ib1,ib1,Xm1);
|
|
ib2 = mltmod(ib2,ib2,Xm2);
|
|
}
|
|
setsd(mltmod(ib1,*(Xcg1+g-1),Xm1),mltmod(ib2,*(Xcg2+g-1),Xm2));
|
|
/*
|
|
NOW, IB1 = A1**K AND IB2 = A2**K
|
|
*/
|
|
#undef numg
|
|
}
|
|
void getsd(long *iseed1,long *iseed2)
|
|
/*
|
|
**********************************************************************
|
|
void getsd(long *iseed1,long *iseed2)
|
|
GET SeeD
|
|
Returns the value of two integer seeds of the current generator
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Get_State from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
iseed1 <- First integer seed of generator G
|
|
iseed2 <- Second integer seed of generator G
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xcg1[],Xcg2[];
|
|
long g;
|
|
long qrgnin;
|
|
/*
|
|
Abort unless random number generator initialized
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(qrgnin) goto S10;
|
|
fprintf(stderr,"%s\n",
|
|
" GETSD called before random number generator initialized -- abort!");
|
|
exit(0);
|
|
S10:
|
|
gscgn(0L,&g);
|
|
*iseed1 = *(Xcg1+g-1);
|
|
*iseed2 = *(Xcg2+g-1);
|
|
#undef numg
|
|
}
|
|
long ignlgi(void)
|
|
/*
|
|
**********************************************************************
|
|
long ignlgi(void)
|
|
GeNerate LarGe Integer
|
|
Returns a random integer following a uniform distribution over
|
|
(1, 2147483562) using the current generator.
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Random from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gssst(long getset,long *qset);
|
|
extern void gscgn(long getset,long *g);
|
|
extern void inrgcm(void);
|
|
extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
|
|
extern long Xqanti[];
|
|
long ignlgi,curntg,k,s1,s2,z;
|
|
long qqssd,qrgnin;
|
|
/*
|
|
IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
|
|
IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
|
|
THIS ROUTINE 2) A CALL TO SETALL.
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(!qrgnin) inrgcm();
|
|
gssst(0,&qqssd);
|
|
if(!qqssd) setall(1234567890L,123456789L);
|
|
/*
|
|
Get Current Generator
|
|
*/
|
|
gscgn(0L,&curntg);
|
|
s1 = *(Xcg1+curntg-1);
|
|
s2 = *(Xcg2+curntg-1);
|
|
k = s1/53668L;
|
|
s1 = Xa1*(s1-k*53668L)-k*12211;
|
|
if(s1 < 0) s1 += Xm1;
|
|
k = s2/52774L;
|
|
s2 = Xa2*(s2-k*52774L)-k*3791;
|
|
if(s2 < 0) s2 += Xm2;
|
|
*(Xcg1+curntg-1) = s1;
|
|
*(Xcg2+curntg-1) = s2;
|
|
z = s1-s2;
|
|
if(z < 1) z += (Xm1-1);
|
|
if(*(Xqanti+curntg-1)) z = Xm1-z;
|
|
ignlgi = z;
|
|
return ignlgi;
|
|
#undef numg
|
|
}
|
|
void initgn(long isdtyp)
|
|
/*
|
|
**********************************************************************
|
|
void initgn(long isdtyp)
|
|
INIT-ialize current G-e-N-erator
|
|
Reinitializes the state of the current generator
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Init_Generator from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
isdtyp -> The state to which the generator is to be set
|
|
isdtyp = -1 => sets the seeds to their initial value
|
|
isdtyp = 0 => sets the seeds to the first value of
|
|
the current block
|
|
isdtyp = 1 => sets the seeds to the first value of
|
|
the next block
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
|
|
long g;
|
|
long qrgnin;
|
|
/*
|
|
Abort unless random number generator initialized
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(qrgnin) goto S10;
|
|
fprintf(stderr,"%s\n",
|
|
" INITGN called before random number generator initialized -- abort!");
|
|
exit(1);
|
|
S10:
|
|
gscgn(0L,&g);
|
|
if(-1 != isdtyp) goto S20;
|
|
*(Xlg1+g-1) = *(Xig1+g-1);
|
|
*(Xlg2+g-1) = *(Xig2+g-1);
|
|
goto S50;
|
|
S20:
|
|
if(0 != isdtyp) goto S30;
|
|
goto S50;
|
|
S30:
|
|
/*
|
|
do nothing
|
|
*/
|
|
if(1 != isdtyp) goto S40;
|
|
*(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
|
|
*(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
|
|
goto S50;
|
|
S40:
|
|
fprintf(stderr,"%s\n","isdtyp not in range in INITGN");
|
|
exit(1);
|
|
S50:
|
|
*(Xcg1+g-1) = *(Xlg1+g-1);
|
|
*(Xcg2+g-1) = *(Xlg2+g-1);
|
|
#undef numg
|
|
}
|
|
void inrgcm(void)
|
|
/*
|
|
**********************************************************************
|
|
void inrgcm(void)
|
|
INitialize Random number Generator CoMmon
|
|
Function
|
|
Initializes common area for random number generator. This saves
|
|
the nuisance of a BLOCK DATA routine and the difficulty of
|
|
assuring that the routine is loaded with the other routines.
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
|
|
extern long Xqanti[];
|
|
long T1;
|
|
long i;
|
|
/*
|
|
V=20; W=30;
|
|
A1W = MOD(A1**(2**W),M1) A2W = MOD(A2**(2**W),M2)
|
|
A1VW = MOD(A1**(2**(V+W)),M1) A2VW = MOD(A2**(2**(V+W)),M2)
|
|
If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
|
|
An efficient way to precompute a**(2*j) MOD m is to start with
|
|
a and square it j times modulo m using the function MLTMOD.
|
|
*/
|
|
Xm1 = 2147483563L;
|
|
Xm2 = 2147483399L;
|
|
Xa1 = 40014L;
|
|
Xa2 = 40692L;
|
|
Xa1w = 1033780774L;
|
|
Xa2w = 1494757890L;
|
|
Xa1vw = 2082007225L;
|
|
Xa2vw = 784306273L;
|
|
for(i=0; i<numg; i++) *(Xqanti+i) = 0;
|
|
T1 = 1;
|
|
/*
|
|
Tell the world that common has been initialized
|
|
*/
|
|
gsrgs(1L,&T1);
|
|
#undef numg
|
|
}
|
|
void setall(long iseed1,long iseed2)
|
|
/*
|
|
**********************************************************************
|
|
void setall(long iseed1,long iseed2)
|
|
SET ALL random number generators
|
|
Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
|
|
initial seeds of the other generators are set accordingly, and
|
|
all generators states are set to these seeds.
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Set_Initial_Seed from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
iseed1 -> First of two integer seeds
|
|
iseed2 -> Second of two integer seeds
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gssst(long getset,long *qset);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
|
|
long T1;
|
|
long g,ocgn;
|
|
long qrgnin;
|
|
T1 = 1;
|
|
/*
|
|
TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
|
|
HAS BEEN CALLED.
|
|
*/
|
|
gssst(1,&T1);
|
|
gscgn(0L,&ocgn);
|
|
/*
|
|
Initialize Common Block if Necessary
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(!qrgnin) inrgcm();
|
|
*Xig1 = iseed1;
|
|
*Xig2 = iseed2;
|
|
initgn(-1L);
|
|
for(g=2; g<=numg; g++) {
|
|
*(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
|
|
*(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
|
|
gscgn(1L,&g);
|
|
initgn(-1L);
|
|
}
|
|
gscgn(1L,&ocgn);
|
|
#undef numg
|
|
}
|
|
void setant(long qvalue)
|
|
/*
|
|
**********************************************************************
|
|
void setant(long qvalue)
|
|
SET ANTithetic
|
|
Sets whether the current generator produces antithetic values. If
|
|
X is the value normally returned from a uniform [0,1] random
|
|
number generator then 1 - X is the antithetic value. If X is the
|
|
value normally returned from a uniform [0,N] random number
|
|
generator then N - 1 - X is the antithetic value.
|
|
All generators are initialized to NOT generate antithetic values.
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Set_Antithetic from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
qvalue -> nonzero if generator G is to generating antithetic
|
|
values, otherwise zero
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xqanti[];
|
|
long g;
|
|
long qrgnin;
|
|
/*
|
|
Abort unless random number generator initialized
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(qrgnin) goto S10;
|
|
fprintf(stderr,"%s\n",
|
|
" SETANT called before random number generator initialized -- abort!");
|
|
exit(1);
|
|
S10:
|
|
gscgn(0L,&g);
|
|
Xqanti[g-1] = qvalue;
|
|
#undef numg
|
|
}
|
|
void setsd(long iseed1,long iseed2)
|
|
/*
|
|
**********************************************************************
|
|
void setsd(long iseed1,long iseed2)
|
|
SET S-ee-D of current generator
|
|
Resets the initial seed of the current generator to ISEED1 and
|
|
ISEED2. The seeds of the other generators remain unchanged.
|
|
This is a transcription from Pascal to Fortran of routine
|
|
Set_Seed from the paper
|
|
L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
|
|
with Splitting Facilities." ACM Transactions on Mathematical
|
|
Software, 17:98-111 (1991)
|
|
Arguments
|
|
iseed1 -> First integer seed
|
|
iseed2 -> Second integer seed
|
|
**********************************************************************
|
|
*/
|
|
{
|
|
#define numg 32L
|
|
extern void gsrgs(long getset,long *qvalue);
|
|
extern void gscgn(long getset,long *g);
|
|
extern long Xig1[],Xig2[];
|
|
long g;
|
|
long qrgnin;
|
|
/*
|
|
Abort unless random number generator initialized
|
|
*/
|
|
gsrgs(0L,&qrgnin);
|
|
if(qrgnin) goto S10;
|
|
fprintf(stderr,"%s\n",
|
|
" SETSD called before random number generator initialized -- abort!");
|
|
exit(1);
|
|
S10:
|
|
gscgn(0L,&g);
|
|
*(Xig1+g-1) = iseed1;
|
|
*(Xig2+g-1) = iseed2;
|
|
initgn(-1L);
|
|
#undef numg
|
|
}
|