1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
#include <cstdlib>
#include <cassert>
#include "numalgo.h"
using namespace std;
Bigint gcd(Bigint a,Bigint b){
while(true){
if(a==0)return b;
if(b==0)return a;
if(a>=b)a=a.divmod(b).second;
else b=b.divmod(a).second;
}
}
Bigint egcd(const Bigint &a,const Bigint &b,Bigint &x,Bigint &y){
Bigint x2(0),y2(1),r(a),r2(b);
x=1; y=0;
//cerr<<x<<"\t * "<<a<<"\t + "<<y<<"\t * "<<b<<"\t = "<<r<<endl;
while(r2!=0){
pair<Bigint,Bigint> dm=r.divmod(r2);
//cerr<<x2<<"\t * "<<a<<"\t + "<<y2<<"\t * "<<b<<"\t = "<<r2<<" (q = "<<dm.first<<')'<<endl;
Bigint xn=x-dm.first*x2;
Bigint yn=y-dm.first*y2;
x=x2; x2=xn;
y=y2; y2=yn;
r=r2; r2=dm.second;
}
return r;
}
Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){
assert(e>=0);
assert(m>=1);
if(m==1)return Bigint::zero;
Bigint res(1);
vector<bool> bits(e.bits());
for(int i=bits.size()-1;i>=0;i--){
res*=res;
if(bits[i])res*=b;
res=res.divmod(m).second;
}
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;
}
}
|
