diff --git a/C/auxi_fun/cardamom_random.h b/C/auxi_fun/cardamom_random.h
new file mode 100644
index 00000000..05bb2773
--- /dev/null
+++ 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
index 92ada834..683ded3f 100644
--- 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
index 51499755..6ad1d693 100644
--- 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
index 01409565..dba3106c 100644
--- 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,7 +160,7 @@ 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);
@@ -167,7 +168,7 @@ for ( ;N.ITER<MCO.nOUT;N.ITER++){
 	}
 
 
-	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
index 17138472..df6f2283 100644
--- 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,7 +159,7 @@ 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);
@@ -166,7 +167,7 @@ for (N.ITER=0;N.ITER<MCO.nOUT;N.ITER++){
 	}
 
 	
-	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
index c0f4d379..4ea6d5b2 100644
--- 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,7 +160,7 @@ 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);
@@ -167,7 +168,7 @@ for ( ;N.ITER<MCO.nOUT;N.ITER++){
 	}
 
 
-	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
index e99a4a08..25a4d477 100644
--- 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
index 6a46f118..556700f7 100644
--- 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
index e8691634..52664d05 100644
--- 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
index bffeca67..70ab0239 100644
--- 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
index f65183b4..bd7a1aa3 100644
--- 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;
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE.c
index bdb7e43c..5665cce5 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE.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_HYBRID_AIDE.c"
@@ -43,7 +44,7 @@ N.ACCRATE=0;
 for (nn=0;nn<NC;nn++){
 for (n=0;n<PI.npars;n++){
 if (MCO.randparini==1 && PI.parfix[n]!=1){
-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;
@@ -71,7 +72,7 @@ for (nn=0;nn<NC;nn++){
 gratio=0;
 withinrange=STEP_HYBRID_AIDE(PARS,pars_new,PI,nn,NC,&gratio);
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 if (withinrange==1){
 wrlocal=wrlocal+1;
 P_new=MODEL_LIKELIHOOD(DATA,pars_new);
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE_DEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE_DEMCMC.c
index 36ded595..88240794 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE_DEMCMC.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/HYBRID_AIDE_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_HYBRID_AIDE.c"
@@ -50,7 +51,7 @@ printf("HYBRID_AIDE_DEMCMC: AIDE phase ends at iteration %d out of %d\n",switch_
 for (nn=0;nn<NC;nn++){
 for (n=0;n<PI.npars;n++){
 if (MCO.randparini==1 && PI.parfix[n]!=1){
-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;
@@ -81,12 +82,12 @@ withinrange=STEP_HYBRID_AIDE(PARS,pars_new,PI,nn,NC,&gratio);
 }else{
 /*Use CARDAMOM's standard DEMCMC gamma mixture before STEP_DEMCMC.*/
 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);
 withinrange=STEP_DEMCMC(PARS,pars_new,PI,nn,NC);
 gratio=0;
 }
 
-lr=log((double)random()/(double)RAND_MAX);
+lr=log(cardaurand());
 if (withinrange==1){
 wrlocal=wrlocal+1;
 P_new=MODEL_LIKELIHOOD(DATA,pars_new);
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/MHMCMC_119.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/MHMCMC_119.c
index 41727e90..6c857c89 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/MHMCMC_119.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/MHMCMC_119.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 "../../../math_fun/covariance.c"
 #include "NORMPARS.c"
@@ -104,7 +105,7 @@ if (MCO.fixedpars!=1){PI.parfix[n]=0;}
 /*BUG IS HERE!!!4-9-2013*/
 if (MCO.randparini==1 && PI.parfix[n]!=1){
 /*random parameter if PI.parini = -9999*/
-PI.parini[n]=nor2par((double)random()/(double)RAND_MAX,PI.parmin[n],PI.parmax[n]);}
+PI.parini[n]=nor2par(cardaurand(),PI.parmin[n],PI.parmax[n]);}
 /*printing parameter values*/
 printf("log10(p%d)=%2.4f ",n+1,log10(PI.parini[n]));
 if ((n+1) % 3==0){printf("\n");}
@@ -153,7 +154,7 @@ while (N.ITER < MCO.nOUT){
 	if (isnan(P)){printf("Warning: MLF generated NaN... treating as -Inf\n");P=log(0);break;}
 
 
-	if (P-P0>log((double)random()/(double)RAND_MAX)){
+	if (P-P0>log(cardaurand())){
 		/*storing accepted solution*/
 		for (n=0;n<PI.npars;n++){
 		PARS0[n]=PARS[n];
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_ADEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_ADEMCMC.c
index fd0cb954..5b31e8e8 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_ADEMCMC.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_ADEMCMC.c
@@ -15,7 +15,7 @@ double rn;
 /*STep 1. Find random sample other than current chain*/
 
 C1=C;
-while (C1==C){C1=ceil((double)random()*NC/((double)RAND_MAX))-1;};
+while (C1==C){C1=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;};
 
 /*Step 2. FInd current chain, reference chain, and step between them*/
 /*allocating memory for three normalized parameter vectors*/
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_AFDEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_AFDEMCMC.c
index ea887c62..3e87a86f 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_AFDEMCMC.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_AFDEMCMC.c
@@ -11,7 +11,7 @@ double STEP_AFDEMCMC(double *PARS, double *pars_new, PARAMETER_INFO PI, int C, i
     /* Step 1. Find random sample (reference chain) other than current chain */
     C1=C;
     while (C1==C){
-        C1 = ceil((double)random() * NC / ((double)RAND_MAX)) - 1;
+        C1 = ceil((double)cardarand() * NC / ((double)CARDAMOM_RAND_MAX)) - 1;
     }
 
     /* Step 2. Allocate memory for normalized parameter vectors */
@@ -28,14 +28,14 @@ double STEP_AFDEMCMC(double *PARS, double *pars_new, PARAMETER_INFO PI, int C, i
 
     /* Step 3. The Affine Invariant Step */
    // double a_gw = 2.0; // The Goodman-Weare tuning parameter
-   // double u = (double)random() / (double)RAND_MAX; // Uniform distribution between 0 and 1
+   // double u = (double)cardarand() / (double)CARDAMOM_RAND_MAX; // Uniform distribution between 0 and 1
  
 // 1. Draw a_gw using a shifted Log-Normal distribution: 1 + exp(N(0,1))
     // (This heavily favors values between 1.1 and 2.5, but occasionally spikes much higher)
     double a_gw = 1.0 + exp(randn()); 
 
     // 2. Draw 'u' for the z transformation
-    double u_z = (double)random() / (double)RAND_MAX; 
+    double u_z = (double)cardarand() / (double)CARDAMOM_RAND_MAX; 
     
     // 3. Transform into the 1/sqrt(z) distribution using our random a_gw
     double z = pow((u_z * (sqrt(a_gw) - 1.0/sqrt(a_gw)) + 1.0/sqrt(a_gw)), 2.0);
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMC.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMC.c
index 3bef5d45..5d053d48 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMC.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMC.c
@@ -28,8 +28,8 @@ int withinlim=1;
 	/*first: pick c1 and c2, make sure that (a) they are between 1 and NC, and (b) they are not C*/
 	c1=C;c2=C;
 	while (c1==C | c2==C | c1==c2){
-	c1=ceil((double)random()*NC/((double)RAND_MAX))-1;
-	c2=ceil((double)random()*NC/((double)RAND_MAX))-1;}
+	c1=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;
+	c2=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;}
 mstep=0;
 	/*SAMPLING PARAMETERS*/
 	for (n=0;n<PI.npars;n++){
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMCZS.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMCZS.c
index 8b78c4fc..e1c58012 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMCZS.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DEMCMCZS.c
@@ -19,12 +19,12 @@ for (n=0;n<PI.npars;n++){npar[n]=par2nor(xi[n],PI.parmin[n],PI.parmax[n]);}
 
 /*picking two distinct rows of the archive*/
 while (r1==r2){
-r1=ceil((double)random()*M/((double)RAND_MAX))-1;
-r2=ceil((double)random()*M/((double)RAND_MAX))-1;}
+r1=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;
+r2=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;}
 
 /*gamma: standard DE-MC scaling, with an occasional large jump (10% of updates)*/
 double gamma_de=2.38/sqrt(2.0*(double)PI.npars);
-if ((double)random()/(double)RAND_MAX<0.1){gamma_de=1;}
+if (cardaurand()<0.1){gamma_de=1;}
 
 int withinlim=1;
 for (n=0;n<PI.npars;n++){
@@ -60,9 +60,9 @@ for (n=0;n<PI.npars;n++){npar[n]=par2nor(xi[n],PI.parmin[n],PI.parmax[n]);}
 
 /*picking three distinct rows: z (the snooker reference line partner) and r1,r2 (for the projection)*/
 while (rz==r1 || rz==r2 || r1==r2){
-rz=ceil((double)random()*M/((double)RAND_MAX))-1;
-r1=ceil((double)random()*M/((double)RAND_MAX))-1;
-r2=ceil((double)random()*M/((double)RAND_MAX))-1;}
+rz=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;
+r1=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;
+r2=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;}
 
 double dnorm2=0;
 for (n=0;n<PI.npars;n++){
@@ -86,7 +86,7 @@ s1=s1+(par2nor(Z[r1*PI.npars+n],PI.parmin[n],PI.parmax[n])-nz[n])*d[n]/dnorm;
 s2=s2+(par2nor(Z[r2*PI.npars+n],PI.parmin[n],PI.parmax[n])-nz[n])*d[n]/dnorm;}
 
 /*gamma_s ~ U[1.2,2.2], the randomised snooker step size ter Braak & Vrugt (2008) use as default*/
-double gammas=1.2+((double)random()/(double)RAND_MAX)*1.0;
+double gammas=1.2+(cardaurand())*1.0;
 
 int withinlim=1;
 double xstarnorm2=0;
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DREAMZS.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DREAMZS.c
index 8dde51ce..7b5890f4 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DREAMZS.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_DREAMZS.c
@@ -24,23 +24,23 @@ for (n=0;n<PI.npars;n++){npar[n]=par2nor(xi[n],PI.parmin[n],PI.parmax[n]);}
 
 /*picking two distinct rows of the archive*/
 while (r1==r2){
-r1=ceil((double)random()*M/((double)RAND_MAX))-1;
-r2=ceil((double)random()*M/((double)RAND_MAX))-1;}
+r1=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;
+r2=ceil((double)cardarand()*M/((double)CARDAMOM_RAND_MAX))-1;}
 
 /*Force at least one updated dimension so low CR values cannot produce a null move.*/
-int force_dim=(int)(((double)random()/((double)RAND_MAX))*PI.npars);
+int force_dim=(int)(((double)cardarand()/((double)CARDAMOM_RAND_MAX))*PI.npars);
 if (force_dim>=PI.npars){force_dim=PI.npars-1;}
 
 int nupdate=0;
 for (n=0;n<PI.npars;n++){
-if (n==force_dim || (double)random()/(double)RAND_MAX<CR){
+if (n==force_dim || cardaurand()<CR){
 update[n]=1;
 nupdate=nupdate+1;
 }}
 
 /*DREAM uses the number of updated dimensions for the DE scale.*/
 double gamma_de=2.38/sqrt(2.0*(double)nupdate);
-if ((double)random()/(double)RAND_MAX<0.1){gamma_de=1;}
+if (cardaurand()<0.1){gamma_de=1;}
 
 int withinlim=1;
 for (n=0;n<PI.npars;n++){
diff --git a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_HYBRID_AIDE.c b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_HYBRID_AIDE.c
index e08ef144..f0d5c246 100644
--- a/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_HYBRID_AIDE.c
+++ b/C/mcmc_fun/MHMCMC/MCMC_FUN/STEP_HYBRID_AIDE.c
@@ -11,9 +11,9 @@ int STEP_HYBRID_AIDE(double *PARS, double *pars_new, PARAMETER_INFO PI, int C, i
 int n,j1=C,j2=C,j3=C;
 double rn;
 
-while (j1==C){j1=ceil((double)random()*NC/((double)RAND_MAX))-1;}
-while (j2==C || j2==j1){j2=ceil((double)random()*NC/((double)RAND_MAX))-1;}
-while (j3==C || j3==j1 || j3==j2){j3=ceil((double)random()*NC/((double)RAND_MAX))-1;}
+while (j1==C){j1=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;}
+while (j2==C || j2==j1){j2=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;}
+while (j3==C || j3==j1 || j3==j2){j3=ceil((double)cardarand()*NC/((double)CARDAMOM_RAND_MAX))-1;}
 
 double *xcur=calloc(PI.npars,sizeof(double));
 double *x1=calloc(PI.npars,sizeof(double));
@@ -30,13 +30,13 @@ x3[n]=par2nor(PARS[j3*PI.npars+n],PI.parmin[n],PI.parmax[n]);
 }
 
 double gamma_de=1 - (1 - 2.38 / sqrt(2.0 * PI.npars) * 0.1) *
-	(double)(((double)random() / (double)RAND_MAX) < 0.9);
+	(double)(((double)cardarand() / (double)CARDAMOM_RAND_MAX) < 0.9);
 
 double a_gw=1.0+exp(randn());
-double u_z=(double)random()/(double)RAND_MAX;
+double u_z=cardaurand();
 double z=pow((u_z*(sqrt(a_gw)-1.0/sqrt(a_gw))+1.0/sqrt(a_gw)),2.0);
 
-int de_first=((double)random()/(double)RAND_MAX>0.5);
+int de_first=(cardaurand()>0.5);
 int withinlim=1;
 
 if (de_first){
diff --git a/C/projects/CARDAMOM_MDF/CARDAMOM_MDF.c b/C/projects/CARDAMOM_MDF/CARDAMOM_MDF.c
index d6dbd736..46841931 100644
--- a/C/projects/CARDAMOM_MDF/CARDAMOM_MDF.c
+++ b/C/projects/CARDAMOM_MDF/CARDAMOM_MDF.c
@@ -1,4 +1,5 @@
 #include <assert.h>
+#include "../../auxi_fun/cardamom_random.h"
 #include "../../auxi_fun/oksofar.c"
 #include "../../auxi_fun/okcheck.c"
 #include "../../auxi_fun/seedrandomnumber.c"