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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
#include <cstdlib>
#include <cassert>
#include "numalgo.h"
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;
if(b==0)return a;
if(a.compareAbs(b)==1)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;
while(r2!=0){
pair<Bigint,Bigint> dm=r.divmod(r2);
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;
if(e==0)return b;
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;
}
int jacobiSymbol(Bigint a,Bigint n){
//https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol
assert(n>0&&n.odd());
int res=1;
while(true){
a=a.divmod(n).second;
if(a<0)a+=n;
int nmod8=n.divmod(Bigint(8)).second.lowdigits();
int two_over_n=nmod8==1||nmod8==7?1:-1;
if(a==0)return 0;
while(a.even()){
res*=two_over_n;
a>>=1;
}
if(a==1)return res;
if(gcd(a,n)!=1)return 0;
if(a.divmod(Bigint(4)).second.lowdigits()==3&&nmod8%4==3)res*=-1;
swap(a,n);
}
}
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))
Bigint x(n);
x>>=n.bitcount()/2;
int iter;
for(iter=0;iter<maxiter;iter++){
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 bigrandom(Rng &rng,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+=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+=rng.get_uniform((uint32_t)1<<rest);
}
if(res<=bound)return res;
}
}
|
