题目：https://www.lydsy.com/JudgeOnline/problem.php?id=3481
推推式子发现：令Q=gcd(P,Q)，ans=Σ(d|Q) d*phi(P/d)。把 d 质因数分解，设 t 为 Q 的指数， w 为 P 的指数，ans变成每个质数的 Σ(i=0~t) p^i * phi( p^(w-i) ) 连乘。
分解质因数用 Pollar Rho 。
注意 Q=0 就是 Q=P，要特判！而且不要以为答案变成 (!x || !y) 了！
d从0到P-1 就是 d从1到P！不要特判 P==Q时给答案减P ！因为算的时候就没算d=0的！
每个质数的那个式子可以化简，把 phi 写开，和 p^i 合并，就不用枚举 i 。但要注意 w-i ==0 时 phi 的式子有些不同了。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<cstdlib> #include<ctime> #define ll long long using namespace std; const int N=650,mod=1e9+7; int n,ans=1,c[N],c2[N],tot; ll p[N],P=1; bool vis[N],flag; ll rdl() { ll ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9') ret=(ret<<3ll)+(ret<<1ll)+ch-'0',ch=getchar(); return fx?ret:-ret; } void upd(ll &x,ll md){x-=(x>=md?md:0);} ll mul(ll a,ll b,ll md) { ll ret=0;//0! while(b) { if(b&1ll)ret=ret+a,upd(ret,md); a=a+a,upd(a,md);b>>=1ll; } return ret; } ll pw(ll x,ll k,ll md) { ll ret=1; while(k) { if(k&1ll)ret=mul(ret,x,md); x=mul(x,x,md);k>>=1ll; } return ret; } int phi(ll p,int k) { if(!k) return 1; return mul(p-1,pw(p,k-1,mod),mod); } void add(ll a,bool flag) { if(flag) { for(int i=1;i<=tot;i++) if(p[i]==a) { c[i]++;return;} p[++tot]=a; c[tot]=1; } else { for(int i=1;i<=tot;i++) if(p[i]==a&&c2[i]<c[i]) c2[i]++;//Q=gcd() } } bool ml_rb(ll n) { if(n<2)return false;if(n==2)return true;if((n&1)==0)return false; ll u=n-1,t=0; while((u&1)==0){u>>=1,t++;} int s=20; while(s--) { ll a=rand()%(n-2)+2;//2~n-1 a=pw(a,u,n); ll pre=a; for(int i=1;i<=t;i++) { a=mul(a,a,n); if(a==1&&pre!=1&&pre!=n-1)return false; pre=a; } if(a!=1) return false; } return true; } ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} ll pl_rho(ll x,ll c) { ll x0=rand()%x,y0=x, k=1,i=0; while(1) { x0=(mul(x0,x0,x)+c)%x; ll g=gcd(abs(x0-y0),x); if(g!=1&&g!=x) return g; if(x0==y0)return x; if((++i)==k){k<<=1;y0=x0;} } } void fd_fac(ll n,bool flag) { if(n<2)return; if(ml_rb(n)) { add(n,flag);return; } ll p=n; while(p==n)p=pl_rho(p,rand()%(n-1)+1); fd_fac(p,flag); fd_fac(n/p,flag); } int main() { srand(time(0)); n=rdl(); for(int i=1;i<=n;i++) { ll a=rdl(); P=mul(P,a,mod); fd_fac(a,1); } for(int i=1;i<=n;i++) { ll a=rdl(); if(!a){flag=1;continue;} if(flag)continue; fd_fac(a,0); } if(flag)for(int i=1;i<=tot;i++)c2[i]=c[i]; for(int i=1,d,tp;i<=tot;i++) { tp=pw(p[i],c[i]-1,mod); d=mul(tp,mul(p[i]-1,c2[i]+1,mod)+(c[i]==c2[i]),mod); ans=mul(ans,d,mod); } printf("%d ",ans); return 0; }