BZOJ2627 JZPKIL

            一句话题意:求$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$,多组数据,模$1e9+7$。

            30%的数据,$x=y$
            另30%的数据,$n<=10^9$, $x, y<=100$
       100%的数据,$T<=100, 1<=n<=10^{18}, 0<=x, y<=3000$

 

       orz出题人顾宇宙!!!orz全网题解的源头ydc!!!感觉ydc的题解已经写得很清楚了,但是为了内部交流还是写一写……下面的公式中$[x]$表示$x==1$。

 

推导部分

Part1 原式整理$$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$$$$=\sum_{i=1}^n(in)^ygcd(i,n)^{x-y}$$$$=n^y\sum_{i=1}^ni^ygcd(i,n)^{x-y}$$$$=n^y\sum_{d|n}d^{x-y}\sum_{k=1,[gcd(k,\frac{n}{d}))]}^{\frac{n}{d}}(kd)^y$$$$=n^y\sum_{d|n}d^x\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$$

 

Part2 套路反演

          这一部分我们来整理上式最后的$\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$

$$\sum_{i=1,[gcd(i,n)]}^ni^y$$$$=\sum_{i=1}^ni^y\sum_{d|gcd(i,n)}\mu(d)$$$$=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\frac{n}{d}}(id)^y$$$$=\sum_{d=1}^n\mu(d)d^y\sum_{i=1}^{\frac{n}{d}}i^y$$

 

Part3 伯努利数

       可以发现Part2的最后我们把Part1的幂之和化为了从1开始的连续一段……这有什么用呢?我们可以套结论了!有一个叫伯努利数的东西(下文推导中用$B_k$表示第$k$个伯努利数),在网上搜一搜的话资料不少,百度百科里基本的公式也给了一些。刚开始看证明说什么从生成函数定义的角度可以证明它和自然数幂和的数学关系,但是我后来看到它就是在求自然数幂和的过程中被发明的就感觉证这个没必要了……求伯努利数的方法放在实现部分里,现在我们继续关注推导过程。

$$\sum_{i=1}^ni^m=\frac{1}{m+1}\sum_{k=0}^mC_{m+1,k}B_kn^{m-k+1}$$

大多数情况下伯努利数都是这么用的……可以看成公式

你可能会发现伯努利数的公式很乱,有人用$n$的幂有人用$n+1$的幂还有人带-1的幂……然而因为不会证我也不能确切地说哪一种是对的。维基百科上大概是说伯努利数有两类,两类伯努利数的区别在于B1。第一类伯努利数,由美国国家标准技术研究所 (NIST)制定,在这标准下B1=-1/2;第二类伯努利数,又被称为是“原始的伯努利数”,在这标准下B1=1/2。$n+1$的式子是关于第一类伯努利数的,而$n$的式子是关于第二类伯努利数的;第一类伯努利数乘上-1的幂即为第二类伯努利数……假装他们都是对的好了编套理论真难啊

设$t_k=\frac{1}{m+1}C_{m+1,k}B_k$,则$\sum_{i=1}^ni^m=\sum_{k=0}^mt_kn^{m-k+1}$

 

Part4 积性函数

       即使化成这步代入原式它也不好处理……不过看一眼$10^{18}$的数据范围,分解质因数来做好像是一种历史必然?so现在我们试一试把积性函数的部分整理到一起会发生什么……

$$n^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y\sum_{k=0}^yt_k(\frac{n}{dd_1})^{y-k+1}$$

$$=\sum_{k=0}^yt_kn^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y(\frac{n}{dd_1})^{y-k+1}$$

       这时候除去$t_k$的部分就是一个积性函数了,$y$很小只有$3000$,我们可以枚举$k$然后分解质因数求后面的部分。那么就变成了对每一个质因子$p_i$的$k_i$幂$p_k$求$p_k^y\sum_{d|p_k}d^x\sum_{d_1|\frac{p_k}{d}}\mu(d_1)d_1^y(\frac{p_k}{dd_1})^{y-k+1}$,又因为$\mu(d_1)$只有在$d_1=1$时为1、$d_1=p_i$时为-1而其他情况下都为0,式子可以直接变成$p_k^y\sum_{d|p_k}d^x((\frac{p_k}{d})^{y-k+1}-p_i^y\frac{p_k}{dp_i}^{y-k+1})$,这就是推式子的终点了。

 

 

实现部分

Part1 伯努利数

         你问我伯努利数怎么求啊……   

$\sum_{k=0}^nC_{n+1,k}B_k=n+1$  

所以杨辉三角求一下组合数就可以$O(k^2)$求出所有伯努利数了。

 

Part2 分解质因数

        给$10^{18}$分解质因数?需要$n^{\frac{1}{4}}$的算法!Miller-Rabin和Pollard-Rho可以解决,关于这两个算法我们的朋友Doggu的blog堪称全网最良心。

 

Part3 快速乘

         $10^{18}*10^{18}$应该怎么做?快速乘显然是必要的。这是ydc改良过的快速乘板子。

1 inline LL mul(LL x,LL y,LL z)
2 {
3     return (x*y-(LL)(((long double)x*y+0.5)/(long double)z)*z+z)%z;
4 }

 

 

好像少了点什么东西

        没有错,这题是有部分分的吧?对于$x=y$的$30$分,式子就是$\sum_{i=1}^n(in)^y$,如果你会用伯努利数就可以拿到这一部分的分。对于$n<=10^9$的$30$分,伯努利数仍然是必要的,不过分解质因数可以$\sqrt{n}$来做。看来出题人只打算给会伯努利数的人分数(显然我不在此之列,所以部分分和正解一样困难……)。想打部分分的同学可以去cogs上提交。

 

 

最后放代码,完结撒花。
  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cstring>
  4 #include<cstdlib>
  5 #include<algorithm>
  6 #define ll long long
  7 using namespace std;
  8 const int mod=1e9+7;
  9 const int sj=3010;
 10 const int ss=110;
 11 inline ll read()
 12 {
 13     ll jg=0;int jk=getchar()-'0';
 14     while(jk<0||jk>9)    jk=getchar()-'0';
 15     while(jk>=0&&jk<=9)  jg*=10,jg+=jk,jk=getchar()-'0';
 16     return jg;
 17 }
 18 int ca,ai,bi,cnt,tot,s[ss],sp,t[sj][sj],c[sj][sj];
 19 ll n,a1,a2,b[sj],pi[ss],pk[ss];
 20 ll p[sj],prime[10]={2,3,5,7,11,13,17,19,23,29},fac[ss];
 21 ll pw[ss][ss],pw1[ss];
 22 bool vi[sj];
 23 ll gcd(ll x,ll y)
 24 {
 25     if(!y)  return x;
 26     return gcd(y,x%y);
 27 }
 28 inline ll mul(ll x,ll y,ll z)
 29 {
 30     return (x*y-(ll)(((long double)x*y+0.5)/(long double)z)*z+z)%z;
 31 }
 32 ll qpow(ll x,ll y,ll z)
 33 {
 34     ll ret=1;x%=z;
 35     while(y)
 36     {
 37         if(y&1)  ret=mul(ret,x,z);
 38         y>>=1,x=mul(x,x,z);   
 39     }
 40     return ret;
 41 }
 42 ll ksm(ll x,int y)
 43 {
 44     ll ret=1;x%=mod;
 45     while(y)
 46     {
 47         if(y&1)  ret=ret*x%mod;
 48         x=x*x%mod,y>>=1;
 49     }
 50     return ret;
 51 }
 52 void pre()
 53 {
 54     for(int i=0;i<=3005;i++)
 55     {
 56         c[i][0]=1;
 57         for(int j=1;j<=i;j++)
 58         {    
 59             c[i][j]=c[i-1][j]+c[i-1][j-1];
 60             if(c[i][j]>=mod)  c[i][j]-=mod;
 61         }
 62     }
 63     for(int i=0;i<=3005;i++)
 64     {
 65         b[i]=i+1;
 66         for(int j=0;j<i;j++)
 67         {
 68             b[i]=(b[i]-c[i+1][j]*b[j])%mod;
 69             if(b[i]<0)  b[i]+=mod;
 70         }
 71         b[i]=b[i]*ksm(c[i+1][i],mod-2)%mod;
 72     }
 73     for(int i=0;i<=3005;i++)
 74     {
 75         a1=ksm(i+1,mod-2);
 76         for(int j=0;j<=i;j++)
 77             t[i][j]=c[i+1][j]*b[j]%mod*a1%mod;
 78     }
 79     for(int i=2;i<=1000;i++) 
 80         for(int j=2;j<=1000;j++)
 81             if(i*j<=1000)
 82                 vi[i*j]=1;
 83     for(int i=2;i<=1000;i++)
 84         if(!vi[i])  p[++cnt]=i;
 85 }
 86 inline ll calc(ll x)
 87 {
 88     ll ret=0,q=x%mod;
 89     ll nt=q;
 90     for(int i=bi;i>=0;i--)
 91     {
 92         ret=(ret+t[bi][i]*nt)%mod;
 93         nt=nt*q%mod;
 94     }
 95     return ret*ksm(q,bi)%mod;
 96 }
 97 bool Miller_Rabin(ll x)
 98 {
 99     if(x==1)  return 0;
100     for(int j=0;j<=9;j++)
101         if(x==prime[j])   
102             return 1;
103     ll y=x-1;
104     int k=0;
105     while((y&1)==0)  k++,y>>=1;
106     for(int i=0;i<10;i++)
107     {
108         ll z=rand()%(x-1)+1;
109         ll c=qpow(z,y,x),d;
110         for(int j=0;j<k;j++)
111         {
112             d=mul(c,c,x);
113             if(d==1&&c!=1&&c!=x-1)  return 0;
114             c=d;
115         }
116         if(d!=1)  return 0;
117     }
118     return 1;
119 }
120 ll pollard_rho(ll x,ll y)
121 {
122     ll i=1,k=2,c=rand()%(x-1)+1;
123     ll d=c;
124     while(1)
125     {
126         i++;
127         c=(mul(c,c,x)+y)%x;
128         ll g=gcd((d-c+x)%x,x);
129         if(g!=1&&g!=x)  return g;
130         if(c==d)  return x;
131         if(i==k)  d=c,k<<=1;
132     }
133 }
134 void find(ll x,int tim)
135 {
136     if(x==1)  return;
137     if(Miller_Rabin(x))
138     {
139         fac[++tot]=x;
140         return;
141     }
142     ll yz=x,mk=tim;
143     while(yz>=x)  yz=pollard_rho(yz,tim--);
144     find(yz,mk),find(x/yz,mk);
145 }
146 ll solve(ll x,ll y,ll z)
147 {
148     ll ret=0,qwq,s1,s2;
149     tot=sp=0;
150     for(int j=1;j<=cnt;j++)
151         while(x%p[j]==0)
152             fac[++tot]=p[j],x/=p[j];
153     find(x,120);
154     sort(fac+1,fac+tot+1);
155     for(int i=1;i<=tot;i++)
156     {
157         if(fac[i]!=fac[i-1])  sp++,pi[sp]=pk[sp]=fac[i],s[sp]=1;
158         else  pk[sp]*=fac[i],s[sp]++;
159     }
160     for(int i=1;i<=sp;i++)
161     {
162         pw[i][0]=1,a1=ksm(pi[i],y);
163         for(int j=1;j<=s[i];j++)
164             pw[i][j]=pw[i][j-1]*a1%mod;
165     }
166     for(int i=0;i<=z;i++)
167     {
168         qwq=1;
169         for(int j=1;j<=sp;j++)
170         {
171             s1=s2=0,pw1[0]=1,a1=ksm(pi[j],z),a2=ksm(pi[j],z-i+1);
172             for(int k=1;k<=s[j];k++)  pw1[k]=pw1[k-1]*a2%mod;
173             for(int k=0;k<=s[j];k++)  s1=(s1+pw[j][k]*pw1[s[j]-k])%mod;
174             for(int k=0;k<s[j];k++)   s2=(s2+pw[j][k]*pw1[s[j]-k-1]%mod*a1)%mod;
175             qwq=qwq*(s1-s2+mod)%mod*ksm(pk[j],z)%mod;
176         }
177         ret=(ret+qwq*t[z][i])%mod;
178     }
179     return ret;
180 }
181 int main()
182 {
183     pre();
184     ca=read();
185     for(int i=1;i<=ca;i++)
186     {
187         n=read(),ai=read(),bi=read();
188         if(ai==bi)  printf("%lld\n",calc(n));
189         else        printf("%lld\n",solve(n,ai,bi));
190     }
191     return 0;
192 }
JZPKIL

 

 

 

posted @ 2018-04-25 07:01  moyiii  阅读(240)  评论(0编辑  收藏