Progress

1 files changed, 60 insertions, 1 deletions

diff --git a/numalgo.cpp b/numalgo.cpp
index 1db0763..f2ebac5 100644
--- a/numalgo.cpp
+++ b/numalgo.cpp
@@ -1,3 +1,4 @@
+#include <cstdlib>
 #include <cassert>
 #include "numalgo.h"
 
@@ -31,7 +32,7 @@ Bigint egcd(const Bigint &a,const Bigint &b,Bigint &x,Bigint &y){
 Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){
 	assert(e>=0);
 	assert(m>=1);
-	if(m==1)return Bigint(0);
+	if(m==1)return Bigint::zero;
 	Bigint res(1);
 	vector<bool> bits(e.bits());
 	for(int i=bits.size()-1;i>=0;i--){
@@ -42,9 +43,67 @@ Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){
 	return res;
 }
 
+Bigint isqrt(const Bigint &n){
+	assert(n>=0);
+	if(n<=1)return n;
+	const int maxiter=20; //empirically, this should happen around 3.5 million bits in n. (niter ~= -1.87+1.45ln(bits))
+	// __asm("int3\n\t");
+	// cout<<"bitcount="<<n.bitcount()<<" maxiter="<<maxiter<<endl;
+	Bigint x(n);
+	x>>=n.bitcount()/2;
+	// cerr<<"isqrt("<<hex<<n<<"); x = "<<x<<dec<<endl;
+	int iter;
+	for(iter=0;iter<maxiter;iter++){
+		// cerr<<"x="<<hex<<x<<dec<<endl;
+		Bigint xnext(x*x);  // xnext = x - (x*x-n)/(2*x)  [Newton's method]
+		xnext-=n;
+		xnext>>=1;
+		xnext=xnext.divmod(x).first;
+		xnext.negate();
+		xnext+=x;
+		if(xnext==x)break;
+		x=xnext;
+	}
+	cerr<<iter<<endl;
+	assert(iter<maxiter);
+	switch((x*x).compare(n)){
+		case -1:{
+			Bigint x1(x);
+			x1+=Bigint::one;
+			assert(x1*x1>n);
+			return x;
+		}
+		case 0: return x;
+		case 1:{
+			x-=Bigint::one;
+			assert(x*x<=n);
+			return x;
+		}
+		default: assert(false);
+	}
+}
+
 int ilog2(uint64_t i){
 	assert(i);
 	int l=0;
 	while(i>>=1)l++;
 	return l;
 }
+
+Bigint cryptrandom_big(const Bigint &bound){
+	const int blocksize=32;
+	int btc=bound.bitcount();
+	int nblocks=btc/blocksize,rest=btc%blocksize;
+	while(true){
+		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
+		}
+		if(rest){
+			res<<=rest;
+			res+=arc4random_uniform((uint32_t)1<<rest);
+		}
+		if(res<=bound)return res;
+	}
+}