数论基础2
文章修改不会置顶看着真难受…
⑤miller-rabin
记n-1=a*2^b。在[1,n)中随机选取一个整数x，如果x^a≡1或x^(a*2^i)≡-1(其中0<=i<b)，那么我们认为n是质数。
正确率大概是((1/4)^随机次数)
下面的代码还带了一个快速乘
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; ll c_f(ll a,ll b,ll c) { if(b==0) return 0; ll hf=c_f(a,b>>1,c); hf=(hf+hf)%c; if(b&1) hf=(hf+a)%c; return hf; } ll cf(ll a,ll b,ll c) { a%=c; b%=c; if(!a||!b) return 0; if(a<0&&b<0) a=-a,b=-b; if(b<0) return c-c_f(a,-b,c); if(a<0) return c-c_f(-a,b,c); return c_f(a,b,c); } ll qp(ll x,ll y,ll m) { if(y==0) return 1; ll hf=qp(x,y>>1,m); hf=cf(hf,hf,m); if(y&1) hf=cf(hf,x%m,m); return hf; } bool isprime(ll p) { int b=0; ll a=p-1; while(!(a&1)) ++b, a>>=1; for(int i=1;i<=10;i++) { ll x=rand()%(p-1)+1; if(qp(x,a,p)==1) goto ct; for(int j=0;j<b;j++) { if(qp(x,a<<j,p)==p-1) goto ct; } return 0; ct:; } return 1; }
⑥pollard-rho
一个比较看脸的算法…
我们先定义一个函数f(x)=(x^2+1)%n。那么如果我一直调用f(x),f(f(x)),f(f(f(x)))…那么就可以生成一个随机数数列对吧。
pollard_rho是一个可以接受一个数n，返回它的一个因数的算法。算法流程是这样的，开始有一个x[1]，在[0,n)里面随机的，我们令x[i]=f(x[i-1])。
我们每次计算gcd(x[k]-x[2k],n)，如果返回n就说明分解失败，需要重新生成一个x[1]。如果返回1就继续把k加1。否则的话我们就找到了n的一个因数。
实现的时候需要注意因为这个随机函数f的一些性质，当n为2的倍数时最好特判掉，否则似乎可能会死循环。
我们来算算复杂度，这个复杂度看起来比较玄学，它可以在O(sqrt(p))的时间复杂度内找到n的一个小因子p，所以复杂度差不多是O($n^{\frac{1}{4}}log$)的（因为还有一个gcd和快速乘(算f的时候由于数据范围的问题一般都会用)）。
那如果要拿来分解质因数，复杂度差不多是O($n^{\frac{1}{4}}log^2$)的…事实上还要写miller-rabin啊，快速乘啊，内存分配啊，常数非常之大…
例4 因数和
给一个比较大的数n，求n的所有因数的和。有多组数据。
t<=50，n<=10^18。
这就是一道pollard-rho裸题。你觉得pollard-rho要跑多久？0.1s？
就是这样，所以pollard-rho的常数平常还是要注意一下…
这份代码用了vector，但是质因子数量不会超过60个，存下来应该没什么太大问题
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; ll c_f(ll a,ll b,ll c) { if(b==0) return 0; ll hf=c_f(a,b>>1,c); hf=(hf+hf)%c; if(b&1) hf=(hf+a)%c; return hf; } ll cf(ll a,ll b,ll c) { a%=c; b%=c; if(!a||!b) return 0; if(a<0&&b<0) a=-a,b=-b; if(b<0) return c-c_f(a,-b,c); if(a<0) return c-c_f(-a,b,c); return c_f(a,b,c); } ll qp(ll x,ll y,ll m) { if(y==0) return 1; ll hf=qp(x,y>>1,m); hf=cf(hf,hf,m); if(y&1) hf=cf(hf,x%m,m); return hf; } bool isprime(ll p) { int b=0; ll a=p-1; while(!(a&1)) ++b, a>>=1; for(int i=1;i<=10;i++) { ll x=rand()%(p-1)+1; if(qp(x,a,p)==1) goto ct; for(int j=0;j<b;j++) { if(qp(x,a*(1<<j),p)==p-1) goto ct; } return 0; ct:; } return 1; } ll f(ll x,ll p) {return (cf(x%p,x%p,p)+1)%p;} ll gcd(ll a,ll b) { if(a<0) a=-a; if(b<0) b=-b; while(b) {ll g=a%b; a=b; b=g;} return a; } ll pollard_rho(ll m) { if(m==1) return 1; if(m%2==0) return 2; if(m%3==0) return 3; if(isprime(m)) return m; s: ll x=rand()%m,y=f(x,m),p=1; while(p==1) { x=f(x,m); y=f(f(y,m),m); p=gcd(x-y,m); } if(p==m) goto s; return p; } #define vll vector<ll> vll merge(vll a,vll b) { vll ans; int as=0,bs=0; while(as<a.size()) ans.push_back(a[as++]); while(bs<b.size()) ans.push_back(b[bs++]); return ans; } vll ysh(ll x) { if(x==1) {return vll();} if(isprime(x)) {vll ans; ans.push_back(x); return ans;} ll ys; vll ans; while((ys=pollard_rho(x))!=1) { vll cur=ysh(ys); while(x%ys==0) x/=ys, ans=merge(ans,cur); } return ans; } ll calc(vll x) { sort(x.begin(),x.end()); int n=x.size(); ll ss=1; for(int i=0;i<n;i++) { ll cur=x[i]; int tot; for(int j=i;j<n;j++) { if(x[j]!=cur) break; tot=j; } int cnt=tot-i+1; ll sum=1,sp=1; for(int i=1;i<=cnt;i++) sp=sp*cur, sum+=sp; ss*=sum; i=tot; } return ss; } int main() { int T; ll x; scanf("%d",&T); while(T--) { scanf("%lld",&x); printf("%lld\n",calc(ysh(x))); } }
⑦中国剩余定理
正常版本：
给n个质数/互质的数，给你一个x除以这些数得到的余数，求最小的x。
比如求一个x使得x mod 3=2，x mod 5=3，x mod 7=5。
我们先求是5、7的倍数且mod 3=2的x，即35p mod 3=2，这个exgcd解一解。然后求是3、7倍数且mod 5=3的，然后是5、7倍数且mod 3=2的，然后加在一起就可以得到满足条件的了。
然后ysy给了一个神奇的公式：
似乎正确性显然。
例5 解方程组
给出n个方程，每个方程形如x mod pi=ai，求最小的满足条件的非负整数x。无解输出-1。
n<=6，pi<=1000（其实就是告诉你答案不会爆long long
我们发现互质的条件不见了，而且还有可能出现相矛盾的情况…
例如x mod 4=2，x mod 2=1就是相矛盾的。
嗯我们似乎只要把每个pi拆分成素数的幂然后就可以转化为普通版本了。
好写吗？ (╯‵□′)╯︵┻━┻ 难写得不行好吗
靠谱一点的做法呢？
我们可以把这些方程合并在一起！
假设我们要合并x mod b1=n1，x mod b2=n2。
设x=b1*k1+n1=b2*k2+n2，那么b1*k1-b2*k2=n2-n1。
我们可以用exgcd解出k1（注意判无解的情况），然后带回得到x=X'，然后我们就可以把两个方程合并为x mod lcm(b1,b2)=X'。
开始令x mod 1=0即可。
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; int n; ll gcd(ll a,ll b) { if(a<0) a=-a; if(b<0) b=-b; while(b) {ll t=a%b; a=b; b=t;} return a; } void ex_gcd(ll a,ll b,ll& x,ll& y) { if(!b) {x=1; y=0; return;} ex_gcd(b,a%b,x,y); ll y_=x-a/b*y; x=y; y=y_; } ll exgcd(ll a,ll b,ll c) { ll gg=gcd(a,b); if(c%gg) return -2333; a/=gg; b/=gg; c/=gg; ll x,y; ex_gcd(a,b,x,y); x=x*c; x%=b; y=(c-x*a)/b; return x; } ll lcm(ll a,ll b) {return a/gcd(a,b)*b;} int main() { ll x=1,y=0,g,h; scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lld%lld",&g,&h); long long p=exgcd(x,g,h-y); if(p==-2333) {puts("-1"); return 0;} long long X=p*x+y; long long lc=lcm(x,g); y=(X%lc+lc)%lc; x=lc; } printf("%lld\n",y); }
