RNG

7 files changed, 262 insertions, 12 deletions

diff --git a/dieharder.log b/dieharder.log
new file mode 100644
index 0000000..521165f
--- /dev/null
+++ b/dieharder.log
@@ -0,0 +1,126 @@
+#=============================================================================#
+#            dieharder version 3.31.1 Copyright 2003 Robert G. Brown          #
+#=============================================================================#
+   rng_name    |rands/second|   Seed   |
+stdin_input_raw|  3.38e+06  |2347499450|
+#=============================================================================#
+        test_name   |ntup| tsamples |psamples|  p-value |Assessment
+#=============================================================================#
+   diehard_birthdays|   0|       100|     100|0.49421610|  PASSED  
+      diehard_operm5|   0|   1000000|     100|0.52934238|  PASSED  
+  diehard_rank_32x32|   0|     40000|     100|0.41152731|  PASSED  
+    diehard_rank_6x8|   0|    100000|     100|0.26737294|  PASSED  
+   diehard_bitstream|   0|   2097152|     100|0.49100139|  PASSED  
+        diehard_opso|   0|   2097152|     100|0.66025098|  PASSED  
+        diehard_oqso|   0|   2097152|     100|0.40375900|  PASSED  
+         diehard_dna|   0|   2097152|     100|0.88608427|  PASSED  
+diehard_count_1s_str|   0|    256000|     100|0.27316362|  PASSED  
+diehard_count_1s_byt|   0|    256000|     100|0.39290164|  PASSED  
+ diehard_parking_lot|   0|     12000|     100|0.20394610|  PASSED  
+    diehard_2dsphere|   2|      8000|     100|0.01965863|  PASSED  
+    diehard_3dsphere|   3|      4000|     100|0.85892974|  PASSED  
+     diehard_squeeze|   0|    100000|     100|0.22082844|  PASSED  
+        diehard_sums|   0|       100|     100|0.04890097|  PASSED  
+        diehard_runs|   0|    100000|     100|0.67362967|  PASSED  
+        diehard_runs|   0|    100000|     100|0.82034207|  PASSED  
+       diehard_craps|   0|    200000|     100|0.02670837|  PASSED  
+       diehard_craps|   0|    200000|     100|0.96939866|  PASSED  
+ marsaglia_tsang_gcd|   0|  10000000|     100|0.16443903|  PASSED  
+ marsaglia_tsang_gcd|   0|  10000000|     100|0.25080776|  PASSED  
+         sts_monobit|   1|    100000|     100|0.00702961|  PASSED  
+            sts_runs|   2|    100000|     100|0.06480441|  PASSED  
+          sts_serial|   1|    100000|     100|0.43938171|  PASSED  
+          sts_serial|   2|    100000|     100|0.74652004|  PASSED  
+          sts_serial|   3|    100000|     100|0.16547027|  PASSED  
+          sts_serial|   3|    100000|     100|0.98663544|  PASSED  
+          sts_serial|   4|    100000|     100|0.61753898|  PASSED  
+          sts_serial|   4|    100000|     100|0.57109674|  PASSED  
+          sts_serial|   5|    100000|     100|0.37997456|  PASSED  
+          sts_serial|   5|    100000|     100|0.60279130|  PASSED  
+          sts_serial|   6|    100000|     100|0.82950080|  PASSED  
+          sts_serial|   6|    100000|     100|0.86583141|  PASSED  
+          sts_serial|   7|    100000|     100|0.81386613|  PASSED  
+          sts_serial|   7|    100000|     100|0.83747318|  PASSED  
+          sts_serial|   8|    100000|     100|0.97078446|  PASSED  
+          sts_serial|   8|    100000|     100|0.70028526|  PASSED  
+          sts_serial|   9|    100000|     100|0.92745231|  PASSED  
+          sts_serial|   9|    100000|     100|0.59325539|  PASSED  
+          sts_serial|  10|    100000|     100|0.04285136|  PASSED  
+          sts_serial|  10|    100000|     100|0.11661089|  PASSED  
+          sts_serial|  11|    100000|     100|0.12385561|  PASSED  
+          sts_serial|  11|    100000|     100|0.91784478|  PASSED  
+          sts_serial|  12|    100000|     100|0.68463426|  PASSED  
+          sts_serial|  12|    100000|     100|0.69381397|  PASSED  
+          sts_serial|  13|    100000|     100|0.99884773|   WEAK   
+          sts_serial|  13|    100000|     100|0.93079127|  PASSED  
+          sts_serial|  14|    100000|     100|0.97999449|  PASSED  
+          sts_serial|  14|    100000|     100|0.99473370|  PASSED  
+          sts_serial|  15|    100000|     100|0.84337184|  PASSED  
+          sts_serial|  15|    100000|     100|0.29415279|  PASSED  
+          sts_serial|  16|    100000|     100|0.81612024|  PASSED  
+          sts_serial|  16|    100000|     100|0.33395808|  PASSED  
+         rgb_bitdist|   1|    100000|     100|0.41943699|  PASSED  
+         rgb_bitdist|   2|    100000|     100|0.49250080|  PASSED  
+         rgb_bitdist|   3|    100000|     100|0.78579363|  PASSED  
+         rgb_bitdist|   4|    100000|     100|0.30229642|  PASSED  
+         rgb_bitdist|   5|    100000|     100|0.73430267|  PASSED  
+         rgb_bitdist|   6|    100000|     100|0.44467027|  PASSED  
+         rgb_bitdist|   7|    100000|     100|0.97087091|  PASSED  
+         rgb_bitdist|   8|    100000|     100|0.87789484|  PASSED  
+         rgb_bitdist|   9|    100000|     100|0.00074196|   WEAK   
+         rgb_bitdist|  10|    100000|     100|0.98411066|  PASSED  
+         rgb_bitdist|  11|    100000|     100|0.31503608|  PASSED  
+         rgb_bitdist|  12|    100000|     100|0.54356466|  PASSED  
+rgb_minimum_distance|   2|     10000|    1000|0.35012021|  PASSED  
+rgb_minimum_distance|   3|     10000|    1000|0.24429438|  PASSED  
+rgb_minimum_distance|   4|     10000|    1000|0.01902755|  PASSED  
+rgb_minimum_distance|   5|     10000|    1000|0.04974865|  PASSED  
+    rgb_permutations|   2|    100000|     100|0.75467096|  PASSED  
+    rgb_permutations|   3|    100000|     100|0.16612371|  PASSED  
+    rgb_permutations|   4|    100000|     100|0.38642879|  PASSED  
+    rgb_permutations|   5|    100000|     100|0.53728843|  PASSED  
+      rgb_lagged_sum|   0|   1000000|     100|0.78799315|  PASSED  
+      rgb_lagged_sum|   1|   1000000|     100|0.37399744|  PASSED  
+      rgb_lagged_sum|   2|   1000000|     100|0.26974480|  PASSED  
+      rgb_lagged_sum|   3|   1000000|     100|0.44903045|  PASSED  
+      rgb_lagged_sum|   4|   1000000|     100|0.44080282|  PASSED  
+      rgb_lagged_sum|   5|   1000000|     100|0.53532929|  PASSED  
+      rgb_lagged_sum|   6|   1000000|     100|0.19790769|  PASSED  
+      rgb_lagged_sum|   7|   1000000|     100|0.62366283|  PASSED  
+      rgb_lagged_sum|   8|   1000000|     100|0.14882099|  PASSED  
+      rgb_lagged_sum|   9|   1000000|     100|0.34536926|  PASSED  
+      rgb_lagged_sum|  10|   1000000|     100|0.97409197|  PASSED  
+      rgb_lagged_sum|  11|   1000000|     100|0.48991082|  PASSED  
+      rgb_lagged_sum|  12|   1000000|     100|0.92693373|  PASSED  
+      rgb_lagged_sum|  13|   1000000|     100|0.57406082|  PASSED  
+      rgb_lagged_sum|  14|   1000000|     100|0.14004367|  PASSED  
+      rgb_lagged_sum|  15|   1000000|     100|0.84545538|  PASSED  
+      rgb_lagged_sum|  16|   1000000|     100|0.96953295|  PASSED  
+      rgb_lagged_sum|  17|   1000000|     100|0.27891981|  PASSED  
+      rgb_lagged_sum|  18|   1000000|     100|0.65386545|  PASSED  
+      rgb_lagged_sum|  19|   1000000|     100|0.39720386|  PASSED  
+      rgb_lagged_sum|  20|   1000000|     100|0.29712584|  PASSED  
+      rgb_lagged_sum|  21|   1000000|     100|0.54751627|  PASSED  
+      rgb_lagged_sum|  22|   1000000|     100|0.15215071|  PASSED  
+      rgb_lagged_sum|  23|   1000000|     100|0.76166879|  PASSED  
+      rgb_lagged_sum|  24|   1000000|     100|0.76841514|  PASSED  
+      rgb_lagged_sum|  25|   1000000|     100|0.24292236|  PASSED  
+      rgb_lagged_sum|  26|   1000000|     100|0.94622491|  PASSED  
+      rgb_lagged_sum|  27|   1000000|     100|0.52948133|  PASSED  
+      rgb_lagged_sum|  28|   1000000|     100|0.51460911|  PASSED  
+      rgb_lagged_sum|  29|   1000000|     100|0.37631497|  PASSED  
+      rgb_lagged_sum|  30|   1000000|     100|0.74393636|  PASSED  
+      rgb_lagged_sum|  31|   1000000|     100|0.52389710|  PASSED  
+      rgb_lagged_sum|  32|   1000000|     100|0.91530195|  PASSED  
+     rgb_kstest_test|   0|     10000|    1000|0.73641402|  PASSED  
+     dab_bytedistrib|   0|  51200000|       1|0.28787304|  PASSED  
+             dab_dct| 256|     50000|       1|0.15114674|  PASSED  
+Preparing to run test 207.  ntuple = 0
+        dab_filltree|  32|  15000000|       1|0.12914597|  PASSED  
+        dab_filltree|  32|  15000000|       1|0.84989818|  PASSED  
+Preparing to run test 208.  ntuple = 0
+       dab_filltree2|   0|   5000000|       1|0.48498152|  PASSED  
+       dab_filltree2|   1|   5000000|       1|0.68790380|  PASSED  
+Preparing to run test 209.  ntuple = 0
+        dab_monobit2|  12|  65000000|       1|0.94661207|  PASSED  
+
diff --git a/numalgo.cpp b/numalgo.cpp
index c3bfa1b..38184b1 100644
--- a/numalgo.cpp
+++ b/numalgo.cpp
@@ -4,6 +4,15 @@
 
 using namespace std;
 
+int64_t gcd(int64_t a,int64_t b){
+	while(true){
+		if(a==0)return b;
+		if(b==0)return a;
+		if(abs(a)>abs(b))a%=b;
+		else b%=a;
+	}
+}
+
 Bigint gcd(Bigint a,Bigint b){
 	while(true){
 		if(a==0)return b;
@@ -113,7 +122,7 @@ int ilog2(uint64_t i){
 	return l;
 }
 
-Bigint cryptrandom_big(const Bigint &bound){
+Bigint bigrandom(Rng &rng,const Bigint &bound){
 	const int blocksize=32;
 	int btc=bound.bitcount();
 	int nblocks=btc/blocksize,rest=btc%blocksize;
@@ -121,11 +130,11 @@ Bigint cryptrandom_big(const Bigint &bound){
 		Bigint res;
 		for(int i=0;i<nblocks;i++){
 			if(i!=0)res<<=blocksize;
-			res+=arc4random_uniform((uint32_t)(((uint64_t)1<<blocksize)-1)); //make sure we don't shift out of our int
+			res+=rng.get_uniform((uint32_t)(((uint64_t)1<<blocksize)-1)); //make sure we don't shift out of our int
 		}
 		if(rest){
 			res<<=rest;
-			res+=arc4random_uniform((uint32_t)1<<rest);
+			res+=rng.get_uniform((uint32_t)1<<rest);
 		}
 		if(res<=bound)return res;
 	}
diff --git a/numalgo.h b/numalgo.h
index 2304e94..9f108e6 100644
--- a/numalgo.h
+++ b/numalgo.h
@@ -2,6 +2,9 @@
 
 #include <cstdint>
 #include "bigint.h"
+#include "rng.h"
+
+int64_t gcd(int64_t a,int64_t b);
 
 Bigint gcd(Bigint a,Bigint b);
 Bigint egcd(const Bigint &a,const Bigint &b,Bigint &x,Bigint &y);
@@ -15,4 +18,4 @@ Bigint isqrt(const Bigint &n);
 
 int ilog2(uint64_t i);
 
-Bigint cryptrandom_big(const Bigint &upperbound); //Return value in [0,upperbound]
+Bigint bigrandom(Rng &rng,const Bigint &upperbound); //Return value in [0,upperbound]
diff --git a/primes.cpp b/primes.cpp
index 934d886..ddb559c 100644
--- a/primes.cpp
+++ b/primes.cpp
@@ -31,7 +31,7 @@ void fillsmallprimes(){
 	//cerr<<endl;
 }
 
-pair<Bigint,Bigint> genprimepair(int nbits){
+pair<Bigint,Bigint> genprimepair(Rng &rng,int nbits){
 	// for x = nbits/2:
 	// (2^x)^2 = 2^(2x)
 	// (2^x + 2^(x-2))^2 = 2^(2x) + 2^(2x-1) + 2^(2x-4)
@@ -41,11 +41,11 @@ pair<Bigint,Bigint> genprimepair(int nbits){
 	int x1=nbits/2-2,x2=(nbits+1)/2+2;
 	assert(x1+x2==nbits);
 	return make_pair(
-		randprime(Bigint::one<<x1,(Bigint::one<<x1)+(Bigint::one<<(x1-2))),
-		randprime(Bigint::one<<x2,(Bigint::one<<x2)+(Bigint::one<<(x2-2))));
+		randprime(rng,Bigint::one<<x1,(Bigint::one<<x1)+(Bigint::one<<(x1-2))),
+		randprime(rng,Bigint::one<<x2,(Bigint::one<<x2)+(Bigint::one<<(x2-2))));
 }
 
-Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
+Bigint randprime(Rng &rng,const Bigint &biglow,const Bigint &bighigh){
 	//https://en.wikipedia.org/wiki/Generating_primes#Large_primes
 	if(!smallprimes_inited)fillsmallprimes();
 
@@ -58,7 +58,7 @@ Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
 		high=bighigh;
 		// cerr<<"low=biglow="<<low<<" high=bighigh="<<high<<endl;
 	} else {
-		high=low=cryptrandom_big(diff-maxrangesize);
+		high=low=bigrandom(rng,diff-maxrangesize);
 		high+=maxrangesize;
 		// cerr<<"low="<<low<<" high="<<high<<endl;
 	}
@@ -102,7 +102,7 @@ Bigint randprime(const Bigint &biglow,const Bigint &bighigh){
 	// cerr<<endl;
 
 	while(maybeprimes.size()){
-		int idx=arc4random_uniform(maybeprimes.size());
+		int idx=rng.get_uniform(maybeprimes.size());
 		int i=maybeprimes[idx];
 		Bigint bi(low+2*i);
 		if(bailliePSW(bi))return bi;
diff --git a/primes.h b/primes.h
index eb817b1..0be8aa7 100644
--- a/primes.h
+++ b/primes.h
@@ -9,11 +9,11 @@ extern std::vector<int> smallprimes;
 void fillsmallprimes();
 
 //for use in RSA (pass target number of bits of N)
-std::pair<Bigint,Bigint> genprimepair(int nbits);
+std::pair<Bigint,Bigint> genprimepair(Rng &rng,int nbits);
 
 //finds random in range [low,high]; throws range_error("No primes") if no prime found
 //Will call fillsmallprimes() if not yet done
-Bigint randprime(const Bigint &low,const Bigint &high);
+Bigint randprime(Rng &rng,const Bigint &low,const Bigint &high);
 
 //checks Fermat [pseudo- or actual] primality
 bool strongPseudoPrime2(const Bigint &n);
diff --git a/rng.cpp b/rng.cpp
new file mode 100644
index 0000000..c9e7f03
--- /dev/null
+++ b/rng.cpp
@@ -0,0 +1,78 @@
+// #include <iostream>
+// #include <bitset>
+#include <cstring>
+#include <cassert>
+#include "rng.h"
+
+using namespace std;
+
+
+//adapted from http://blog.regehr.org/archives/1063 (rotl32c version)
+inline uint64_t rotl64(uint64_t x,uint32_t n){
+	assert(n<64);
+	return (x<<n)|(x>>(-n&63));
+}
+inline uint64_t rotr64(uint64_t x,uint32_t n){
+	assert(n<64);
+	return (x>>n)|(x<<(-n&63));
+}
+
+KeyRng::KeyRng(const char *key_,int keylen_)
+		:keylen(keylen_),idx(0),state(0){
+	assert(keylen>0);
+	assert(key_);
+	key=new uint8_t[keylen];
+	memcpy(key,key_,keylen);
+	stir();
+}
+
+KeyRng::KeyRng(const string &key_)
+		:keylen(key_.size()),idx(0),state(0){
+	assert(keylen>0);
+	key=new uint8_t[keylen];
+	memcpy(key,key_.data(),keylen);
+	stir();
+}
+
+KeyRng::~KeyRng(){
+	delete[] key;
+}
+
+void KeyRng::stir(){
+	for(int i=0;i<10;i++)get();
+}
+
+uint32_t KeyRng::get(){
+	//Progressively and cyclicly mixes in the key with the state.
+	//Mostly own creation; the "tempering" part is derived from a (part of a) Mersenne Twister implementation.
+	//In terms of distribution of values, this is similar to arc4random().
+	//This passes DieHarder 3.31.1 (tested 2016/10/06) with 2 "weak" results. The initial key "wachtwoord" was used.
+	//At least n subsequent values should be independent, where n is the length of the initial key.
+	state^=key[idx];
+	state+=17;
+	key[idx]+=13;
+	idx=(idx+1)%keylen;
+	state^=rotr64(state,11); //tempering
+	state^=rotl64(state,7)&0x9d2c5680;
+	state^=rotr64(state,18);
+	// cerr<<bitset<64>(state)<<endl;
+	return state>>32;
+}
+
+uint32_t KeyRng::get_uniform(uint32_t upbound){
+	if(upbound<=1)return 0;
+	uint32_t min=((uint64_t)1<<32)%upbound; //this is the amount of unusable RNG outputs when avoiding modulo bias
+	while(true){
+		uint32_t v=get();
+		if(v>=min)return v%upbound;
+	}
+}
+
+
+uint32_t Arc4Rng::get(){
+	return arc4random();
+}
+
+uint32_t Arc4Rng::get_uniform(uint32_t upbound){
+	return arc4random_uniform(upbound);
+}
diff --git a/rng.h b/rng.h
new file mode 100644
index 0000000..105e390
--- /dev/null
+++ b/rng.h
@@ -0,0 +1,34 @@
+#pragma once
+
+#include <string>
+#include <cstdint>
+
+class Rng{
+public:
+	virtual uint32_t get()=0;
+	virtual uint32_t get_uniform(uint32_t upbound)=0;
+};
+
+class KeyRng : public Rng{
+	uint8_t *key;
+	int keylen;
+	int idx;
+	uint64_t state;
+
+	void stir();
+
+public:
+	KeyRng(const char *key,int keylen);
+	explicit KeyRng(const std::string &key);
+	KeyRng(const Rng&)=delete; //just keep it at one KeyRng please
+	~KeyRng();
+
+	uint32_t get();
+	uint32_t get_uniform(uint32_t upbound);
+};
+
+class Arc4Rng : public Rng{
+public:
+	uint32_t get();
+	uint32_t get_uniform(uint32_t upbound);
+};