【HDU6051】if the starlight never fade

題意：

給出m，P，求$\sum\limits_{i=1}^{P-1}if(i)$，其中$f(i)=\sum\limits_{x=1}^{P-1}\sum\limits_{y=1}^{m}[(x+y)^i=(x^i)(\mod P)]$

$1\leq m\leq P-1\leq 10^9+6$

題解：

神仙數學題。。。

因爲$P$是奇質數，因此確定存在原根$g$，那麼確定惟一存在$a$，$b$知足$1\leq a,b\leq P-1$使得$x=g^a$，$y=g^b$，因此有：

$(x+y)^i≡x^i(\mod P)$

$(g^{a}+g^{b})^i≡g^{ai}(\mod P)$

$(1+g^{b-a})^i≡1(\mod P)$

因爲$1+g^{b-a}\geq 2$，因此一定存在惟一的$k$知足$1\leq k\leq P-1$使得$g^k≡1+g^{b-a}(\mod P)$；

顯然$g^{ki}≡1(\mod P)$，因此$ki≡0(\mod P-1)$；

那麼$k=a\times\frac{P-1}{gcd(P-1,i)}$，其中$a=1,2,3,...,gcd(P-1,i)-1$；

因此知足條件的$k$有$gcd(P-1,i)-1$個，固定$y$，則有$x≡y(g^{k}-1)^{-1}(\mod P)$；

易知若$y(g^{k_1}-1)^{-1}≡y(g^{k_2}-1)^{-1}(\mod P)$，則$g^{k_1}≡g^{k_2}(\mod P)$，即$k_1=k_2$；

因此對於肯定的$y$，$x$會有$gcd(P-1,i)-1$種不一樣取值。

至此式子化爲$f(i)=m(gcd(P-1,i))$

而後大力推一波式子：

$\sum\limits_{i=1}^{P-1}if(i)=m\sum\limits_{i=1}^{P-1}i((gcd(P-1),i)-1)=m\sum\limits_{i=1}^{P-1}i(gcd(P-1),i)-m\frac{P\times(P-1)}{2}$

把前半部分拿出來：

$m\sum\limits_{i=1}^{P-1}i(gcd(P-1),i)$

$=\sum\limits_{d|(p-1)}d\sum\limits_{d|i,1\leq i\leq P-1}i[gcd(P-1,i)=d]$

$=\sum\limits_{d|(p-1)}d^2\sum\limits_{i=1}^{\frac{P-1}{d}}i[gcd(\frac{P-1}{d},i)=1]$

$=\sum\limits_{d|(p-1)}d^2\frac{\frac{P-1}{d}\varphi(\frac{P-1}{d})+[\frac{P-1}{d}=1]}{2}$（引理）

這樣複雜度就是$O(\sqrt{n})$的，能夠直接作。

引理：$\sum\limits_{i=1}^{n}i[gcd(i,n)=1]=\frac{1}{2}\sum\limits_{i=1}^{n}(i+(n-i))[gcd(i,n)=1]=\frac{1}{2}n\sum\limits_{i=1}^{n}[gcd(i,n)=1]$

$=\frac{n\varphi(n)+[n=1]}{2}$

證實：若$gcd(i,n)=1$，那麼$gcd(n-i,n)=1$。

代碼：

 1 #include<algorithm>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cstdio>
 5 #include<cmath>
 6 #define eps 1e-4
 7 #define mod 1000000007
 8 #define S(n) (ll)(n*(n+1)/2)
 9 using namespace std; 10 typedef long long ll; 11 ll phi(int x){ 12     ll ret=x; 13     for(ll i=2;i*i<=x;i++){ 14         if(!(x%i)){ 15             ret=ret/i*(i-1); 16             while(!(x%i))x/=i; 17  } 18  } 19     if(x>1)ret=ret/x*(x-1); 20     return ret; 21 } 22 ll f(ll x,ll d){ 23     return d*d%mod*((x*phi(x)+(x==1))/2)%mod; 24 } 25 ll calc(ll x){ 26     ll ret=0; 27     for(ll i=1;i*i<=x;i++){ 28         if(!(x%i)){ 29             ret=(ret+f(x/i,i))%mod; 30             if(i*i!=x)ret=(ret+f(i,x/i))%mod; 31  } 32  } 33     ret=(ret-x*(x+1)/2%mod+mod)%mod; 34     return ret; 35 } 36 ll m,p; 37 int main(){ 38     scanf("%lld%lld",&m,&p); 39     printf("%lld",calc(p-1)*m%mod); 40     return 0; 41 }