BZOJ 2627 JZPKIL

Posted 2020-10-10 avancent

题面：

2627: JZPKIL

Time Limit: 50 Sec  Memory Limit: 128 MB
Submit: 347  Solved: 104
[Submit][Status][Discuss]

Description

Input

第一行，询问个数T。
　　下面T行，每行三个整数，n, x, y。

Output

　　T行，每行一个整数，表示相应的询问的答案

Sample Input

5
6 0 0
6 0 1
6 1 0
6 1 1
1000000000 50 50

Sample Output

6
66
15
126
393442025

HINT

数据规模和约定
30%的数据，x=y
另30%的数据，n<=10^9, x, y<=100
100%的数据，T<=100, 1<=n<=10^18, 0<=x, y<=3000

 

这题是个50sec卡评测反演好题。

$\quad \sum_{i=1}^{n}(i,n)^{x}[i,n]^{y}$

$=\sum_{i=1}^{n}(i,n)^{x-y}(i\cdot n)^{y}$

$=n^{y}\sum_{i=1}^{n}(i,n)^{x-y}i^{y}$

$=n^{y}\sum_{d|n}d^{x}\sum_{i=1}^{\frac{n}{d}}\sum_{k|i,k|\frac{n}{d}}\mu(k)i^{y}$

$=n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}\sum_{i=1}^{\frac{n}{k\cdot d}i^{y}}$

$=n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}\sum_{i=0}^{y}\frac{1}{y+1}\dbinom{y+1}{i}B_{i}(\frac{n}{k\cdot d})^{y+1-i}$

$=\sum_{i=0}^{y}t_{i}n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}(\frac{n}{k\cdot d})^{y+1-i}$

$t_{i}=\frac{1}{y+1}\dbinom{y+1}{i}B_{i}$ ($B_{i}$为伯努利数)

将$n$分解质因数$n=p_{1}^{k_1}\cdot p_{2}^{k_2}\cdot p_{3}^{k_3}\cdots $

令$f(p_{j}^{k_j},i)=\sum_{d|p_{j}^{k_j}}d^{p_{j}^{k_j}}\sum_{k|\frac{p_{j}^{k_j}}{d}}\mu(k)k^{y}(\frac{p_{j}^{k_j}}{k\cdot d})^{y+1-i}$

$f$为积性函数,可以用每个质因子的答案相乘得到$n$的答案。

但是由于$n$很大,我们只能用$Pollard_rho$分解大数因子。

时间复杂度为$O(y^{2}+T\cdot n^{\frac{1}{4}}+T\cdot y\cdot log\;n)$

  1 #include <iostream>
  2 #include <stdio.h>
  3 #include <algorithm>
  4 using namespace std;
  5 #define LL long long
  6 template<typename __>
  7 inline void read(__ &s)
  8 {
  9     char ch=getchar();
 10     for(s=0;ch<‘0‘||ch>‘9‘;ch=getchar());
 11     for(;ch>=‘0‘&&ch<=‘9‘;s=s*10+ch-‘0‘,ch=getchar());
 12 }
 13 namespace Pollard_rho
 14 {
 15     inline LL mul(LL x,LL y,LL mod)
 16     {
 17         return (x*y-(LL)(x/(long double)mod*y+0.5)*mod+mod)%mod;
 18     }
 19     inline LL qpow(LL x,LL y,LL mod)
 20     {
 21         LL ans=1;
 22         for(;y;y>>=1,x=mul(x,x,mod))
 23             if(y&1)
 24                 ans=mul(ans,x,mod);
 25         return ans;
 26     }
 27     inline bool judge(LL n,LL p)
 28     {
 29         if(n<=p)
 30             return 1;
 31         if(n%p==0)
 32             return 0;
 33         LL x=n-1;
 34         int s=-1;
 35         while(!(x&1))
 36             x>>=1,++s;
 37         x=qpow(p,x,n);
 38         if(x==1||x==n-1)
 39             return 1;
 40         while(s--)
 41         {
 42             x=mul(x,x,n);
 43             if(x==1)
 44                 return 0;
 45             if(x==n-1)
 46                 return 1;
 47         }
 48         return 0;
 49     }
 50     inline bool miller_rabin(LL n)
 51     {
 52         return judge(n,2)&&judge(n,3)&&judge(n,5)&&judge(n,7)&&judge(n,11)&&judge(n,13)&&judge(n,17)&&judge(n,19)&&judge(n,23);
 53     }
 54     LL *P;
 55     int Pcnt;
 56     inline LL pollard_rho(LL n)
 57     {
 58         LL x=rand()%(n-1)+1,xx=mul(x,x,n)+1,t=1;
 59         while(t==1)
 60         {
 61             x=mul(x,x,n)+1;
 62             xx=mul(xx,xx,n)+1;
 63             xx=mul(xx,xx,n)+1;
 64             t=__gcd(x-xx+n,n);
 65         }
 66         return t;
 67     }
 68     inline void divide(LL n)
 69     {
 70         if(n==1)
 71             return ;
 72         else
 73             if(n%2==0)
 74                 P[++Pcnt]=2,divide(n/2);
 75             else
 76                 if(miller_rabin(n))
 77                     P[++Pcnt]=n;
 78                 else
 79                 {
 80                     LL x=pollard_rho(n);
 81                     divide(x);
 82                     divide(n/x);
 83                 }
 84     }
 85     inline void init(LL n,LL *p,int &cnt)
 86     {
 87         P=p;
 88         Pcnt=0;
 89         divide(n);
 90         cnt=Pcnt;
 91     }
 92 }
 93 #define mod 1000000007
 94 inline int mul(int x,int y)
 95 {
 96     return ((LL)x*y)%mod;
 97 }
 98 inline int add(int x,int y)
 99 {
100     x+=y;
101     if(x>=mod)
102         x-=mod;
103     return x;
104 }
105 inline int dec(int x,int y)
106 {
107     x-=y;
108     if(x<0)
109         x+=mod;
110     return x;
111 }
112 inline int qpow(int x,int y)
113 {
114     x%=mod;
115     int ans=1;
116     for(;y;y>>=1,x=mul(x,x))
117         if(y&1)
118             ans=mul(ans,x);
119     return ans;
120 }
121 #define maxn 3002
122 int B[maxn],C[maxn][maxn],a[maxn],pw[20][maxn],iv[20][maxn];
123 int p[20],c[20],x,y,N;
124 LL n;
125 int tot;
126 inline void init()
127 {
128     int i,j;
129     C[0][0]=1;
130     for(i=1;i<=3001;i++)
131     {
132         C[i][0]=1;
133         for(j=1;j<=i;j++)
134             C[i][j]=add(C[i-1][j],C[i-1][j-1]);  
135     }
136     B[0]=1;
137     for(i=1;i<=3001;i++)
138     {
139         for(j=0;j<i;++j)
140             B[i]=dec(B[i],mul(C[i+1][j],B[j]));
141         B[i]=mul(B[i],qpow(i+1,mod-2));
142     }
143 }
144 inline void pre()
145 {
146     for(int i=0;i<=y+1;i++)
147         a[i]=0;
148     for(int i=0;i<=y;i++)
149         a[y+1-i]=mul(C[y+1][i],B[i]);
150     int tmp=qpow(y+1,mod-2);
151     for(int i=0;i<=y+1;i++)
152         a[i]=mul(a[i],tmp);
153     a[y]=add(a[y],1);
154     if(!y)
155         a[0]=dec(a[0],1);
156 }
157 inline void get_frac()
158 {
159     static LL frac[97];
160     int cnt=0;
161     tot=0;
162     Pollard_rho::init(n,frac,cnt);
163     sort(frac+1,frac+cnt+1);
164     for(int i=1;i<=cnt;i++)
165         if(frac[i]!=frac[i-1])
166             p[++tot]=frac[i]%mod,c[tot]=1;
167         else
168             ++c[tot];
169     for(int i=1;i<=tot;i++)
170     {
171         pw[i][0]=iv[i][0]=1;
172         int t=max(max(x,y),c[i])+1;
173         for(int j=1;j<=t;j++)
174             pw[i][j]=mul(pw[i][j-1],p[i]);
175         int v=qpow(p[i],mod-2);
176         for(int j=1;j<=t;j++)
177             iv[i][j]=mul(iv[i][j-1],v);
178     }
179 }
180 inline int geo(int i,int x)
181 {
182     int c=::c[i];
183     int p;
184     if(x<0)
185         p=iv[i][-x];
186     else
187         p=pw[i][x];
188     if(p==1)
189         return c;
190     int ans=0,t=p;
191     for(int j=1;j<=c;j++)
192     {
193         ans=add(ans,t);
194         t=mul(t,p);
195     }
196     return ans;
197 }
198 inline void input()
199 {
200     read(n);
201     read(x);
202     read(y);
203     unsigned int rand_seed=18230742;
204     rand_seed+=(n+(x<<16)+y)+(rand_seed<<2);
205     srand(rand_seed);
206 }
207 inline int calc(int x,int y)
208 {
209     int ans=x<0?qpow(qpow(N,mod-2),-x):qpow(N,x);
210     for(int i=1;i<=tot;i++)
211     {
212         int t1=geo(i,-x);
213         int t2=dec(1,y>=0?pw[i][y]:qpow(p[i],mod-2));
214         ans=mul(ans,add(1,mul(t1,t2)));
215     }
216     return ans;
217 }
218 void solve()
219 {
220     pre();
221     get_frac();
222     N=n%mod;
223     if(!N)
224         return (void)(puts("0"));
225     int ans=0;
226     int t=1,tmp;
227     for(int i=0;i<=y+1;i++)
228     {
229         tmp=calc(x-i,y-i);
230         ans=add(ans,mul(mul(a[i],t),tmp));
231         t=mul(t,N);
232     }
233     ans=mul(ans,qpow(N,y));
234     printf("%d\n",ans);
235 }
236 int main()
237 {
238     init();
239     int T;
240     read(T);
241     while(T--)
242     {
243         input();
244         solve();
245     }
246 }

BZOJ 2627

PS：$Pollard_rho$写丑了会TLE,需要大力卡一波常数才能过。