CARDAMOM-framework · aantonyb · Sep 14, 2026 · Sep 14, 2026 · Sep 14, 2026 · Sep 14, 2026
diff --git a/C/auxi_fun/cardamom_random.h b/C/auxi_fun/cardamom_random.h
@@ -0,0 +1,50 @@
+#ifndef CARDAMOM_RANDOM_H
+#define CARDAMOM_RANDOM_H
+
+#include <stdint.h>
+
+typedef struct {
+    uint64_t state;
+    uint64_t inc;
+} pcg32_random_t;
+
+#define CARDAMOM_RAND_MAX 2147483647L
+
+/* Global state */
+static pcg32_random_t pcg32_global = {0x853c49e6748fea9bULL, 0xda3e39cb94b95bdbULL};
+
+/* PCG32 core algorithm */
+static inline uint32_t pcg32_random_r(pcg32_random_t* rng) {
+    uint64_t oldstate = rng->state;
+    rng->state = oldstate * 6364136223846793005ULL + rng->inc;
+    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
+    uint32_t rot = oldstate >> 59u;
+    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
+}
+
+static inline void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq) {
+    rng->state = 0U;
+    rng->inc = (initseq << 1u) | 1u;
+    pcg32_random_r(rng);
+    rng->state += initstate;
+    pcg32_random_r(rng);
+}
+
+/* Public API */
+static inline void cardarand_seed(uint64_t seed) {
+    pcg32_srandom_r(&pcg32_global, seed, 0xda3e39cb94b95bdbULL);
+}
+
+static inline long cardarand(void) {
+    return (long)(pcg32_random_r(&pcg32_global) & 0x7FFFFFFFL);
+}
+
+static inline double cardaurand(void) {
+    return (double)cardarand() / (double)CARDAMOM_RAND_MAX;
+}
+
+static inline double cardarand_uniform(void) {
+    return (double)pcg32_random_r(&pcg32_global) / (double)0x100000000ULL;
+}
+
+#endif
diff --git a/C/auxi_fun/seedrandomnumber.c b/C/auxi_fun/seedrandomnumber.c
@@ -1,11 +1,12 @@
 
 #pragma once
-#include <stdlib.h>  // rand(), srand()
+#include <stdlib.h>
 #include <math.h>
-#include <time.h> //for time
+#include <time.h>
 #include <stdio.h>
 #include <string.h>
-#include <stdint.h>  // For explicit integer types
+#include <stdint.h>
+#include "cardamom_random.h"
 /*triple-seeding in simple C code*/
 
 
@@ -59,13 +60,14 @@ int charseed_filename(const char *charinput){
 int seedrandomnumber(const char *charinput){
 
     long seed1=(long)charseed(charinput);
-    
+
     // unsigned int seed1 = (unsigned int)charseed(charinput);
     // unsigned int seed2 = (unsigned int)time(NULL);
     // unsigned int seed3 = (unsigned int)(uintptr_t)&charinput;
 
     // srand(seed1 + seed2 + seed3);
+    // srandom(seed1);
 
-    srandom(seed1);
+    cardarand_seed((uint64_t)seed1);
 
 return 0;}
diff --git a/C/math_fun/randn.c b/C/math_fun/randn.c
@@ -1,9 +1,10 @@
 #pragma once
+#include "../auxi_fun/cardamom_random.h"
 double randn(){
 
 double pi=3.141592653589793;
-double r1=(double)random()/(double)RAND_MAX;
-double r2=(double)random()/(double)RAND_MAX;
+double r1=cardaurand();
+double r2=cardaurand();
 
 
 double rn=sqrt(-2*log(r1)) * cos(2*pi*r2);

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/ADEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/ADEMCMC.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DEMCMC.c"
@@ -100,7 +101,7 @@ for (n=0;n<PI.npars;n++){
 
 if (MCO.randparini==1 && PI.parfix[n]!=1){
 /*random parameter if PI.parini = -9999*/
-PARS[n + nn * PI.npars] = nor2par((double)random() / (double)RAND_MAX, PI.parmin[n], PI.parmax[n]);}
+PARS[n + nn * PI.npars] = nor2par((double)cardarand() / (double)CARDAMOM_RAND_MAX, PI.parmin[n], PI.parmax[n]);}
 else{par=PI.parini[n+nn*PI.npars];
 PARS[n+nn*PI.npars]=par;
 if (par>PI.parmax[n] | par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range");}
@@ -159,15 +160,15 @@ for ( ;N.ITER<MCO.nOUT;N.ITER++){
 	else {
 	/*Step size is 1 wigth 10% prob iterations*/
         PI.stepsize[0] = 1 - (1 - 2.38 / sqrt(2 * PI.npars) * 0.1) * 
-                 (double)(( (double)random() / (double)RAND_MAX ) < 0.9);
+                 (double)(( (double)cardarand() / (double)CARDAMOM_RAND_MAX ) < 0.9);
 	/*take a step (DE-MCMC style)*/
 	PI.stepsize[0]=PI.stepsize[0];
 	withinrange=STEP_DEMCMC(PARS,pars_new,PI,nn,NC);
 	gratio=0;
 	}
 
 
-	lr = log((double)random() / (double)RAND_MAX);
+	lr = log((double)cardarand() / (double)CARDAMOM_RAND_MAX);
 	/*p(x) = 0 if parameters outside bounds*/
 	if (withinrange==1 & -P[nn]+gratio>lr){
 wrlocal=wrlocal+1;

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/ADEMCMC_301.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/ADEMCMC_301.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DEMCMC.c"
@@ -101,7 +102,7 @@ for (n=0;n<PI.npars;n++){
 
 if (MCO.randparini==1 && PI.parfix[n]!=1){
 /*random parameter if PI.parini = -9999*/
-PARS[n+nn*PI.npars]=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+PARS[n+nn*PI.npars]=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 else{par=PI.parini[n+nn*PI.npars];
 PARS[n+nn*PI.npars]=par;
 if (par>PI.parmax[n] | par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range");}
@@ -158,15 +159,15 @@ for (N.ITER=0;N.ITER<MCO.nOUT;N.ITER++){
 	//Standard DEMCMC
 	else {
 	/*Step size is 1 wigth 10% prob iterations*/
-        PI.stepsize[0]=1 - (1-2.38/sqrt(2*PI.npars)/10)*(double)((double)(random()/(double)RAND_MAX)<0.9);
+        PI.stepsize[0]=1 - (1-2.38/sqrt(2*PI.npars)/10)*(double)((double)(cardarand()/(double)CARDAMOM_RAND_MAX)<0.9);
 	/*take a step (DE-MCMC style)*/
 	PI.stepsize[0]=PI.stepsize[0];
 	withinrange=STEP_DEMCMC(PARS0,pars_new,PI,nn,NC);
 	gratio=0;
 	}
 
 
-	lr=log((double)random()/(double)RAND_MAX);
+	lr=log(cardaurand());
 	/*p(x) = 0 if parameters outside bounds*/
 	if (withinrange==1 & -P[nn]+gratio>lr){
 wrlocal=wrlocal+1;

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/AFDEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/AFDEMCMC.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DEMCMC.c"
@@ -100,7 +101,7 @@ for (n=0;n<PI.npars;n++){
 
 if (MCO.randparini==1 && PI.parfix[n]!=1){
 /*random parameter if PI.parini = -9999*/
-PARS[n + nn * PI.npars] = nor2par((double)random() / (double)RAND_MAX, PI.parmin[n], PI.parmax[n]);}
+PARS[n + nn * PI.npars] = nor2par((double)cardarand() / (double)CARDAMOM_RAND_MAX, PI.parmin[n], PI.parmax[n]);}
 else{par=PI.parini[n+nn*PI.npars];
 PARS[n+nn*PI.npars]=par;
 if (par>PI.parmax[n] | par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range");}
@@ -159,15 +160,15 @@ for ( ;N.ITER<MCO.nOUT;N.ITER++){
 	else {
 	/*Step size is 1 wigth 10% prob iterations*/
         PI.stepsize[0] = 1 - (1 - 2.38 / sqrt(2 * PI.npars) * 0.1) * 
-                 (double)(( (double)random() / (double)RAND_MAX ) < 0.9);
+                 (double)(( (double)cardarand() / (double)CARDAMOM_RAND_MAX ) < 0.9);
 	/*take a step (DE-MCMC style)*/
 	PI.stepsize[0]=PI.stepsize[0];
 	withinrange=STEP_DEMCMC(PARS,pars_new,PI,nn,NC);
 	gratio=0;
 	}
 
 
-	lr = log((double)random() / (double)RAND_MAX);
+	lr = log((double)cardarand() / (double)CARDAMOM_RAND_MAX);
 	/*p(x) = 0 if parameters outside bounds*/
 	if (withinrange==1 & -P[nn]+gratio>lr){
 wrlocal=wrlocal+1;

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/AFDEMCMCZS.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/AFDEMCMCZS.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_AFDEMCMC.c"
@@ -63,7 +64,7 @@ double par;
 for (nn=0;nn<NC;nn++){
 for (n=0;n<PI.npars;n++){
 if (MCO.randparini==1 && PI.parfix[n]!=1){
-par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);
+par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);
 }else{
 par=PI.parini[n+nn*PI.npars];
 if (par>PI.parmax[n] || par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range\n");}
@@ -76,7 +77,7 @@ Z[nn*PI.npars+n]=par;
 for (nn=NC;nn<M0;nn++){
 for (n=0;n<PI.npars;n++){
 if (PI.parfix[n]==1){par=PI.parini[n];}
-else{par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+else{par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 Z[nn*PI.npars+n]=par;
 }}
 
@@ -163,14 +164,14 @@ gratio=0;
 if (N.ITER<switch_iter){
 withinrange=STEP_AFDEMCMC(PARS,pars_new,PI,nn,activeNC,&gratio);
 }else{
-if ((double)random()/(double)RAND_MAX<psnooker){
+if (cardaurand()<psnooker){
 withinrange=STEP_DEMCMCZ_SNOOKER(&PARS[nn*PI.npars],Z,M,pars_new,PI,&gratio);
 }else{
 withinrange=STEP_DEMCMCZ_PARALLEL(&PARS[nn*PI.npars],Z,M,pars_new,PI);
 }
 }
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 
 if (withinrange==1){
 wrlocal=wrlocal+1;

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMC.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DEMCMC.c"
@@ -97,7 +98,7 @@ for (n=0;n<PI.npars;n++){
 
 if (MCO.randparini==1 && PI.parfix[n]!=1){
 /*random parameter if PI.parini = -9999*/
-PARS[n + nn * PI.npars] = nor2par((double)random() / (double)RAND_MAX, PI.parmin[n], PI.parmax[n]);}
+PARS[n + nn * PI.npars] = nor2par((double)cardarand() / (double)CARDAMOM_RAND_MAX, PI.parmin[n], PI.parmax[n]);}
 else
 
 /*{PARS[n+nn*PI.npars]=PI.parini[n+nn*PI.npars];}}}
@@ -141,7 +142,7 @@ for (N.ITER=0;N.ITER<MCO.nOUT;N.ITER++){
 
 	/*Step size is 1 wigth 10% prob iterations*/
         PI.stepsize[0] = 1 - (1 - 2.38 / sqrt(2 * PI.npars) * 0.1) * 
-                 (double)(( (double)random() / (double)RAND_MAX ) < 0.9);
+                 (double)(( (double)cardarand() / (double)CARDAMOM_RAND_MAX ) < 0.9);
 	/*take a step (DE-MCMC style)*/
 	//PI.stepsize[0]=PI.stepsize[0]/10;
 	withinrange=STEP_DEMCMC(PARS,pars_new,PI,nn,NC);
@@ -157,7 +158,7 @@ wrlocal=wrlocal+1;
 	*/
 	/*treating nans as -inf*/
 	if (isnan(P_new)){P_new=log(0);}
-	if (P_new - P[nn] > log((double)random() / (double)RAND_MAX)) { N.ACC = N.ACC + 1;
+	if (P_new - P[nn] > log((double)cardarand() / (double)CARDAMOM_RAND_MAX)) { N.ACC = N.ACC + 1;
 	if (isinf(P_new)==0 && isinf(P[nn])){printf("pnew = %2.1f, p = %2.1f, (P_new-P[nn]) = %2.1f\n",P_new,P[nn],P_new-P[nn]);}
 
 

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMCZS.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMCZS.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DEMCMCZS.c"
@@ -60,7 +61,7 @@ double par;
 for (nn=0;nn<NC;nn++){
 for (n=0;n<PI.npars;n++){
 if (MCO.randparini==1 && PI.parfix[n]!=1){
-par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);
+par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);
 }else{
 par=PI.parini[n+nn*PI.npars];
 if (par>PI.parmax[n] || par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range\n");}
@@ -74,7 +75,7 @@ X[nn*PI.npars+n]=par;
 for (nn=NC;nn<M0;nn++){
 for (n=0;n<PI.npars;n++){
 if (PI.parfix[n]==1){par=PI.parini[n];}
-else{par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+else{par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 Z[nn*PI.npars+n]=par;
 }}
 
@@ -102,13 +103,13 @@ for (N.ITER=0;N.ITER<MCO.nOUT;N.ITER++){
 for (nn=0;nn<NC;nn++){
 
 gratio=0;
-if ((double)random()/(double)RAND_MAX<psnooker){
+if (cardaurand()<psnooker){
 withinrange=STEP_DEMCMCZ_SNOOKER(&X[nn*PI.npars],Z,M,pars_new,PI,&gratio);
 }else{
 withinrange=STEP_DEMCMCZ_PARALLEL(&X[nn*PI.npars],Z,M,pars_new,PI);
 }
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 
 if (withinrange==1){
 wrlocal=wrlocal+1;

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMCZS_WARMUP.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/DEMCMCZS_WARMUP.c
@@ -41,7 +41,7 @@ X[nn*PI.npars+n]=parini[n+nn*PI.npars];
 for (nn=NC;nn<M0;nn++){
 for (n=0;n<PI.npars;n++){
 if (PI.parfix[n]==1){Z[nn*PI.npars+n]=parini[n];}
-else{Z[nn*PI.npars+n]=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+else{Z[nn*PI.npars+n]=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 }}
 
 for (nn=0;nn<NC;nn++){
@@ -57,13 +57,13 @@ for (iter=0;iter<niter;iter++){
 for (nn=0;nn<NC;nn++){
 
 gratio=0;
-if ((double)random()/(double)RAND_MAX<psnooker){
+if (cardaurand()<psnooker){
 withinrange=STEP_DEMCMCZ_SNOOKER(&X[nn*PI.npars],Z,M,pars_new,PI,&gratio);
 }else{
 withinrange=STEP_DEMCMCZ_PARALLEL(&X[nn*PI.npars],Z,M,pars_new,PI);
 }
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 
 if (withinrange==1){
 P_new=MODEL_LIKELIHOOD(DATA,pars_new);

diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/DREAMZS.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/DREAMZS.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "../../../auxi_fun/cardamom_random.h"
 #include "../../../math_fun/std.c"
 #include "NORMPARS.c"
 #include "STEP_DREAMZS.c"
@@ -55,7 +56,7 @@ double par;
 for (nn=0;nn<NC;nn++){
 for (n=0;n<PI.npars;n++){
 if (MCO.randparini==1 && PI.parfix[n]!=1){
-par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);
+par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);
 }else{
 par=PI.parini[n+nn*PI.npars];
 if (par>PI.parmax[n] || par<PI.parmin[n]){printf("Warning, prescribed initial parameters are out of range\n");}
@@ -68,7 +69,7 @@ X[nn*PI.npars+n]=par;
 for (nn=NC;nn<M0;nn++){
 for (n=0;n<PI.npars;n++){
 if (PI.parfix[n]==1){par=PI.parini[n];}
-else{par=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+else{par=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 Z[nn*PI.npars+n]=par;
 }}
 
@@ -99,18 +100,18 @@ for ( ;N.ITER<MCO.nOUT;N.ITER++){
 for (nn=0;nn<NC;nn++){
 
 gratio=0;
-if ((double)random()/(double)RAND_MAX<psnooker){
+if (cardaurand()<psnooker){
 withinrange=STEP_DEMCMCZ_SNOOKER(&X[nn*PI.npars],Z,M,pars_new,PI,&gratio);
 nupdate=PI.npars;
 }else{
-int cridx=(int)(((double)random()/((double)RAND_MAX))*nCR);
+int cridx=(int)(((double)cardarand()/((double)CARDAMOM_RAND_MAX))*nCR);
 if (cridx>=nCR){cridx=nCR-1;}
 CR=CRvals[cridx];
 withinrange=STEP_DREAMZS_PARALLEL(&X[nn*PI.npars],Z,M,pars_new,PI,CR,&nupdate);
 }
 totalupdates=totalupdates+nupdate;
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 
 if (withinrange==1){
 wrlocal=wrlocal+1;