Bzoj3481 DZY Loves Math III

Posted 2020-09-22 SilverNebula

Time Limit: 5 Sec  Memory Limit: 64 MB
Submit: 310  Solved: 65

Description

Input

Output

Sample Input

3
1
2
3
2
4
2

Sample Output

6

HINT

 

1<=N<=10,0<=Qi<=10^18,1<=Pi<=10^18,P>=2

本题仅四组数据。

 

Source

By Jc

 

数学问题 欧拉函数 Miller-Rabin Pollard-rho

花了一整晚来怼这道题……在long long的领域遇到了许多问题。

 

假设我们有充足的时间枚举每一个x，那么在x确定的情况下，原式变成了一个模方程。
根据裴蜀定理，我们知道只有当$ gcd(x,P)|Q $ 的情况下方程有 $ gcd(x,P) $ 个解。
现在我们可以愉快地枚举每一个gcd，那么答案就是：
$$ans=\sum_{d|P,d|Q} d * \sum_{x}[gcd(x,\frac{P}{d})==1]$$
后面那个sum显然是欧拉函数
那么$$ ans=\sum_{d|P,d|Q} d · \varphi(\frac{P}{d}) $$
分(an)析(zhao)一(tao)波(lu)，发现这是一个积性函数，所以我们可以分别考虑每个质因子的贡献，再把它们累乘起来。
我们现在考虑一个质因子p，它在P中有$q$个，在Q中有$q‘$个：
它对答案的贡献是
$$ \sum_{i=0}^{q‘} p^i * \varphi(p^{q-i})$$
我们知道
$$\varphi(p^{q})=p^{q}·\frac{p-1}{p}$$
所以上面的sum可以化成：
$$p^{q-1}·(p-1)·(q‘+1)$$
但是有一个特殊情况，当i=q的时候,$\varphi(1)=1$，不能表示成$\frac{p-1}{p}$的形式，而我们却把它用这种形式算进去了。
也就是说我们把一个$p^q$ 的贡献算成了 $(p-1)p^{q-1}$，特判一下消除即可。

 

  1 #include<iostream>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<cstdio>
  5 #include<cmath>
  6 #include<map>
  7 #define LL long long
  8 using namespace std;
  9 const int mxn=100010;
 10 LL read(){
 11     LL x=0,f=1;char ch=getchar();
 12     while(ch<‘0‘ || ch>‘9‘){if(ch==‘-‘)f=-1;ch=getchar();}
 13     while(ch>=‘0‘ && ch<=‘9‘){x=x*10-‘0‘+ch;ch=getchar();}
 14     return x*f;
 15 }
 16 LL Rand(LL n){return (((LL)rand()<<31)|rand())%(n-1)+1;}
 17 //
 18 map<LL,int>primap;
 19 LL pri[mxn*10];int cnt=0;
 20 bool vis[mxn];
 21 void init(){
 22     int mxn=200;
 23     for(int i=2;i<mxn;i++){
 24         if(!vis[i]){
 25             pri[++cnt]=i;
 26             primap[i]=cnt;
 27         }
 28         for(int j=1;j<=cnt && pri[j]*i<mxn;j++){
 29             vis[pri[j]*i]=1;
 30             if(i%pri[j]==0)break;
 31         }
 32     }
 33     return;
 34 }
 35 //
 36 LL mul(LL x,LL k,LL P)
 37 {
 38     LL res=0;
 39     while(k){
 40         if(k&1)res=(res+x)%P;
 41         x=(x+x)%P;
 42         k>>=1;
 43     }
 44     return res;
 45 }
 46 /*
 47 LL mul(LL a,LL b,LL P){
 48     LL d=(long double)a/P*b+1e-8;
 49     
 50     a=a*b-d*P;
 51     a=(a<0)?a+P:a;
 52     printf("%lld \n",a);
 53     return a;
 54 }*/
 55 LL ksm(LL a,LL k,LL P){
 56     LL res=1;
 57     while(k){
 58         if(k&1)res=mul(res,a,P);
 59         a=mul(a,a,P);
 60         k>>=1;
 61     }
 62     return res;
 63 }
 64 ///
 65 bool check(LL a,LL n,LL r,LL s){
 66     LL res=ksm(a,r,n);LL b=res;
 67     for(int i=1;i<=s;i++){
 68         res=mul(res,res,n);
 69         if(res==1 && b!=1 && b!=n-1)return 1;
 70         b=res;
 71     }
 72     return (res!=1);
 73 }
 74 bool MR(LL n){//素数测试 
 75     if(n<=1)return 0;
 76     if(n==2)return 1;
 77     if(~n&1)return 0;
 78     LL r=n-1,s=0;
 79     while(!(r&1))r>>=1,s++;
 80     for(int i=1;i<=10;i++)
 81         if(check(Rand(n),n,r,s))return 0;
 82     return 1;
 83 }
 84 ///
 85 LL p[mxn],q[mxn];
 86 void addpri_P(LL x){
 87     if(primap.count(x)){
 88         ++p[primap[x]];return;
 89     }
 90     pri[++cnt]=x;
 91     primap[x]=cnt;
 92     p[cnt]=1;
 93     return;
 94 }
 95 void addpri_Q(LL x){
 96     if(primap.count(x)){
 97         int t=primap[x];
 98         if(q[t]<p[t])++q[t];
 99     }
100     return;
101 }
102 LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
103 LL Rho(LL n,LL c){
104     if(n<0)while(1);
105     LL k=2,x=Rand(n),y=x,p=1;
106     for(LL i=1;p==1;i++){
107         x=(mul(x,x,n)+c)%n;
108         p=(y>x)?(y-x):(x-y);
109         p=gcd(n,p);
110         if(i==k)y=x,k<<=1;
111     }
112     return p;
113 }
114 void solve(LL x,bool flag){//分解质因数 
115     if(x==1)return;
116     if(MR(x)){
117         if(!flag)addpri_P(x);//P
118         else addpri_Q(x);//Q
119         return;
120     }
121     LL t=x;
122     while(t==x)t=Rho(x,Rand(x));
123     solve(t,flag);
124     solve(x/t,flag);
125     return;
126 }
127 //
128 const int mod=1e9+7;
129 int n;
130 LL P[50],Q[50];
131 void Calc(){
132     LL ans=1;
133     for(int i=1;i<=cnt;i++){
134         if(!p[i])continue;
135         LL R=mul((q[i]+1),(pri[i]-1),mod);
136         if(p[i]==q[i])R++;
137         R=mul(R,ksm(pri[i],p[i]-1,mod),mod);
138         ans=mul(ans,R,mod);
139     }
140     printf("%lld\n",ans);
141     return;
142 }
143 int main(){
144     srand(19260817);
145     init();
146     int i,j;
147     n=read();
148     for(i=1;i<=n;i++){
149         P[i]=read();
150         solve(P[i],0);
151     }
152     for(i=1;i<=n;i++){
153         Q[i]=read();
154         if(!Q[i]){//特判0 
155             for(j=1;j<=cnt;j++)q[j]=p[j];
156             break;
157         }
158         else solve(Q[i],1);
159     }
160     Calc();
161     return 0;
162 }