cdqz2017-test1-數論 （BSGS + 二次剩餘 + CRT）

若m=0，

就是求n^2n ≡ x mod p （x--）

由於必定優解，因此x必定是p的二次剩餘

令g爲p的1個原根，且g^k ≡ x mod p

則k是偶數，證實k是偶數：

假設

g1^k1 ≡ x mod p

g2^k2 ≡ x mod p，k2是偶數

g1^k3 ≡ g2 mod p

那麼 g1^k3k2 ≡ x ≡ g1^k1 mod p

由歐拉定理可得，k3k2 ≡ k1 mod p-1   

∴ k1是偶數

因此對於任意g，k是偶數

因此等價於求 n^n ≡ g^(k/2) mod p

顯然

知足 n≡ g mod  p  且  n ≡ k/2 mod p-1 的n 是一個可行解

又由於p與p-1必定互質，因此用CRT便可求得n

 

若m≠0   

求n^2n+n^m ≡ x mod p

在上面m=0 的時候，咱們是令 n ≡ g mod p

即 n^m ≡ g^m mod p

能解出合法的n的條件是 x-n^m 是p的二次剩餘

因此嘗試枚舉g，判斷x-g^m 是不是p的二次剩餘

判斷方法：利用歐拉準則計算勒讓德符號，即 判斷（x-n^m）^ ((p-1)/2)  mod p 是否等於1

若是x-g^m 是 p的二次剩餘

方程變成 n^2n ≡ x-g^m mod p

令g^k ≡ x-g^m mod p

用BSGS求出一個知足上述條件的k

若k是偶數

那麼方程就變成了 n^n ≡ g^(k/2) mod p

知足 n≡ g mod  p  且  n ≡ k/2 mod p-1 的n 是一個可行解

又由於p與p-1必定互質，因此用CRT便可求得n

 

#include<map>
#include<cmath>
#include<cstdio>
#include<iostream>

using namespace std;

map<int,int>mp;

int Pow(int a,int b,int p)
{
    int res=1;
    for(;b;a=1LL*a*a%p,b>>=1)
        if(b&1) res=1LL*res*a%p;
    return res;
}

long long Mul(long long a,int b,long long p)
{
    long long res=0;
    while(b)
    {
        if(b&1) res+=a,res%=p;
        b>>=1; a+=a; a%=p;
    }
    return res;
}        

bool Legendre(int n,int p)
{
    return Pow(n,p-1>>1,p)+1!=p;
}

int bsgs(int a,int b,int p)
{
    mp.clear();
    int m=sqrt(p);
    mp[b]=0;
    for(int i=1;i<=m;++i)
    {
        b=1LL*a*b%p;
        mp[b]=i;
    }
    int am=Pow(a,m,p);
    int mul=1;
    for(int i=1;i<=m;++i)
    {
        mul=1LL*mul*am%p;
        if(mp.find(mul)!=mp.end()) return i*m-mp[mul];
    }
    return -1;
}

int inv(int a,int p)
{
    return Pow(a,p-2,p);
}

int main()
{
    freopen("theory.in","r",stdin);
    freopen("theory.out","w",stdout);
    int x,m,p;
    scanf("%d%d%d",&x,&m,&p);
    if(p==2) printf("1");
    int y;
    long long ans;
    for(int g=1;;++g)
    {
        if(!Legendre(x-Pow(g,m,p),p)) continue; 
        y=bsgs(g,(x-Pow(g,m,p)+p)%p,p);
        if(y==-1 || (y&1)) continue;
        long long P=1LL*p*(p-1);
        ans=Mul(1LL*(p-1)*inv(p-1,p)%P,g,P)+1LL*p*(y/2)%P;
        ans%=P;
        cout<<ans;
        return 0;
    }
}