Random number generator: replace a bunch of usage of the global random number

generator with a local one. It's not thread safe and will not give repeatable
results, so in most cases it should not be used.

Also fixes #34992 where the noise texture of a displacement modifier was not
properly random in opengl animation render, because the seed got reset to a
fixed value by an unrelated function while for final render it changed each
frame.

31 files changed, 247 insertions, 158 deletions

diff --git a/source/blender/blenkernel/BKE_boids.h b/source/blender/blenkernel/BKE_boids.h
index bb724c6632e..582cd0cef8d 100644
--- a/source/blender/blenkernel/BKE_boids.h
+++ b/source/blender/blenkernel/BKE_boids.h
@@ -36,6 +36,8 @@
 
 #include "DNA_boid_types.h"
 
+struct RNG;
+
 typedef struct BoidBrainData {
 	struct ParticleSimulationData *sim;
 	struct ParticleSettings *part;
@@ -47,6 +49,8 @@ typedef struct BoidBrainData {
 	float goal_co[3];
 	float goal_nor[3];
 	float goal_priority;
+
+	struct RNG *rng;
 } BoidBrainData;
 
 void boids_precalc_rules(struct ParticleSettings *part, float cfra);
diff --git a/source/blender/blenkernel/BKE_particle.h b/source/blender/blenkernel/BKE_particle.h
index 2b753cba098..846fe68c2b3 100644
--- a/source/blender/blenkernel/BKE_particle.h
+++ b/source/blender/blenkernel/BKE_particle.h
@@ -61,6 +61,7 @@ struct SurfaceModifierData;
 struct BVHTreeRay;
 struct BVHTreeRayHit; 
 struct EdgeHash;
+struct RNG;
 
 #define PARTICLE_P              ParticleData * pa; int p
 #define LOOP_PARTICLES  for (p = 0, pa = psys->particles; p < psys->totpart; p++, pa++)
@@ -82,6 +83,7 @@ typedef struct ParticleSimulationData {
 	struct ParticleSystem *psys;
 	struct ParticleSystemModifierData *psmd;
 	struct ListBase *colliders;
+	struct RNG *rng;
 	/* Courant number. This is used to implement an adaptive time step. Only the
 	 * maximum value per time step is important. Only sph_integrate makes use of
 	 * this at the moment. Other solvers could, too. */
diff --git a/source/blender/blenkernel/intern/anim.c b/source/blender/blenkernel/intern/anim.c
index 6fa33e8d3c5..1f54da402ab 100644
--- a/source/blender/blenkernel/intern/anim.c
+++ b/source/blender/blenkernel/intern/anim.c
@@ -1263,6 +1263,7 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 	ChildParticle *cpa = NULL;
 	ParticleKey state;
 	ParticleCacheKey *cache;
+	RNG *rng = NULL;
 	float ctime, pa_time, scale = 1.0f;
 	float tmat[4][4], mat[4][4], pamat[4][4], vec[3], size = 0.0;
 	float (*obmat)[4], (*oldobmat)[4];
@@ -1293,14 +1294,9 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 	totpart = psys->totpart;
 	totchild = psys->totchild;
 
-	BLI_srandom(31415926 + psys->seed);
-
 	if ((psys->renderdata || part->draw_as == PART_DRAW_REND) && ELEM(part->ren_as, PART_DRAW_OB, PART_DRAW_GR)) {
 		ParticleSimulationData sim = {NULL};
-		sim.scene = scene;
-		sim.ob = par;
-		sim.psys = psys;
-		sim.psmd = psys_get_modifier(par, psys);
+
 		/* make sure emitter imat is in global coordinates instead of render view coordinates */
 		invert_m4_m4(par->imat, par->obmat);
 
@@ -1332,6 +1328,12 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 
 		psys_check_group_weights(part);
 
+		sim.scene = scene;
+		sim.ob = par;
+		sim.psys = psys;
+		sim.psmd = psys_get_modifier(par, psys);
+		sim.rng = BLI_rng_new(0);
+
 		psys->lattice = psys_get_lattice(&sim);
 
 		/* gather list of objects or single object */
@@ -1378,6 +1380,8 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 			obcopy = *ob;
 		}
 
+		rng = BLI_rng_new_srandom(31415926 + psys->seed);
+
 		if (totchild == 0 || part->draw & PART_DRAW_PARENT)
 			a = 0;
 		else
@@ -1417,7 +1421,7 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 
 				/* for groups, pick the object based on settings */
 				if (part->draw & PART_DRAW_RAND_GR)
-					b = BLI_rand() % totgroup;
+					b = BLI_rng_get_int(rng) % totgroup;
 				else
 					b = a % totgroup;
 
@@ -1561,6 +1565,8 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 		}
 		else
 			*ob = obcopy;
+
+		BLI_rng_free(sim.rng);
 	}
 
 	/* clean up */
@@ -1573,6 +1579,9 @@ static void new_particle_duplilist(ListBase *lb, ID *id, Scene *scene, Object *p
 		end_latt_deform(psys->lattice);
 		psys->lattice = NULL;
 	}
+
+	if (rng)
+		BLI_rng_free(rng);
 }
 
 static Object *find_family_object(Object **obar, char *family, char ch)
diff --git a/source/blender/blenkernel/intern/boids.c b/source/blender/blenkernel/intern/boids.c
index b6f1b88c912..dfffb7c795e 100644
--- a/source/blender/blenkernel/intern/boids.c
+++ b/source/blender/blenkernel/intern/boids.c
@@ -234,9 +234,9 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
 			if (dot_v3v3(col.pce.nor, pa->prev_state.ave) < -0.99f) {
 				/* don't know why, but uneven range [0.0, 1.0] */
 				/* works much better than even [-1.0, 1.0] */
-				bbd->wanted_co[0] = BLI_frand();
-				bbd->wanted_co[1] = BLI_frand();
-				bbd->wanted_co[2] = BLI_frand();
+				bbd->wanted_co[0] = BLI_rng_get_float(bbd->rng);
+				bbd->wanted_co[1] = BLI_rng_get_float(bbd->rng);
+				bbd->wanted_co[2] = BLI_rng_get_float(bbd->rng);
 			}
 			else {
 				copy_v3_v3(bbd->wanted_co, col.pce.nor);
@@ -558,9 +558,9 @@ static int rule_average_speed(BoidRule *rule, BoidBrainData *bbd, BoidValues *va
 
 	if (asbr->wander > 0.0f) {
 		/* abuse pa->r_ave for wandering */
-		bpa->wander[0] += asbr->wander * (-1.0f + 2.0f * BLI_frand());
-		bpa->wander[1] += asbr->wander * (-1.0f + 2.0f * BLI_frand());
-		bpa->wander[2] += asbr->wander * (-1.0f + 2.0f * BLI_frand());
+		bpa->wander[0] += asbr->wander * (-1.0f + 2.0f * BLI_rng_get_float(bbd->rng));
+		bpa->wander[1] += asbr->wander * (-1.0f + 2.0f * BLI_rng_get_float(bbd->rng));
+		bpa->wander[2] += asbr->wander * (-1.0f + 2.0f * BLI_rng_get_float(bbd->rng));
 
 		normalize_v3(bpa->wander);
 
@@ -586,9 +586,9 @@ static int rule_average_speed(BoidRule *rule, BoidBrainData *bbd, BoidValues *va
 
 		/* may happen at birth */
 		if (dot_v2v2(bbd->wanted_co, bbd->wanted_co)==0.0f) {
-			bbd->wanted_co[0] = 2.0f*(0.5f - BLI_frand());
-			bbd->wanted_co[1] = 2.0f*(0.5f - BLI_frand());
-			bbd->wanted_co[2] = 2.0f*(0.5f - BLI_frand());
+			bbd->wanted_co[0] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
+			bbd->wanted_co[1] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
+			bbd->wanted_co[2] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
 		}
 		
 		/* leveling */
@@ -663,7 +663,7 @@ static int rule_fight(BoidRule *rule, BoidBrainData *bbd, BoidValues *val, Parti
 
 		/* attack if in range */
 		if (closest_dist <= bbd->part->boids->range + pa->size + enemy_pa->size) {
-			float damage = BLI_frand();
+			float damage = BLI_rng_get_float(bbd->rng);
 			float enemy_dir[3];
 
 			normalize_v3_v3(enemy_dir, bbd->wanted_co);
@@ -1164,9 +1164,9 @@ void boid_body(BoidBrainData *bbd, ParticleData *pa)
 			/* choose random direction to turn if wanted velocity */
 			/* is directly behind regardless of z-coordinate */
 			if (dot_v2v2(old_dir2, wanted_dir2) < -0.99f) {
-				wanted_dir[0] = 2.0f*(0.5f - BLI_frand());
-				wanted_dir[1] = 2.0f*(0.5f - BLI_frand());
-				wanted_dir[2] = 2.0f*(0.5f - BLI_frand());
+				wanted_dir[0] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
+				wanted_dir[1] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
+				wanted_dir[2] = 2.0f*(0.5f - BLI_rng_get_float(bbd->rng));
 				normalize_v3(wanted_dir);
 			}
 
diff --git a/source/blender/blenkernel/intern/effect.c b/source/blender/blenkernel/intern/effect.c
index 9e8693e957e..3fbbc484169 100644
--- a/source/blender/blenkernel/intern/effect.c
+++ b/source/blender/blenkernel/intern/effect.c
@@ -619,6 +619,7 @@ int get_effector_data(EffectorCache *eff, EffectorData *efd, EffectedPoint *poin
 			sim.scene= eff->scene;
 			sim.ob= eff->ob;
 			sim.psys= eff->psys;
+			sim.rng= NULL;
 
 			/* TODO: time from actual previous calculated frame (step might not be 1) */
 			state.time = cfra - 1.0f;
diff --git a/source/blender/blenkernel/intern/ocean.c b/source/blender/blenkernel/intern/ocean.c
index 695ac7da792..94dbe7b87c2 100644
--- a/source/blender/blenkernel/intern/ocean.c
+++ b/source/blender/blenkernel/intern/ocean.c
@@ -144,12 +144,12 @@ typedef struct Ocean {
 
 
 
-static float nextfr(float min, float max)
+static float nextfr(RNG *rng, float min, float max)
 {
-	return BLI_frand() * (min - max) + max;
+	return BLI_rng_get_float(rng) * (min - max) + max;
 }
 
-static float gaussRand(void)
+static float gaussRand(RNG *rng)
 {
 	/* Note: to avoid numerical problems with very small numbers, we make these variables singe-precision floats,
 	 * but later we call the double-precision log() and sqrt() functions instead of logf() and sqrtf().
@@ -765,6 +765,7 @@ void BKE_init_ocean(struct Ocean *o, int M, int N, float Lx, float Lz, float V,
                     float alignment, float depth, float time, short do_height_field, short do_chop, short do_normals,
                     short do_jacobian, int seed)
 {
+	RNG *rng;
 	int i, j, ii;
 
 	BLI_rw_mutex_lock(&o->oceanmutex, THREAD_LOCK_WRITE);
@@ -825,12 +826,12 @@ void BKE_init_ocean(struct Ocean *o, int M, int N, float Lx, float Lz, float V,
 			o->_k[i * (1 + o->_N / 2) + j] = sqrt(o->_kx[i] * o->_kx[i] + o->_kz[j] * o->_kz[j]);
 
 	/*srand(seed);*/
-	BLI_srand(seed);
+	rng = BLI_rng_new(seed);
 
 	for (i = 0; i < o->_M; ++i) {
 		for (j = 0; j < o->_N; ++j) {
-			float r1 = gaussRand();
-			float r2 = gaussRand();
+			float r1 = gaussRand(rng);
+			float r2 = gaussRand(rng);
 
 			fftw_complex r1r2;
 			init_complex(r1r2, r1, r2);
@@ -890,6 +891,7 @@ void BKE_init_ocean(struct Ocean *o, int M, int N, float Lx, float Lz, float V,
 
 	set_height_normalize_factor(o);
 
+	BLI_rng_free(rng);
 }
 
 void BKE_free_ocean_data(struct Ocean *oc)
@@ -1188,13 +1190,14 @@ void BKE_bake_ocean(struct Ocean *o, struct OceanCache *och, void (*update_cb)(v
 	int res_x = och->resolution_x;
 	int res_y = och->resolution_y;
 	char string[FILE_MAX];
+	//RNG *rng;
 
 	if (!o) return;
 
 	if (o->_do_jacobian) prev_foam = MEM_callocN(res_x * res_y * sizeof(float), "previous frame foam bake data");
 	else prev_foam = NULL;
 
-	BLI_srand(0);
+	//rng = BLI_rng_new(0);
 
 	/* setup image format */
 	imf.imtype = R_IMF_IMTYPE_OPENEXR;
@@ -1232,7 +1235,7 @@ void BKE_bake_ocean(struct Ocean *o, struct OceanCache *och, void (*update_cb)(v
 						pr = prev_foam[res_x * y + x];
 					}
 
-					/* r = BLI_frand(); */ /* UNUSED */ /* randomly reduce foam */
+					/* r = BLI_rng_get_float(rng); */ /* UNUSED */ /* randomly reduce foam */
 
 					/* pr = pr * och->foam_fade; */		/* overall fade */
 
@@ -1311,10 +1314,12 @@ void BKE_bake_ocean(struct Ocean *o, struct OceanCache *och, void (*update_cb)(v
 
 		if (cancel) {
 			if (prev_foam) MEM_freeN(prev_foam);
+			//BLI_rng_free(rng);
 			return;
 		}
 	}
 
+	//BLI_rng_free(rng);
 	if (prev_foam) MEM_freeN(prev_foam);
 	och->baked = 1;
 }
diff --git a/source/blender/blenkernel/intern/particle.c b/source/blender/blenkernel/intern/particle.c
index 4e78ab52499..b0ef6715382 100644
--- a/source/blender/blenkernel/intern/particle.c
+++ b/source/blender/blenkernel/intern/particle.c
@@ -264,11 +264,12 @@ static void psys_create_frand(ParticleSystem *psys)
 {
 	int i;
 	float *rand = psys->frand = MEM_callocN(PSYS_FRAND_COUNT * sizeof(float), "particle randoms");
-
-	BLI_srandom(psys->seed);
+	RNG *rng = BLI_rng_new_srandom(psys->seed);
 
 	for (i = 0; i < 1024; i++, rand++)
-		*rand = BLI_frand();
+		*rand = BLI_rng_get_float(rng);
+
+	BLI_rng_free(rng);
 }
 int psys_check_enabled(Object *ob, ParticleSystem *psys)
 {
@@ -3302,8 +3303,11 @@ void psys_cache_edit_paths(Scene *scene, Object *ob, PTCacheEdit *edit, float cf
 		sim.ob = ob;
 		sim.psys = psys;
 		sim.psmd = psys_get_modifier(ob, psys);
+		sim.rng = BLI_rng_new(0);
 
 		psys_cache_child_paths(&sim, cfra, 1);
+
+		BLI_rng_free(sim.rng);
 	}
 
 	/* clear recalc flag if set here */
@@ -4660,6 +4664,7 @@ void psys_apply_hair_lattice(Scene *scene, Object *ob, ParticleSystem *psys)
 	sim.ob = ob;
 	sim.psys = psys;
 	sim.psmd = psys_get_modifier(ob, psys);
+	sim.rng = BLI_rng_new(0);
 
 	psys->lattice = psys_get_lattice(&sim);
 
@@ -4687,4 +4692,6 @@ void psys_apply_hair_lattice(Scene *scene, Object *ob, ParticleSystem *psys)
 		/* protect the applied shape */
 		psys->flag |= PSYS_EDITED;
 	}
+
+	BLI_rng_free(sim.rng);
 }
diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c
index 5124cb1da26..f131e7297ed 100644
--- a/source/blender/blenkernel/intern/particle_system.c
+++ b/source/blender/blenkernel/intern/particle_system.c
@@ -423,7 +423,7 @@ void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
 	}
 }
 
-static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, ParticleSystem *psys)
+static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, ParticleSystem *psys, RNG *rng)
 {
 	ChildParticle *cpa = NULL;
 	int i, p;
@@ -440,9 +440,9 @@ static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *fi
 					
 			/* create even spherical distribution inside unit sphere */
 			while (length>=1.0f) {
-				cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
-				cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
-				cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
+				cpa->fuv[0]=2.0f*BLI_rng_get_float(rng)-1.0f;
+				cpa->fuv[1]=2.0f*BLI_rng_get_float(rng)-1.0f;
+				cpa->fuv[2]=2.0f*BLI_rng_get_float(rng)-1.0f;
 				length=len_v3(cpa->fuv);
 			}
 
@@ -872,7 +872,7 @@ static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, Ch
 						pa->foffset *= ctx->jit[p % (2 * ctx->jitlevel)];
 						break;
 					case PART_DISTR_RAND:
-						pa->foffset *= BLI_frand();
+						pa->foffset *= BLI_rng_get_float(ctx->sim.rng);
 						break;
 				}
 			}
@@ -1065,15 +1065,15 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
 	if (from == PART_FROM_CHILD) {
 		/* Simple children */
 		if (part->childtype != PART_CHILD_FACES) {
-			BLI_srandom(31415926 + psys->seed + psys->child_seed);
-			distribute_simple_children(scene, ob, finaldm, psys);
+			BLI_rng_srandom(ctx->sim.rng, 31415926 + psys->seed + psys->child_seed);
+			distribute_simple_children(scene, ob, finaldm, psys, ctx->sim.rng);
 			return 0;
 		}
 	}
 	else {
 		/* Grid distribution */
 		if (part->distr==PART_DISTR_GRID && from != PART_FROM_VERT) {
-			BLI_srandom(31415926 + psys->seed);
+			BLI_rng_srandom(ctx->sim.rng, 31415926 + psys->seed);
 			dm= CDDM_from_mesh((Mesh*)ob->data, ob);
 			DM_ensure_tessface(dm);
 			distribute_grid(dm,psys);
@@ -1085,7 +1085,7 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
 	/* Create trees and original coordinates if needed */
 	if (from == PART_FROM_CHILD) {
 		distr=PART_DISTR_RAND;
-		BLI_srandom(31415926 + psys->seed + psys->child_seed);
+		BLI_rng_srandom(ctx->sim.rng, 31415926 + psys->seed + psys->child_seed);
 		dm= finaldm;
 
 		/* BMESH ONLY */
@@ -1108,7 +1108,7 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
 	}
 	else {
 		distr = part->distr;
-		BLI_srandom(31415926 + psys->seed);
+		BLI_rng_srandom(ctx->sim.rng, 31415926 + psys->seed);
 		
 		dm= CDDM_from_mesh((Mesh*)ob->data, ob);
 
@@ -1264,7 +1264,7 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
 
 		for (p=0; p<totpart; p++) {
 			/* In theory element_sum[totelem] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
-			pos= BLI_frand() * element_sum[totelem];
+			pos= BLI_rng_get_float(ctx->sim.rng) * element_sum[totelem];
 			particle_element[p] = distribute_binary_search(element_sum, totelem, pos);
 			particle_element[p] = MIN2(totelem-1, particle_element[p]);
 			jitter_offset[particle_element[p]] = pos;
@@ -2916,9 +2916,9 @@ static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, flo
 
 	/* brownian force */
 	if (part->brownfac != 0.0f) {
-		force[0] += (BLI_frand()-0.5f) * part->brownfac;
-		force[1] += (BLI_frand()-0.5f) * part->brownfac;
-		force[2] += (BLI_frand()-0.5f) * part->brownfac;
+		force[0] += (BLI_rng_get_float(sim->rng)-0.5f) * part->brownfac;
+		force[1] += (BLI_rng_get_float(sim->rng)-0.5f) * part->brownfac;
+		force[2] += (BLI_rng_get_float(sim->rng)-0.5f) * part->brownfac;
 	}
 
 	if (part->flag & PART_ROT_DYN && epoint.ave)
@@ -3497,7 +3497,7 @@ static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRay
 
 	return hit->index >= 0;
 }
-static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, int kill, int dynamic_rotation)
+static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, int kill, int dynamic_rotation, RNG *rng)
 {
 	ParticleCollisionElement *pce = &col->pce;
 	PartDeflect *pd = col->hit->pd;
@@ -3506,7 +3506,7 @@ static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeR
 	float f = col->f + x * (1.0f - col->f);				/* time factor of collision between timestep */
 	float dt1 = (f - col->f) * col->total_time;			/* time since previous collision (in seconds) */
 	float dt2 = (1.0f - f) * col->total_time;			/* time left after collision (in seconds) */
-	int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
+	int through = (BLI_rng_get_float(rng) < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
 
 	/* calculate exact collision location */
 	interp_v3_v3v3(co, col->co1, col->co2, x);
@@ -3531,8 +3531,8 @@ static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeR
 		float v0_tan[3];/* tangential component of v0 */
 		float vc_tan[3];/* tangential component of collision surface velocity */
 		float v0_dot, vc_dot;
-		float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
-		float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
+		float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_rng_get_float(rng) - 0.5f);
+		float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_rng_get_float(rng) - 0.5f);
 		float distance, nor[3], dot;
 
 		CLAMP(damp,0.0f, 1.0f);
@@ -3740,7 +3740,7 @@ static void collision_check(ParticleSimulationData *sim, int p, float dfra, floa
 
 			if (collision_count == COLLISION_MAX_COLLISIONS)
 				collision_fail(pa, &col);
-			else if (collision_response(pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN)==0)
+			else if (collision_response(pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN, sim->rng)==0)
 				return;
 		}
 		else
@@ -4100,6 +4100,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 {
 	ParticleSystem *psys = sim->psys;
 	ParticleSettings *part=psys->part;
+	RNG *rng;
 	BoidBrainData bbd;
 	ParticleTexture ptex;
 	PARTICLE_P;
@@ -4126,7 +4127,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 		return;
 	}
 
-	BLI_srandom(31415926 + (int)cfra + psys->seed);
+	rng = BLI_rng_new_srandom(31415926 + (int)cfra + psys->seed);
 
 	psys_update_effectors(sim);
 
@@ -4143,6 +4144,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 			bbd.cfra = cfra;
 			bbd.dfra = dfra;
 			bbd.timestep = timestep;
+			bbd.rng = rng;
 
 			psys_update_particle_tree(psys, cfra);
 
@@ -4326,6 +4328,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 	}
 
 	free_collider_cache(&sim->colliders);
+	BLI_rng_free(rng);
 }
 static void update_children(ParticleSimulationData *sim)
 {
@@ -4822,6 +4825,8 @@ void particle_system_update(Scene *scene, Object *ob, ParticleSystem *psys)
 	if (!sim.psmd->dm)
 		return;
 
+	sim.rng = BLI_rng_new(0);
+
 	if (part->from != PART_FROM_VERT) {
 		DM_ensure_tessface(sim.psmd->dm);
 	}
@@ -4962,5 +4967,7 @@ void particle_system_update(Scene *scene, Object *ob, ParticleSystem *psys)
 	/* save matrix for duplicators, at rendertime the actual dupliobject's matrix is used so don't update! */
 	if (psys->renderdata==0)
 		invert_m4_m4(psys->imat, ob->obmat);
+	
+	BLI_rng_free(sim.rng);
 }
 
diff --git a/source/blender/blenkernel/intern/smoke.c b/source/blender/blenkernel/intern/smoke.c
index 6b8f0348e16..ba6e8fe0a3a 100644
--- a/source/blender/blenkernel/intern/smoke.c
+++ b/source/blender/blenkernel/intern/smoke.c
@@ -969,6 +969,7 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 		sim.scene = scene;
 		sim.ob = flow_ob;
 		sim.psys = psys;
+		sim.rng = BLI_rng_new(sim.rng);
 
 		if (psys->part->type == PART_HAIR)
 		{
@@ -1055,6 +1056,8 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 			MEM_freeN(particle_pos);
 		if (particle_vel)
 			MEM_freeN(particle_vel);
+
+		BLI_rng_free(sim.rng);
 	}
 }
 
diff --git a/source/blender/blenlib/BLI_rand.h b/source/blender/blenlib/BLI_rand.h
index 3bc9a3bdb32..378beff3aa0 100644
--- a/source/blender/blenlib/BLI_rand.h
+++ b/source/blender/blenlib/BLI_rand.h
@@ -33,14 +33,15 @@
  *  \brief Random number functions.
  */
 
-/** RNG is just an abstract random number generator
- * type that avoids using globals, otherwise identical
- * to BLI_rand functions below.
+/* RNG is an abstract random number generator type that avoids using globals.
+ * Always use this instead of the global RNG unless you have a good reason,
+ * the global RNG is not thread safe and will not give repeatable results.
  */
 struct RNG;
 typedef struct RNG RNG;
 
 struct RNG *BLI_rng_new(unsigned int seed);
+struct RNG *BLI_rng_new_srandom(unsigned int seed);
 void        BLI_rng_free(struct RNG *rng);
 
 void        BLI_rng_seed(struct RNG *rng, unsigned int seed);
@@ -53,27 +54,17 @@ void        BLI_rng_shuffle_array(struct RNG *rng, void *data, int elemSize, int
 /** Note that skipping is as slow as generating n numbers! */
 void        BLI_rng_skip(struct RNG *rng, int n);
 
-/** Seed the random number generator */
-void    BLI_srand(unsigned int seed);
-
-/** Better seed for the random number generator, using noise.c hash[] */
+/** Seed for the random number generator, using noise.c hash[] */
 void    BLI_srandom(unsigned int seed);
 
 /** Return a pseudo-random number N where 0<=N<(2^31) */
 int     BLI_rand(void);
 
-/** Return a pseudo-random number N where 0.0<=N<1.0 */
-double  BLI_drand(void);
-
 /** Return a pseudo-random number N where 0.0f<=N<1.0f */
 float   BLI_frand(void);
 
-/** Fills a block of memory starting at \a addr
- * and extending \a len bytes with pseudo-random
- * contents. This routine does not use nor modify
- * the state of the BLI random number generator.
- */
-void    BLI_fillrand(void *addr, int len);
+/** Return a pseudo-random (hash) float from an integer value */
+float	BLI_hash_frand(unsigned int seed);
 
 /** Shuffle an array randomly using the given seed.
  * contents. This routine does not use nor modify
diff --git a/source/blender/blenlib/intern/jitter.c b/source/blender/blenlib/intern/jitter.c
index 3fe0ef158df..7141bedcf05 100644
--- a/source/blender/blenlib/intern/jitter.c
+++ b/source/blender/blenlib/intern/jitter.c
@@ -139,6 +139,7 @@ void BLI_jitterate2(float *jit1, float *jit2, int num, float rad2)
 void BLI_jitter_init(float *jitarr, int num)
 {
 	float *jit2, x, rad1, rad2, rad3;
+	RNG *rng;
 	int i;
 
 	if (num == 0) return;
@@ -148,15 +149,18 @@ void BLI_jitter_init(float *jitarr, int num)
 	rad2 = 1.0f / ((float)num);
 	rad3 = sqrtf((float)num) / ((float)num);
 
-	BLI_srand(31415926 + num);
+	rng = BLI_rng_new(31415926 + num);
+
 	x = 0;
 	for (i = 0; i < 2 * num; i += 2) {
-		jitarr[i] = x + rad1 * (float)(0.5 - BLI_drand());
-		jitarr[i + 1] = ((float)i / 2) / num + rad1 * (float)(0.5 - BLI_drand());
+		jitarr[i] = x + rad1 * (float)(0.5 - BLI_rng_get_double(rng));
+		jitarr[i + 1] = ((float)i / 2) / num + rad1 * (float)(0.5 - BLI_rng_get_double(rng));
 		x += rad3;
 		x -= floorf(x);
 	}
 
+	BLI_rng_free(rng);
+
 	for (i = 0; i < 24; i++) {
 		BLI_jitterate1(jitarr, jit2, num, rad1);
 		BLI_jitterate1(jitarr, jit2, num, rad1);
diff --git a/source/blender/blenlib/intern/rand.c b/source/blender/blenlib/intern/rand.c
index 76d17f34b5e..743d910e418 100644
--- a/source/blender/blenlib/intern/rand.c
+++ b/source/blender/blenlib/intern/rand.c
@@ -73,6 +73,15 @@ RNG *BLI_rng_new(unsigned int seed)
 	return rng;
 }
 
+RNG *BLI_rng_new_srandom(unsigned int seed)
+{
+	RNG *rng = MEM_mallocN(sizeof(*rng), "rng");
+
+	BLI_rng_srandom(rng, seed);
+
+	return rng;
+}
+
 void BLI_rng_free(RNG *rng)
 {
 	MEM_freeN(rng);
@@ -145,13 +154,8 @@ void BLI_rng_skip(RNG *rng, int n)
 
 /***/
 
-static RNG theBLI_rng = {0};
-
-/* note, this one creates periodical patterns */
-void BLI_srand(unsigned int seed)
-{
-	BLI_rng_seed(&theBLI_rng, seed);
-}
+/* initialize with some non-zero seed */
+static RNG theBLI_rng = {611330372042337130};
 
 /* using hash table to create better seed */
 void BLI_srandom(unsigned int seed)
@@ -164,23 +168,27 @@ int BLI_rand(void)
 	return BLI_rng_get_int(&theBLI_rng);
 }
 
-double BLI_drand(void)
-{
-	return BLI_rng_get_double(&theBLI_rng);
-}
-
 float BLI_frand(void)
 {
 	return BLI_rng_get_float(&theBLI_rng);
 }
 
-void BLI_fillrand(void *addr, int len)
+float BLI_hash_frand(unsigned int seed)
 {
-	RNG rng;
-	unsigned char *p = addr;
+	r_uint64 X;
+
+	seed = seed + hash[seed & 255];
+	X = (((r_uint64) seed) << 16) | LOWSEED;
+	seed = (int)(((MULTIPLIER * X + ADDEND) & MASK) >> 17);
+
+	seed = seed + hash[seed & 255];
+	X = (((r_uint64) seed) << 16) | LOWSEED;
+	X = (int)(((MULTIPLIER * X + ADDEND) & MASK) >> 17);
+
+	seed = seed + hash[seed & 255];
+	X = (((r_uint64) seed) << 16) | LOWSEED;
 
-	BLI_rng_seed(&rng, (unsigned int) (PIL_check_seconds_timer() * 0x7FFFFFFF));
-	while (len--) *p++ = BLI_rng_get_int(&rng) & 0xFF;
+	return (int)(((MULTIPLIER * X + ADDEND) & MASK) >> 17);
 }
 
 void BLI_array_randomize(void *data, int elemSize, int numElems, unsigned int seed)
diff --git a/source/blender/bmesh/operators/bmo_edgenet.c b/source/blender/bmesh/operators/bmo_edgenet.c
index a4af570ded0..b6d6af1e6da 100644
--- a/source/blender/bmesh/operators/bmo_edgenet.c
+++ b/source/blender/bmesh/operators/bmo_edgenet.c
@@ -362,6 +362,7 @@ static void init_rotsys(BMesh *bm, EdgeData *edata, VertData *vdata)
 	BMEdge **edges = NULL;
 	BLI_array_staticdeclare(edges, BM_DEFAULT_NGON_STACK_SIZE);
 	BMVert *v;
+	RNG *rng;
 	/* BMVert **verts = NULL; */
 	/* BLI_array_staticdeclare(verts, BM_DEFAULT_NGON_STACK_SIZE); */ /* UNUSE */
 	int i;
@@ -427,7 +428,7 @@ static void init_rotsys(BMesh *bm, EdgeData *edata, VertData *vdata)
 			copy_v3_v3(vdata[BM_elem_index_get(v2)].sco, vec1);
 		}
 
-		BLI_srandom(0);
+		rng = BLI_rng_new_srandom(0);
 
 		/* first, ensure no 0 or 180 angles between adjacent
 		 * (and that adjacent's adjacent) edges */
@@ -471,7 +472,7 @@ static void init_rotsys(BMesh *bm, EdgeData *edata, VertData *vdata)
 				copy_v3_v3(cent, v->co);
 
 				for (j = 0; j < 3; j++) {
-					float fac = (BLI_frand() - 0.5f) * size;
+					float fac = (BLI_rng_get_float(rng) - 0.5f) * size;
 					cent[j] += fac;
 				}
 
@@ -488,6 +489,8 @@ static void init_rotsys(BMesh *bm, EdgeData *edata, VertData *vdata)
 			}
 		}
 
+		BLI_rng_free(rng);
+
 		copy_v3_v3(vdata[BM_elem_index_get(v)].offco, cent);
 		//copy_v3_v3(v->co, cent);
 
diff --git a/source/blender/bmesh/operators/bmo_subdivide.c b/source/blender/bmesh/operators/bmo_subdivide.c
index aaa2914e6ae..8df52c3a6ae 100644
--- a/source/blender/bmesh/operators/bmo_subdivide.c
+++ b/source/blender/bmesh/operators/bmo_subdivide.c
@@ -831,11 +831,13 @@ void bmo_subdivide_edges_exec(BMesh *bm, BMOperator *op)
 	params.origkey = skey;
 
 	if (params.use_fractal) {
-		BLI_srandom(seed);
+		RNG *rng = BLI_rng_new_srandom(seed);
 
-		params.fractal_ofs[0] = (float)BLI_drand() * 200.0f;
-		params.fractal_ofs[1] = (float)BLI_drand() * 200.0f;
-		params.fractal_ofs[2] = (float)BLI_drand() * 200.0f;
+		params.fractal_ofs[0] = BLI_rng_get_float(rng) * 200.0f;
+		params.fractal_ofs[1] = BLI_rng_get_float(rng) * 200.0f;
+		params.fractal_ofs[2] = BLI_rng_get_float(rng) * 200.0f;
+
+		BLI_rng_free(rng);
 	}
 	
 	BMO_slot_map_to_flag(bm, op->slots_in, "custom_patterns",
diff --git a/source/blender/editors/curve/editcurve.c b/source/blender/editors/curve/editcurve.c
index 5504b431c50..08cae94fc93 100644
--- a/source/blender/editors/curve/editcurve.c
+++ b/source/blender/editors/curve/editcurve.c
@@ -5402,8 +5402,6 @@ static void selectrandom_curve(ListBase *editnurb, float randfac)
 	BPoint *bp;
 	int a;
 	
-	BLI_srand(BLI_rand()); /* random seed */
-	
 	for (nu = editnurb->first; nu; nu = nu->next) {
 		if (nu->type == CU_BEZIER) {
 			bezt = nu->bezt;
diff --git a/source/blender/editors/gpencil/gpencil_edit.c b/source/blender/editors/gpencil/gpencil_edit.c
index ad32c5fcccf..f12e6da59b2 100644
--- a/source/blender/editors/gpencil/gpencil_edit.c
+++ b/source/blender/editors/gpencil/gpencil_edit.c
@@ -562,7 +562,7 @@ static void gp_timing_data_add_point(tGpTimingData *gtd, double stroke_inittime,
 #define MIN_TIME_DELTA 0.02f
 
 /* Loop over next points to find the end of the stroke, and compute */
-static int gp_find_end_of_stroke_idx(tGpTimingData *gtd, int idx, int nbr_gaps, int *nbr_done_gaps,
+static int gp_find_end_of_stroke_idx(tGpTimingData *gtd, RNG *rng, int idx, int nbr_gaps, int *nbr_done_gaps,
                                      float tot_gaps_time, float delta_time, float *next_delta_time)
 {
 	int j;
@@ -598,7 +598,7 @@ static int gp_find_end_of_stroke_idx(tGpTimingData *gtd, int idx, int nbr_gaps,
 						/* Clamp max between [0.0, gap_randomness], with lower delta giving higher max */
 						max = gtd->gap_randomness - delta;
 						CLAMP(max, 0.0f, gtd->gap_randomness);
-						*next_delta_time += gtd->gap_duration + (BLI_frand() * (max - min)) + min;
+						*next_delta_time += gtd->gap_duration + (BLI_rng_get_float(rng) * (max - min)) + min;
 					}
 				}
 				else {
@@ -613,7 +613,7 @@ static int gp_find_end_of_stroke_idx(tGpTimingData *gtd, int idx, int nbr_gaps,
 	return j - 1;
 }
 
-static void gp_stroke_path_animation_preprocess_gaps(tGpTimingData *gtd, int *nbr_gaps, float *tot_gaps_time)
+static void gp_stroke_path_animation_preprocess_gaps(tGpTimingData *gtd, RNG *rng, int *nbr_gaps, float *tot_gaps_time)
 {
 	int i;
 	float delta_time = 0.0f;
@@ -637,12 +637,12 @@ static void gp_stroke_path_animation_preprocess_gaps(tGpTimingData *gtd, int *nb
 		printf("%f, %f, %f, %d\n", gtd->tot_time, delta_time, *tot_gaps_time, *nbr_gaps);
 	}
 	if (gtd->gap_randomness > 0.0f) {
-		BLI_srandom(gtd->seed);
+		BLI_rng_srandom(rng, gtd->seed);
 	}
 }
 
 static void gp_stroke_path_animation_add_keyframes(ReportList *reports, PointerRNA ptr, PropertyRNA *prop, FCurve *fcu,
-                                                   Curve *cu, tGpTimingData *gtd, float time_range,
+                                                   Curve *cu, tGpTimingData *gtd, RNG *rng, float time_range,
                                                    int nbr_gaps, float tot_gaps_time)
 {
 	/* Use actual recorded timing! */
@@ -669,7 +669,7 @@ static void gp_stroke_path_animation_add_keyframes(ReportList *reports, PointerR
 			start_stroke_idx = i;
 			delta_time = next_delta_time;
 			/* find end of that new stroke */
-			end_stroke_idx = gp_find_end_of_stroke_idx(gtd, i, nbr_gaps, &nbr_done_gaps,
+			end_stroke_idx = gp_find_end_of_stroke_idx(gtd, rng, i, nbr_gaps, &nbr_done_gaps,
 			                                           tot_gaps_time, delta_time, &next_delta_time);
 			/* This one should *never* be negative! */
 			end_stroke_time = time_start + ((gtd->times[end_stroke_idx] + delta_time) / gtd->tot_time * time_range);
@@ -777,6 +777,7 @@ static void gp_stroke_path_animation(bContext *C, ReportList *reports, Curve *cu
 	}
 	else {
 		/* Use actual recorded timing! */
+		RNG *rng = BLI_rng_new(0);
 		float time_range;
 		
 		/* CustomGaps specific */
@@ -784,7 +785,7 @@ static void gp_stroke_path_animation(bContext *C, ReportList *reports, Curve *cu
 		
 		/* Pre-process gaps, in case we don't want to keep their original timing */
 		if (gtd->mode == GP_STROKECONVERT_TIMING_CUSTOMGAP) {
-			gp_stroke_path_animation_preprocess_gaps(gtd, &nbr_gaps, &tot_gaps_time);
+			gp_stroke_path_animation_preprocess_gaps(gtd, rng, &nbr_gaps, &tot_gaps_time);
 		}
 		
 		if (gtd->realtime) {
@@ -798,8 +799,11 @@ static void gp_stroke_path_animation(bContext *C, ReportList *reports, Curve *cu
 			printf("GP Stroke Path Conversion: Starting keying!\n");
 		}
 		
-		gp_stroke_path_animation_add_keyframes(reports, ptr, prop, fcu, cu, gtd, time_range,
+		gp_stroke_path_animation_add_keyframes(reports, ptr, prop, fcu, cu, gtd,
+		                                       rng, time_range,
 		                                       nbr_gaps, tot_gaps_time);
+
+		BLI_rng_free(rng);
 	}
 	
 	/* As we used INSERTKEY_FAST mode, we need to recompute all curve's handles now */
diff --git a/source/blender/editors/mesh/editmesh_knife.c b/source/blender/editors/mesh/editmesh_knife.c
index 69deafabc42..171cbb6859c 100644
--- a/source/blender/editors/mesh/editmesh_knife.c
+++ b/source/blender/editors/mesh/editmesh_knife.c
@@ -1860,12 +1860,12 @@ typedef struct facenet_entry {
 	KnifeEdge *kfe;
 } facenet_entry;
 
-static void rnd_offset_co(float co[3], float scale)
+static void rnd_offset_co(RNG *rng, float co[3], float scale)
 {
 	int i;
 
 	for (i = 0; i < 3; i++) {
-		co[i] += (BLI_frand() - 0.5) * scale;
+		co[i] += (BLI_rng_get_float(rng) - 0.5) * scale;
 	}
 }
 
@@ -1966,6 +1966,7 @@ static void knifenet_fill_faces(KnifeTool_OpData *kcd)
 	BMFace **faces = MEM_callocN(sizeof(BMFace *) * bm->totface, "faces knife");
 	MemArena *arena = BLI_memarena_new(1 << 16, "knifenet_fill_faces");
 	SmallHash shash;
+	RNG *rng;
 	int i, j, k = 0, totface = bm->totface;
 
 	BMO_push(bm, NULL);
@@ -2065,7 +2066,7 @@ static void knifenet_fill_faces(KnifeTool_OpData *kcd)
 		}
 	}
 
-	BLI_srand(0);
+	rng = BLI_rng_new(0);
 
 	for (i = 0; i < totface; i++) {
 		SmallHash *hash = &shash;
@@ -2086,7 +2087,7 @@ static void knifenet_fill_faces(KnifeTool_OpData *kcd)
 			if (!BLI_smallhash_haskey(hash, (intptr_t)entry->kfe->v1)) {
 				sf_vert = BLI_scanfill_vert_add(&sf_ctx, entry->kfe->v1->v->co);
 				sf_vert->poly_nr = 0;
-				rnd_offset_co(sf_vert->co, rndscale);
+				rnd_offset_co(rng, sf_vert->co, rndscale);
 				sf_vert->tmp.p = entry->kfe->v1->v;
 				BLI_smallhash_insert(hash, (intptr_t)entry->kfe->v1, sf_vert);
 			}
@@ -2094,7 +2095,7 @@ static void knifenet_fill_faces(KnifeTool_OpData *kcd)
 			if (!BLI_smallhash_haskey(hash, (intptr_t)entry->kfe->v2)) {
 				sf_vert = BLI_scanfill_vert_add(&sf_ctx, entry->kfe->v2->v->co);
 				sf_vert->poly_nr = 0;
-				rnd_offset_co(sf_vert->co, rndscale);
+				rnd_offset_co(rng, sf_vert->co, rndscale);
 				sf_vert->tmp.p = entry->kfe->v2->v;
 				BLI_smallhash_insert(hash, (intptr_t)entry->kfe->v2, sf_vert);
 			}
@@ -2201,6 +2202,7 @@ static void knifenet_fill_faces(KnifeTool_OpData *kcd)
 	if (faces)
 		MEM_freeN(faces);
 	BLI_memarena_free(arena);
+	BLI_rng_free(rng);
 
 	BMO_error_clear(bm); /* remerge_faces sometimes raises errors, so make sure to clear them */
 
diff --git a/source/blender/editors/mesh/editmesh_select.c b/source/blender/editors/mesh/editmesh_select.c
index 486401eeec9..8d3086c66a3 100644
--- a/source/blender/editors/mesh/editmesh_select.c
+++ b/source/blender/editors/mesh/editmesh_select.c
@@ -3227,8 +3227,6 @@ static int edbm_select_random_exec(bContext *C, wmOperator *op)
 	BMIter iter;
 	const float randfac =  RNA_float_get(op->ptr, "percent") / 100.0f;
 
-	BLI_srand(BLI_rand()); /* random seed */
-	
 	if (!RNA_boolean_get(op->ptr, "extend"))
 		EDBM_flag_disable_all(em, BM_ELEM_SELECT);
 
diff --git a/source/blender/editors/mesh/editmesh_tools.c b/source/blender/editors/mesh/editmesh_tools.c
index 33c665576c6..d43120f998e 100644
--- a/source/blender/editors/mesh/editmesh_tools.c
+++ b/source/blender/editors/mesh/editmesh_tools.c
@@ -3249,7 +3249,7 @@ static void sort_bmelem_flag(Scene *scene, Object *ob,
 		if (totelem[0]) {
 			/* Re-init random generator for each element type, to get consistent random when
 			 * enabling/disabling an element type. */
-			BLI_srandom(seed);
+			RNG *rng = BLI_rng_new_srandom(seed);
 			pb = pblock[0] = MEM_callocN(sizeof(char) * totelem[0], "sort_bmelem vert pblock");
 			sb = sblock[0] = MEM_callocN(sizeof(BMElemSort) * totelem[0], "sort_bmelem vert sblock");
 
@@ -3257,16 +3257,18 @@ static void sort_bmelem_flag(Scene *scene, Object *ob,
 				if (BM_elem_flag_test(ve, flag)) {
 					pb[i] = false;
 					sb[affected[0]].org_idx = i;
-					sb[affected[0]++].srt = BLI_frand();
+					sb[affected[0]++].srt = BLI_rng_get_float(rng);
 				}
 				else {
 					pb[i] = true;
 				}
 			}
+
+			BLI_rng_free(rng);
 		}
 
 		if (totelem[1]) {
-			BLI_srandom(seed);
+			RNG *rng = BLI_rng_new_srandom(seed);
 			pb = pblock[1] = MEM_callocN(sizeof(char) * totelem[1], "sort_bmelem edge pblock");
 			sb = sblock[1] = MEM_callocN(sizeof(BMElemSort) * totelem[1], "sort_bmelem edge sblock");
 
@@ -3274,16 +3276,18 @@ static void sort_bmelem_flag(Scene *scene, Object *ob,
 				if (BM_elem_flag_test(ed, flag)) {
 					pb[i] = false;
 					sb[affected[1]].org_idx = i;
-					sb[affected[1]++].srt = BLI_frand();
+					sb[affected[1]++].srt = BLI_rng_get_float(rng);
 				}
 				else {
 					pb[i] = true;
 				}
 			}
+
+			BLI_rng_free(rng);
 		}
 
 		if (totelem[2]) {
-			BLI_srandom(seed);
+			RNG *rng = BLI_rng_new_srandom(seed);
 			pb = pblock[2] = MEM_callocN(sizeof(char) * totelem[2], "sort_bmelem face pblock");
 			sb = sblock[2] = MEM_callocN(sizeof(BMElemSort) * totelem[2], "sort_bmelem face sblock");
 
@@ -3291,12 +3295,14 @@ static void sort_bmelem_flag(Scene *scene, Object *ob,
 				if (BM_elem_flag_test(fa, flag)) {
 					pb[i] = false;
 					sb[affected[2]].org_idx = i;
-					sb[affected[2]++].srt = BLI_frand();
+					sb[affected[2]++].srt = BLI_rng_get_float(rng);
 				}
 				else {
 					pb[i] = true;
 				}
 			}
+
+			BLI_rng_free(rng);
 		}
 	}
 
diff --git a/source/blender/editors/metaball/mball_edit.c b/source/blender/editors/metaball/mball_edit.c
index 3fe8b93ada3..de0af1d1bb8 100644
--- a/source/blender/editors/metaball/mball_edit.c
+++ b/source/blender/editors/metaball/mball_edit.c
@@ -190,7 +190,6 @@ static int select_random_metaelems_exec(bContext *C, wmOperator *op)
 		return OPERATOR_CANCELLED;
 	
 	ml = mb->editelems->first;
-	BLI_srand(BLI_rand());  /* Random seed */
 	
 	/* Stupid version of random selection. Should be improved. */
 	while (ml) {
diff --git a/source/blender/editors/physics/particle_edit.c b/source/blender/editors/physics/particle_edit.c
index 4825416f6e9..f12df7aa3dd 100644
--- a/source/blender/editors/physics/particle_edit.c
+++ b/source/blender/editors/physics/particle_edit.c
@@ -3254,12 +3254,13 @@ static int brush_add(PEData *data, short number)
 	short size= pset->brush[PE_BRUSH_ADD].size;
 	short size2= size*size;
 	DerivedMesh *dm=0;
+	RNG *rng;
 	invert_m4_m4(imat, ob->obmat);
 
 	if (psys->flag & PSYS_GLOBAL_HAIR)
 		return 0;
 
-	BLI_srandom(psys->seed+data->mval[0]+data->mval[1]);
+	rng = BLI_rng_new_srandom(psys->seed+data->mval[0]+data->mval[1]);
 
 	sim.scene= scene;
 	sim.ob= ob;
@@ -3281,8 +3282,8 @@ static int brush_add(PEData *data, short number)
 
 			/* rejection sampling to get points in circle */
 			while (dmx*dmx + dmy*dmy > size2) {
-				dmx= (2.0f*BLI_frand() - 1.0f)*size;
-				dmy= (2.0f*BLI_frand() - 1.0f)*size;
+				dmx= (2.0f*BLI_rng_get_float(rng) - 1.0f)*size;
+				dmy= (2.0f*BLI_rng_get_float(rng) - 1.0f)*size;
 			}
 		}
 		else {
@@ -3458,6 +3459,8 @@ static int brush_add(PEData *data, short number)
 	if (!psmd->dm->deformedOnly)
 		dm->release(dm);
 	
+	BLI_rng_free(rng);
+	
 	return n;
 }
 
diff --git a/source/blender/editors/transform/transform_generics.c b/source/blender/editors/transform/transform_generics.c
index 43930e3c063..444410cf94e 100644
--- a/source/blender/editors/transform/transform_generics.c
+++ b/source/blender/editors/transform/transform_generics.c
@@ -1759,7 +1759,6 @@ void calculatePropRatio(TransInfo *t)
 						td->factor = (float)sqrt(2 * dist - dist * dist);
 						break;
 					case PROP_RANDOM:
-						BLI_srand(BLI_rand()); /* random seed */
 						td->factor = BLI_frand() * dist;
 						break;
 					default:
diff --git a/source/blender/editors/uvedit/uvedit_unwrap_ops.c b/source/blender/editors/uvedit/uvedit_unwrap_ops.c
index 087e8eed0c9..69f6daeeb19 100644
--- a/source/blender/editors/uvedit/uvedit_unwrap_ops.c
+++ b/source/blender/editors/uvedit/uvedit_unwrap_ops.c
@@ -231,6 +231,7 @@ static ParamHandle *construct_param_handle(Scene *scene, Object *ob, BMEditMesh
 	BMFace *efa;
 	BMLoop *l;
 	BMEdge *eed;
+	RNG *rng;
 	BMIter iter, liter;
 	
 	const int cd_loop_uv_offset = CustomData_get_offset(&bm->ldata, CD_MLOOPUV);
@@ -250,7 +251,7 @@ static ParamHandle *construct_param_handle(Scene *scene, Object *ob, BMEditMesh
 	/* we need the vert indices */
 	BM_mesh_elem_index_ensure(em->bm, BM_VERT);
 
-	BLI_srand(0);
+	rng = BLI_rng_new(0);
 	
 	BM_ITER_MESH (efa, &iter, em->bm, BM_FACES_OF_MESH) {
 		ScanFillVert *sf_vert = NULL, *sf_vert_last, *sf_vert_first;
@@ -311,7 +312,7 @@ static ParamHandle *construct_param_handle(Scene *scene, Object *ob, BMEditMesh
 				
 				/* add small random offset */
 				for (i = 0; i < 3; i++) {
-					sf_vert->co[i] += (BLI_frand() - 0.5f) * FLT_EPSILON * 50;
+					sf_vert->co[i] += (BLI_rng_get_float(rng) - 0.5f) * FLT_EPSILON * 50;
 				}
 				
 				sf_vert->tmp.p = l;
@@ -362,6 +363,7 @@ static ParamHandle *construct_param_handle(Scene *scene, Object *ob, BMEditMesh
 	}
 
 	param_construct_end(handle, fill, implicit);
+	BLI_rng_free(rng);
 
 	return handle;
 }
diff --git a/source/blender/modifiers/intern/MOD_explode.c b/source/blender/modifiers/intern/MOD_explode.c
index 638f8f0ae01..f0eb113e46f 100644
--- a/source/blender/modifiers/intern/MOD_explode.c
+++ b/source/blender/modifiers/intern/MOD_explode.c
@@ -104,6 +104,7 @@ static void createFacepa(ExplodeModifierData *emd,
 	MVert *mvert = NULL;
 	ParticleData *pa;
 	KDTree *tree;
+	RNG *rng;
 	float center[3], co[3];
 	int *facepa = NULL, *vertpa = NULL, totvert = 0, totface = 0, totpart = 0;
 	int i, p, v1, v2, v3, v4 = 0;
@@ -114,7 +115,7 @@ static void createFacepa(ExplodeModifierData *emd,
 	totvert = dm->getNumVerts(dm);
 	totpart = psmd->psys->totpart;
 
-	BLI_srandom(psys->seed);
+	rng = BLI_rng_new_srandom(psys->seed);
 
 	if (emd->facepa)
 		MEM_freeN(emd->facepa);
@@ -136,7 +137,7 @@ static void createFacepa(ExplodeModifierData *emd,
 		if (dvert) {
 			const int defgrp_index = emd->vgroup - 1;
 			for (i = 0; i < totvert; i++, dvert++) {
-				float val = BLI_frand();
+				float val = BLI_rng_get_float(rng);
 				val = (1.0f - emd->protect) * val + emd->protect * 0.5f;
 				if (val < defvert_find_weight(dvert, defgrp_index))
 					vertpa[i] = -1;
@@ -182,6 +183,8 @@ static void createFacepa(ExplodeModifierData *emd,
 
 	if (vertpa) MEM_freeN(vertpa);
 	BLI_kdtree_free(tree);
+
+	BLI_rng_free(rng);
 }
 
 static int edgecut_get(EdgeHash *edgehash, unsigned int v1, unsigned int v2)
diff --git a/source/blender/modifiers/intern/MOD_particleinstance.c b/source/blender/modifiers/intern/MOD_particleinstance.c
index 63dcb5d942d..8a6a83b65b2 100644
--- a/source/blender/modifiers/intern/MOD_particleinstance.c
+++ b/source/blender/modifiers/intern/MOD_particleinstance.c
@@ -313,8 +313,7 @@ static DerivedMesh *applyModifier(ModifierData *md, Object *ob,
 			{
 				float ran = 0.0f;
 				if (pimd->random_position != 0.0f) {
-					BLI_srandom(psys->seed + p);
-					ran = pimd->random_position * BLI_frand();
+					ran = pimd->random_position * BLI_hash_frand(psys->seed + p);
 				}
 
 				if (pimd->flag & eParticleInstanceFlag_KeepShape) {
diff --git a/source/blender/modifiers/intern/MOD_weightvg_util.c b/source/blender/modifiers/intern/MOD_weightvg_util.c
index 7181b5bbb44..a5e63a7832f 100644
--- a/source/blender/modifiers/intern/MOD_weightvg_util.c
+++ b/source/blender/modifiers/intern/MOD_weightvg_util.c
@@ -59,7 +59,7 @@
  * vertex index (in case the weight tables do not cover the whole vertices...).
  * cmap might be NULL, in which case curve mapping mode will return unmodified data.
  */
-void weightvg_do_map(int num, float *new_w, short falloff_type, CurveMapping *cmap)
+void weightvg_do_map(int num, float *new_w, short falloff_type, CurveMapping *cmap, RNG *rng)
 {
 	int i;
 
@@ -100,8 +100,7 @@ void weightvg_do_map(int num, float *new_w, short falloff_type, CurveMapping *cm
 				fac = (float)sqrt(2 * fac - fac * fac);
 				break;
 			case MOD_WVG_MAPPING_RANDOM:
-				BLI_srand(BLI_rand()); /* random seed */
-				fac = BLI_frand() * fac;
+				fac = BLI_rng_get_float(rng) * fac;
 				break;
 			case MOD_WVG_MAPPING_STEP:
 				fac = (fac >= 0.5f) ? 1.0f : 0.0f;
diff --git a/source/blender/modifiers/intern/MOD_weightvg_util.h b/source/blender/modifiers/intern/MOD_weightvg_util.h
index 209b784d573..c7509a937b6 100644
--- a/source/blender/modifiers/intern/MOD_weightvg_util.h
+++ b/source/blender/modifiers/intern/MOD_weightvg_util.h
@@ -32,6 +32,8 @@
 #ifndef __MOD_WEIGHTVG_UTIL_H__
 #define __MOD_WEIGHTVG_UTIL_H__
 
+#include "BLI_rand.h"
+
 /* so modifier types match their defines */
 #include "MOD_modifiertypes.h"
 
@@ -65,7 +67,7 @@ struct Scene;
  * vertex index (in case the weight tables do not cover the whole vertices...).
  * cmap might be NULL, in which case curve mapping mode will return unmodified data.
  */
-void weightvg_do_map(int num, float *new_w, short mode, struct CurveMapping *cmap);
+void weightvg_do_map(int num, float *new_w, short mode, struct CurveMapping *cmap, RNG *rng);
 
 /* Applies new_w weights to org_w ones, using either a texture, vgroup or constant value as factor.
  * Return values are in org_w.
diff --git a/source/blender/modifiers/intern/MOD_weightvgedit.c b/source/blender/modifiers/intern/MOD_weightvgedit.c
index 779ac6b5ecf..dd84dce05c2 100644
--- a/source/blender/modifiers/intern/MOD_weightvgedit.c
+++ b/source/blender/modifiers/intern/MOD_weightvgedit.c
@@ -29,6 +29,7 @@
  */
 
 #include "BLI_utildefines.h"
+#include "BLI_ghash.h"
 #include "BLI_math.h"
 #include "BLI_string.h"
 
@@ -239,7 +240,15 @@ static DerivedMesh *applyModifier(ModifierData *md, Object *ob, DerivedMesh *der
 
 	/* Do mapping. */
 	if (wmd->falloff_type != MOD_WVG_MAPPING_NONE) {
-		weightvg_do_map(numVerts, new_w, wmd->falloff_type, wmd->cmap_curve);
+		RNG *rng = NULL;
+
+		if (wmd->falloff_type == MOD_WVG_MAPPING_RANDOM)
+			rng = BLI_rng_new_srandom(BLI_ghashutil_strhash(ob->id.name));
+
+		weightvg_do_map(numVerts, new_w, wmd->falloff_type, wmd->cmap_curve, rng);
+
+		if (rng)
+			BLI_rng_free(rng);
 	}
 
 	/* Do masking. */
diff --git a/source/blender/modifiers/intern/MOD_weightvgproximity.c b/source/blender/modifiers/intern/MOD_weightvgproximity.c
index 71d6d4880ad..9a6b46f81b8 100644
--- a/source/blender/modifiers/intern/MOD_weightvgproximity.c
+++ b/source/blender/modifiers/intern/MOD_weightvgproximity.c
@@ -30,9 +30,10 @@
 
 #define DO_PROFILE 0
 
+#include "BLI_utildefines.h"
+#include "BLI_ghash.h"
 #include "BLI_math.h"
 #include "BLI_string.h"
-#include "BLI_utildefines.h"
 
 #if DO_PROFILE
 	#include "PIL_time.h"
@@ -185,7 +186,7 @@ static float get_ob2ob_distance(const Object *ob, const Object *obr)
 /**
  * Maps distances to weights, with an optional "smoothing" mapping.
  */
-static void do_map(float *weights, const int nidx, const float min_d, const float max_d, short mode)
+static void do_map(Object *ob, float *weights, const int nidx, const float min_d, const float max_d, short mode)
 {
 	const float range_inv = 1.0f / (max_d - min_d); /* invert since multiplication is faster */
 	unsigned int i = nidx;
@@ -210,7 +211,15 @@ static void do_map(float *weights, const int nidx, const float min_d, const floa
 	}
 
 	if (!ELEM(mode, MOD_WVG_MAPPING_NONE, MOD_WVG_MAPPING_CURVE)) {
-		weightvg_do_map(nidx, weights, mode, NULL);
+		RNG *rng = NULL;
+
+		if (mode == MOD_WVG_MAPPING_RANDOM)
+			rng = BLI_rng_new_srandom(BLI_ghashutil_strhash(ob->id.name));
+
+		weightvg_do_map(nidx, weights, mode, NULL, rng);
+
+		if (rng)
+			BLI_rng_free(rng);
 	}
 }
 
@@ -500,7 +509,7 @@ static DerivedMesh *applyModifier(ModifierData *md, Object *ob, DerivedMesh *der
 	}
 
 	/* Map distances to weights. */
-	do_map(new_w, numIdx, wmd->min_dist, wmd->max_dist, wmd->falloff_type);
+	do_map(ob, new_w, numIdx, wmd->min_dist, wmd->max_dist, wmd->falloff_type);
 
 	/* Do masking. */
 	weightvg_do_mask(numIdx, indices, org_w, new_w, ob, dm, wmd->mask_constant,
diff --git a/source/blender/render/intern/source/convertblender.c b/source/blender/render/intern/source/convertblender.c
index 6b033a4eb3c..76179b25871 100644
--- a/source/blender/render/intern/source/convertblender.c
+++ b/source/blender/render/intern/source/convertblender.c
@@ -176,6 +176,7 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 	Scene *scene;
 	Object *camera;
 	Camera *cam;
+	RNG *rng;
 	double dblrand, hlfrand;
 	float vec[4], fx, fy, fz;
 	float fac, starmindist, clipend;
@@ -244,14 +245,16 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 	if (re) /* add render object for stars */
 		obr= RE_addRenderObject(re, NULL, NULL, 0, 0, 0);
 	
+	rng = BLI_rng_new(0);
+	
 	for (x = sx, fx = sx * stargrid; x <= ex; x++, fx += stargrid) {
 		for (y = sy, fy = sy * stargrid; y <= ey; y++, fy += stargrid) {
 			for (z = sz, fz = sz * stargrid; z <= ez; z++, fz += stargrid) {
 
-				BLI_srand((hash[z & 0xff] << 24) + (hash[y & 0xff] << 16) + (hash[x & 0xff] << 8));
-				vec[0] = fx + (hlfrand * BLI_drand()) - dblrand;
-				vec[1] = fy + (hlfrand * BLI_drand()) - dblrand;
-				vec[2] = fz + (hlfrand * BLI_drand()) - dblrand;
+				BLI_rng_seed(rng, (hash[z & 0xff] << 24) + (hash[y & 0xff] << 16) + (hash[x & 0xff] << 8));
+				vec[0] = fx + (hlfrand * BLI_rng_get_double(rng)) - dblrand;
+				vec[1] = fy + (hlfrand * BLI_rng_get_double(rng)) - dblrand;
+				vec[2] = fz + (hlfrand * BLI_rng_get_double(rng)) - dblrand;
 				vec[3] = 1.0;
 				
 				if (vertexfunc) {
@@ -281,7 +284,7 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 					
 					
 					if (alpha != 0.0f) {
-						fac = force * BLI_drand();
+						fac = force * BLI_rng_get_double(rng);
 						
 						har = initstar(re, obr, vec, fac);
 						
@@ -290,9 +293,9 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 							har->add= 255;
 							har->r = har->g = har->b = 1.0;
 							if (maxjit) {
-								har->r += ((maxjit * BLI_drand()) ) - maxjit;
-								har->g += ((maxjit * BLI_drand()) ) - maxjit;
-								har->b += ((maxjit * BLI_drand()) ) - maxjit;
+								har->r += ((maxjit * BLI_rng_get_double(rng)) ) - maxjit;
+								har->g += ((maxjit * BLI_rng_get_double(rng)) ) - maxjit;
+								har->b += ((maxjit * BLI_rng_get_double(rng)) ) - maxjit;
 							}
 							har->hard = 32;
 							har->lay= -1;
@@ -321,6 +324,8 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 
 	if (obr)
 		re->tothalo += obr->tothalo;
+
+	BLI_rng_free(rng);
 }
 
 
diff --git a/source/blender/render/intern/source/rayshade.c b/source/blender/render/intern/source/rayshade.c
index 88c0719fa59..387004527c0 100644
--- a/source/blender/render/intern/source/rayshade.c
+++ b/source/blender/render/intern/source/rayshade.c
@@ -935,15 +935,14 @@ void init_jitter_plane(LampRen *lar)
 	
 	/* if 1 sample, we leave table to be zero's */
 	if (tot>1) {
+		/* set per-lamp fixed seed */
+		RNG *rng = BLI_rng_new_srandom(tot);
 		int iter=12;
 
-		/* set per-lamp fixed seed */
-		BLI_srandom(tot);
-		
 		/* fill table with random locations, area_size large */
 		for (x=0; x<tot; x++, fp+=2) {
-			fp[0]= (BLI_frand()-0.5f)*lar->area_size;
-			fp[1]= (BLI_frand()-0.5f)*lar->area_sizey;
+			fp[0]= (BLI_rng_get_float(rng)-0.5f)*lar->area_size;
+			fp[1]= (BLI_rng_get_float(rng)-0.5f)*lar->area_sizey;
 		}
 		
 		while (iter--) {
@@ -952,6 +951,8 @@ void init_jitter_plane(LampRen *lar)
 				DP_energy(lar->jitter, fp, tot, lar->area_size, lar->area_sizey);
 			}
 		}
+
+		BLI_rng_free(rng);
 	}
 	/* create the dithered tables (could just check lamp type!) */
 	jitter_plane_offset(lar->jitter, lar->jitter+2*tot, tot, lar->area_size, lar->area_sizey, 0.5f, 0.0f);
@@ -1727,12 +1728,12 @@ static int UNUSED_FUNCTION(ray_trace_shadow_rad)(ShadeInput *ship, ShadeResult *
 }
 
 /* aolight: function to create random unit sphere vectors for total random sampling */
-static void RandomSpherical(float v[3])
+static void RandomSpherical(RNG *rng, float v[3])
 {
 	float r;
-	v[2] = 2.f*BLI_frand()-1.f;
+	v[2] = 2.f*BLI_rng_get_float(rng)-1.f;
 	if ((r = 1.f - v[2]*v[2])>0.f) {
-		float a = 6.283185307f*BLI_frand();
+		float a = 6.283185307f*BLI_rng_get_float(rng);
 		r = sqrt(r);
 		v[0] = r * cosf(a);
 		v[1] = r * sinf(a);
@@ -1770,6 +1771,8 @@ static void DS_energy(float *sphere, int tot, float vec[3])
 /* and allocates threadsafe memory */
 void init_ao_sphere(World *wrld)
 {
+	/* fixed random */
+	RNG *rng;
 	float *fp;
 	int a, tot, iter= 16;
 
@@ -1777,14 +1780,12 @@ void init_ao_sphere(World *wrld)
 	tot= 2*wrld->aosamp*wrld->aosamp;
 	
 	wrld->aosphere= MEM_mallocN(3*tot*sizeof(float), "AO sphere");
-	
-	/* fixed random */
-	BLI_srandom(tot);
+	rng = BLI_rng_new_srandom(tot);
 	
 	/* init */
 	fp= wrld->aosphere;
 	for (a=0; a<tot; a++, fp+= 3) {
-		RandomSpherical(fp);
+		RandomSpherical(rng, fp);
 	}
 	
 	while (iter--) {
@@ -1795,6 +1796,8 @@ void init_ao_sphere(World *wrld)
 	
 	/* tables */
 	wrld->aotables= MEM_mallocN(BLENDER_MAX_THREADS*3*tot*sizeof(float), "AO tables");
+
+	BLI_rng_free(rng);
 }
 
 /* give per thread a table, we have to compare xs ys because of way OSA works... */
@@ -1823,17 +1826,20 @@ static float *sphere_sampler(int type, int resol, int thread, int xs, int ys, in
 	tot= 2*resol*resol;
 
 	if (type & WO_AORNDSMP) {
+		/* total random sampling. NOT THREADSAFE! (should be removed, is not useful) */
+		RNG *rng = BLI_rng_new(BLI_thread_rand(thread));
 		float *sphere;
 		int a;
 		
 		/* always returns table */
 		sphere= threadsafe_table_sphere(0, thread, xs, ys, tot);
 
-		/* total random sampling. NOT THREADSAFE! (should be removed, is not useful) */
 		vec= sphere;
 		for (a=0; a<tot; a++, vec+=3) {
-			RandomSpherical(vec);
+			RandomSpherical(rng, vec);
 		}
+
+		BLI_rng_free(rng);
 		
 		return sphere;
 	}