Commit: ffb6fce2383e062b18e2fa7ddfb23a4b010c0c1c
Parent: efaa778d96864408dc9dd999f9c4f2ddf660460d
Author: Randy Palamar
Date: Tue, 26 Mar 2024 12:25:32 -0600
replace slow garbage rng with xoroshiro128+
Diffstat:
| M | Photon.cpp | | | 35 | +++++++---------------------------- |
| M | main.cpp | | | 2 | ++ |
| M | utils.cpp | | | 110 | ++++++++++++++++++++----------------------------------------------------------- |
| M | utils.h | | | 11 | +++++------ |
4 files changed, 41 insertions(+), 117 deletions(-)
diff --git a/Photon.cpp b/Photon.cpp
@@ -58,14 +58,9 @@ void Photon::Launch()
void Photon::CalStep()
{
- double ran;
- if (s<ZERO)
- {
- do
- {
- ran=RandomNum();
- } while(ran<=0||ran>1);
- s=-log(ran); //dimensionless,
+ if (s < ZERO) {
+ double ran = random_uniform();
+ s = -log(ran + EPS);
}
}
@@ -119,10 +114,7 @@ void Photon::Scatter() // calculate next propagation direction
if (!dead)
{
double costheta,fei,ran;
- do
- {
- ran=RandomNum();
- } while(ran<0||ran>=1);
+ ran = random_uniform();
if (g!=0)
{
double aa=(1-g*g)/(1-g+2*g*ran);
@@ -136,10 +128,7 @@ void Photon::Scatter() // calculate next propagation direction
}//(3.28)
double sintheta,sinfei,cosfei;
sintheta=sqrt(1-costheta*costheta);
- do
- {
- ran=RandomNum();
- } while(ran<0||ran>=1);
+ ran = random_uniform();
fei=2*PI*ran;//(3.29)
cosfei=cos(fei);
if (fei<PI)
@@ -193,12 +182,7 @@ void Photon::ReflectTransmit() // deal with photon across boundary
B=B*B;
ralphai=0.5*(2*A-A*A-A*B)/(B*(1-A));//(3.36)
}
- double ran;
- do
- {
- ran=RandomNum();
- }
- while(ran<0||ran>=1);
+ double ran = random_uniform();
if (ran<=ralphai) //rebound to current layer
{
uz=-uz;
@@ -228,11 +212,7 @@ void Photon::Termination() //Russian roulette survial test
double m=10;
if (w<=wth)
{
- double e;
- do
- {
- e=RandomNum();
- }while(e<0||e>1);
+ double e = random_uniform();
if (e>(1/m))
{
w=0;
@@ -245,4 +225,3 @@ void Photon::Termination() //Russian roulette survial test
}
}
}
-
diff --git a/main.cpp b/main.cpp
@@ -30,6 +30,8 @@ int main(void)
cout<<"*****************************************"<<endl;
cout<<endl;
+ random_init();
+
cout<<"File name(Input:*.mci; Output:*.mco): ";
string sFilename, sInfilename, sOutfilename;
cin>>sFilename;
diff --git a/utils.cpp b/utils.cpp
@@ -1,93 +1,37 @@
+#include <stdint.h>
+#include <time.h>
+
+typedef uint64_t u64;
+
#include "utils.h"
-/***********************************************************
-* A random number generator from Numerical Recipes in C. ****/
-#define MBIG 1000000000
-#define MSEED 161803398
-#define MZ 0
-#define FAC 1.0E-9
-float ran3(int &idum)
+static u64 rand_state[2];
+
+void
+random_init(void)
{
- static int inext,inextp;
- static long ma[56];
- static int iff=0;
- long mj,mk;
- int i,ii,k;
- if (idum < 0 || iff == 0)
- {
- iff=1;
- mj=MSEED-(idum < 0 ? -idum : idum);
- mj %= MBIG;
- ma[55]=mj;
- mk=1;
- for (i=1;i<=54;i++)
- {
- ii=(21*i) % 55;
- ma[ii]=mk;
- mk=mj-mk;
- if (mk < MZ) mk += MBIG;
- mj=ma[ii];
- }
- for (k=1;k<=4;k++)
- for (i=1;i<=55;i++)
- {
- ma[i] -= ma[1+(i+30) % 55];
- if (ma[i] < MZ) ma[i] += MBIG;
- }
- inext=0;
- inextp=31;
- idum=1;
- }
- if (++inext == 56)
- inext=1;
- if (++inextp == 56)
- inextp=1;
- mj=ma[inext]-ma[inextp];
- if (mj < MZ) mj += MBIG;
- ma[inext]=mj;
- return mj*FAC;
+ struct timespec ts;
+ clock_gettime(CLOCK_REALTIME, &ts);
+ rand_state[0] = (intptr_t)&printf ^ ts.tv_sec
+ ^ ((u64)ts.tv_nsec * 0xAC5533CD);
+ rand_state[1] = (intptr_t)&malloc ^ ts.tv_sec
+ ^ ((u64)ts.tv_nsec * 0xAC5533CD);
}
-#undef MBIG
-#undef MSEED
-#undef MZ
-#undef FAC
-
-/***********************************************************
- * Generate a random number between 0 and 1. Take a
- * number as seed the first time entering the function.
- * The seed is limited to 1<<15.
- * Wang et al found that when idum is too large, ran3 may return
- * numbers beyond 0 and 1.
- ****/
-//#define STANDARDTEST 0 /* testing program using fixed rnd seed. */
-double RandomNum(void)
+static u64
+xoroshiro128plus(u64 s[2])
{
- static bool first_time=1;
- static int idum; /* seed for ran3. */
- if(first_time)
- {
- #if STANDARDTEST /* Use fixed seed to test the program. */
- idum = - 1;
- #else
- idum = -(int)time(NULL)%(1<<15); /* use 16-bit integer as the seed. */
- #endif
- ran3(idum);
- first_time = 0;
- idum = 1;
- }
- return( (double)ran3(idum) );
+ u64 s0 = s[0];
+ u64 s1 = s[1];
+ u64 result = s0 + s1;
+ s1 ^= s0;
+ s[0] = ((s0 << 24) | (s0 >> 40)) ^ s1 ^ (s1 << 16);
+ s[1] = (s1 << 37) | (s1 >> 27);
+ return result;
}
-/////////////////////////////////////////////
-int min(int d1,int d2)
+double
+random_uniform(void)
{
- if (d1<d2)
- {
- return d1;
- }
- else
- {
- return d2;
- }
+ return xoroshiro128plus(rand_state) / (double)UINT64_MAX;
}
diff --git a/utils.h b/utils.h
@@ -21,6 +21,8 @@
#define BIGVAL 1e6//big value
#endif
+#define EPS 1.0e-9
+
//global variables
#ifdef GLOBALORIGIN
#define GLOBAL
@@ -46,13 +48,10 @@ GLOBAL double x_offset;
GLOBAL double y_offset;
//generate random data
-float ran3(int *idum);
-double RandomNum(void);
-
+void random_init(void);
+double random_uniform(void);
-//other auxilliary functions
-int min(int d1, int d2);
-#define SIGN(x) ((x)>=0?1:-1)
+#define SIGN(x) ((x) >= 0 ? 1 : -1)
////////allocte dynamic array/////////////////////////