HDU 6051 If the starlight never fade（原根+推式子）

Posted 2021-02-02 xu-daxia

题目大意：

设(f(i))为使((x+y)^i equiv x^i (mod p))成立的(x,y)的对数。其中(1 leq x leq p-1 , 1leq yleq m),m，p给定且p是一个质数。求(sum_{i=1}^{p-1}i*f(i)),p<=1e9+7,m<=p-1

思路

我们考虑用原根去代换x，y。
设g为p的一个原根，(g^aequiv x(mod p),g^b equiv y(mod p))。
然后我们用(g)去代换(x,y)。((g^a+g^b)^iequiv g^{a*i}(mod p))然后我们在式子两边同时除一个(g^{a*i})。得到((1+g^{b-a})^iequiv 1(mod p))
设(1+g^{b-a}equiv g^k(mod p))(因为原根的性质所以我们一定可以找到这样的k)。
此时原式为(g^{k*i} equiv 1(mod p))由费马小定理可以到到(p-1 mid k*i)。
这就要求k为(frac{p-1}{gcd(p-1,i)})的倍数(即至少包含p-1所特有的因子)。由于(0<k<p-1),为什么k不能取0呢？因为$g^{b-a}不会为0，所以(1+g^{b-1}> 1)。
所以可以得到这样的k的数量是(frac{p-1}{frac{p-1}{gcd(p-1,i)}}-1=gcd(p-1,i)-1)，这里因为k不能取0和p-1，所以要减1。
又因为(1+g^{b-a} equiv g^k(mod p))得(g^b equiv(g^k-1)*g^a (mod p))，(yequiv x*(g^k-1)(mod p))。每有一个k，y对应一个x。所以对于一个y有(gcd(p-1,i)-1)个x对应。
所以可以得出(f(i)=m*(gcd(p-1,i)-1))。
(sum_{i=1}^{p-1}i*f(i))
(=msum_{i=1}^{p-1}i*gcd(p-1,i)-m*frac{(p-1)*p}{2})
重点是如何求(sum_{i=1}^{p-1}i*gcd(p-1,i))。
(sum_{i=1}^{p-1}i*gcd(p-1,i))
显然d是p-1的约数
(=sum_{dmid p-1}dsum_{i=1}^{p-1}i*[gcd(p-1,i)==d])
(=sum_{dmid p-1}d^2sum_{i=1}^{frac{p-1}{d}}i*[gcd(frac{p-1}{d},i)==1])
然后有一个神奇的变换。
(sum_{i=1}^{t}i*[gcd(t,i)==1]=frac{t*varphi(t)+[t==1]}{2})
为什么？
令(t=frac{p-1}{d})就是(sum_{i=1}^{t}i*[gcd(t,i)==1])
其实就是求1到t中与t互质的数的和。
由更相减损术得若(gcd(n,i)=1)则(gcd(n,n-i)=1)
所以一个与(t)互质的数(x)，(t-x)也与(t)互质。
所以与t互质的数成对出现，设这一对数为a,b，有(a+b=t)。
所以(sum_{i=1}^{t}i*[gcd(t,i)==1]=frac{t*varphi(t)+[t==1]}{2})
故原式可以化为
(sum_{dmid p-1}d^2sum_{i=1}^{frac{p-1}{d}}i*[gcd(frac{p-1}{d},i)==1])
(=sum_{dmid p-1}d^2*frac{frac{p-1}{d}*varphi(frac{p-1}{d})+[frac{p-1}{d}==1]}{2})
带回去求解即可。
复杂度？(O(能过))，因为求约数可以先扫出质数来优化，所以(sqrt{frac{p}{lnp}}*sqrt{p}=sqrt{frac{p^2}{lnp}})差不多两亿，不是很大。而且求欧拉函数是求p-1的约数的约数，一定比(sqrt{p})小，所以能过？

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define int long long
const int N=1001000;
const int mod=1e9+7;
bool book[N];
int prime[N],cnt,inv,T;
int read(){
    int sum=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
    return sum*f;
}
void pre_work(){
    for(int i=2;i<=1000000;i++){
        if(book[i]==0)prime[++cnt]=i;
        for(int j=1;j<=cnt&&prime[j]*i<=1000000;j++){
            book[prime[j]*i]=1;
            if(i%prime[j]==0)break;
        }
    }
}
int ksm(int x,int b){
    int tmp=1;
    while(b){
        if(b&1)tmp=tmp*x%mod;
        b>>=1;
        x=x*x%mod;
    }
    return tmp;
}
int phi(int x){
    int tmp=x;
    int ans=x;
    for(int i=1;i<=cnt&&prime[i]*prime[i]<=x;i++){
        if(tmp%prime[i]==0){
            ans=ans/prime[i]*(prime[i]-1);
            while(tmp%prime[i]==0)tmp/=prime[i];
        }
    }
    if(tmp>1)ans=ans/tmp*(tmp-1);
    return ans;
}
int work(int x){
    int tmp=0;
    for(int i=1;i*i<=x;i++){
        if(x%i==0){
            int a=i*i%mod;
            int b=x/i;
            int c=phi(b);
            int d=(b==1);
            tmp=(tmp+a*(b*c%mod+d)%mod*inv%mod)%mod;
            if(x/i>i){
                int hh=x/i;
                int a=hh*hh%mod;
                int b=x/hh;
                int c=phi(b);
                int d=(b==1);
                tmp=(tmp+a*(b*c%mod+d)%mod*inv%mod)%mod;
            }
        }
    }
    return tmp;
}
signed main(){
    T=read();
    pre_work();
    inv=ksm(2,mod-2);
    int now=0;
    while(T--){
        now++;
        int m=read(),p=read();
        int tmp=((m*work(p-1)%mod-m*(p-1ll)%mod*p%mod*inv%mod)%mod+mod)%mod;
        printf("Case #%lld: %lld
",now,tmp);
    }
    return 0;
}