题解：P5668 【模板】N 次剩余

· · 题解

Upd：用递推解模 2^q 的情况居然被叉飞了，惊恐地发现本题解幸存。

此题解将使用复杂度正确的方法解决模 2^q 的情况。记 \{0,1,\dots,m-1\} 为 \Z_m。

解关于整数 x 的方程，x^n\equiv k\pmod m，其中 x\in\Z_m。

首先按照惯例，分解 m=\prod p_i^{q_i}，接下来只需要考虑素数的幂。2 的幂比较特殊，没有原根，分开考虑。

回想 m 是素数时如何求解，特判掉 0 的情况，首先找到 m 的一个原根 g，接着求解方程 g^a\equiv k\pmod m。显然 x 一定可以表示为 g^t，得到方程 g^{tn}\equiv g^a\pmod m，两侧同时取离散对数，有 tn\equiv a\pmod{(m-1)}，扩欧求得 t，快速幂得 x 即可。同时知道此方程若有解，则有 \gcd(n,m-1) 个解。

同样的方法可以应用到此题。于是分情况：

首先是 k=0。则有 x^n\bmod p^q=0，满足时当且仅当 x\bmod p^{\lceil\frac qn\rceil}=0。

然后当 k 和 p 不互质就不太行，所以将 k 分解为 u\cdot p^t(u\bmod p\ne0)，重新解 x^n\equiv u\pmod{p^{q-t}} 即可。

接下来是 p_i=2。有 x^n\equiv k\pmod{2^q}。 可以发现每个数可以唯一地表示为 (-1)^a5^b\bmod 2^q，其中 a\in\Z_2,b\in\Z_{2^{q-2}}。两个维度是独立的，于是套用 m 是素数的解法。

最后是 p_i 为奇素数。套用 m 是素数的写法即可。

下面是详细证明。

引理 1（中国剩余定理）\quad 设 m=\prod p_i^{q_i}，p_i 为素数且 q_i>0，则 \Z_m 和 \Z_{p_1^{q_1}}\times\Z_{p_2^{q_2}}\times\cdots 可以建立一个唯一的双射，满足这个双射的一侧的映射 f:\Z_m\to\Z_{p_1^{q_1}}\times\Z_{p_2^{q_2}}\times\cdots，表示为 f(x)=(x\bmod p_1^{q_1},x\bmod p_2^{q_2},\dots)。

证明：显然映射 f 成立，考虑证明 f 映射到的两两不同。使用反证法，假设 x_1\ne x_2，同时 x_1\bmod p_1^{q_1}=x_2\bmod p_1^{q_1},x_1\bmod p_2^{q_2}=x_2\bmod p_2^{q_2},\dots。可以发现 x_1-x_2 必然是所有 p_i^{q_i} 的倍数，又因为 p_i^{q_i} 都是两两互素，所以得到 x_1-x_2 是 m 的倍数，矛盾。证毕。

于是可以分解原方程为若干互相独立的方程，并逐个求解，最后合并。

引理 2 \quad 若关于整数 x 的方程 x^n\equiv k\pmod m 有解，当且仅当 k^{\frac{\varphi(m)}{\gcd(n,\varphi(m))}}\equiv1\pmod m。其中 m 有原根且 \gcd(k,m)=1。（注意，此处 \gcd(0,x)=x。）

证明：

充分性：则须证明对于所有 x，(x^n)^{\frac{\varphi(m)}{\gcd(n,\varphi(m))}}\equiv1\pmod m。显然可化为 x^{\operatorname{lcm}(n,\varphi(m))}\equiv1\pmod m，由欧拉定理可证。

必要性：设 m 的一个原根为 g，显然存在唯一 k'\in\Z_{\varphi(m)} 满足 g^{k'}\equiv k\pmod m，亦存在唯一 x'\in\Z_{\varphi(m)} 满足 g^{x'}\equiv x\pmod m。则条件变为 k'\cdot\frac{\varphi(m)}{\gcd(n,\varphi(m))}\equiv0\pmod{\varphi(m)}，可得 k' 必然是 \gcd(n,\varphi(m)) 的倍数。另一方面，要证明原方程有解，则要证明 x'n\equiv k'\pmod{\varphi(m)} 有解，而此方程有解的条件即为 k' 是 \gcd(n,\varphi(m)) 的倍数（根据 Bézout's Lemma）。得证。

然后就可以判别对于一个 p^q 有没有解了。不幸的是，当 q>2 时，2^q 并没有原根，所以另外处理。先证明几个基础的。

引理 3 \quad 5^{2^{q-3}}\equiv2^{q-1}+1\pmod{2^q}，其中 q>2。

证明：使用归纳法。显然 q=3 时成立。假设 q=i 成立，有 5^{2^{i-3}}\equiv2^{i-1}+1\pmod{2^i}，则 5^{2^{i-3}}\equiv2^{i-1}+1\pmod{2^{i+1}} 或者 5^{2^{i-3}}\equiv2^i+2^{i-1}+1\pmod{2^{i+1}}。两侧同时平方，前一个得到 5^{2^{i-2}}\equiv2^{2i-2}+2^i+1\equiv2^i+1\pmod{2^{i+1}}，后一个得到 5^{2^{i-2}}\equiv2^{2i+1}+2^{2i-2}+2^{i+1}+2^i+1\equiv2^i+1\pmod{2^{i+1}}。于是 q=i+1 成立。根据归纳，证毕。

引理 4 \quad 5 在模 2^q 意义下的阶为 2^{q-2}，其中 q>2。

证明：易知 5 的阶必定是 \varphi(2^q)=2^{q-1} 的因数，也就是 2^i(i\in\Z_q)。而如果 5^{2^i}\equiv1\pmod{2^q}，则 5^{2^j}\equiv1\pmod{2^q}(j\ge i)，于是根据引理 3，5 的阶必然 >2^{q-3}。另一方面，引理 3 两侧同时平方，有 5^{2^{q-2}}\equiv2^{2q-2}+2^q+1\equiv1\pmod{2^q}，所以 5 的阶必然 \le2^{q-2}。综上有 5 的阶为 2^{q-2}，证毕。

引理 5 \quad 关于 a,b 的方程 (-1)^a5^b\equiv x\pmod{2^q}，其中 q>2,a\in\Z_2,b\in\Z_{2^{q-2}},\gcd(x,2)=1，有唯一解。

证明：首先，根据引理 4，5 的阶是 2^{q-2}，所以 \{5^0\bmod2^q,5^1\bmod2^q,\dots,5^{2^{q-2}-1}\bmod2^q\} 中的元素互不相同。证明没有二解，使用反证法，假设 (-1)^{a_1}5^{b_1}\equiv(-1)^{a_2}5^{b_2}\pmod{2^q}，由上得 a_1\ne a_2。于是有 5^{b_1}+5^{b_2}\equiv0\pmod{2^q}，但是左侧模 4 为 2，右侧模 4 为 0，矛盾。从而得证。证明不可能无解，则由上有所有 (-1)^a5^b 互不相同，而一共有 2^{q-1} 个 (-1)^a5^b，于是一定可以和所有 x 对应，证毕。

发现 x=(-1)^a5^b\pmod{2^q} 有唯一解，所以考虑改写方程为 (-1)^{a_xn}5^{b_xn}\equiv(-1)^{a_k}5^{b_k}\pmod{2^q}。a_k 和 b_k 可由两次离散对数得到。然后因为相互独立，拆解为两个方程 a_xn\equiv a_k\pmod2 和 b_xn\equiv b_k\pmod{2^{p-2}}。扩欧求解即可。

至此完整地给出了 x^n\equiv k\pmod{p^q}（p 是素数且 \gcd(k,p)=1）的情况的证明。接下来的问题是：如果 k=0 或是 \gcd(k,p^q)>1 怎么办？

先解决 k=0。

引理 6 \quad 关于整数 x 的方程 x^n\equiv0\pmod{p^q} 的解为 x\equiv0\pmod{p^{\lceil\frac qn\rceil}}。

证明：反证。若 x 含 p 少于 \lceil\frac qn\rceil，则 x^n 含 p 少于 \lceil\frac qn\rceil\cdot n，也就少于 q，矛盾。证毕。

然后是不互素。

引理 7 \quad 关于整数 x 的方程 x^n\equiv u\cdot p^t\pmod{p^q}，其中 \gcd(u,p)=1，有解当且仅当 t\bmod n=0。

证明：直接构造。设 x^n\equiv u\pmod{p^{q-t}} 的解为 x_0，则原方程解为 (x_0+z\cdot p^{q-t})\cdot p^{\frac tn}。证毕。

至此全部情况梳理完毕。时间复杂度为 O(\sqrt m+\tau\log m+s)，s 为解数，\tau 为离散对数时间复杂度。

代码用了下 pb_ds 的哈希表（用来加速 BSGS）。处理 p=2 时马蜂较为阴间，且含有压行，建议不阅读。

#include<stdio.h>
#include<vector>
#include<ext/pb_ds/assoc_container.hpp>
typedef long long ll;
bool f;
int lp,le,e[2000001];
ll n,m,k,p[21],pw[21],c[21];
std::vector<ll>r[21];
__gnu_pbds::gp_hash_table<ll,ll>h;
ll gcd(ll x,ll y){return y?gcd(y,x%y):x;}
ll Pow(ll a,ll b,ll m)
{
    ll s=1;
    for(a%=m;b;a=a*a%m,b>>=1)
        if(b&1)
            s=s*a%m;
    return s;
}
ll BSGS(ll a,ll b,ll m)
{
    a%=m,b%=m,h.clear();
    if(b==1)
        return 0;
    if(!a&&b)
        return -1;
    ll s=sqrt(m)+1,aa=1;
    for(ll i=0;i<s;i++)
        h[b*aa%m]=i,aa=aa*a%m;
    for(ll i=1,j=aa;i<=s;i++,j=j*aa%m)
    {
        int x=h[j];
        if(x||j==b)
            return i*s-x;
    }
    return -1;
}
ll exgcd(ll a,ll b,ll&x,ll&y)
{
    if(!b)
        return x=1,y=0,a;
    int g=exgcd(b,a%b,y,x);
    return y-=a/b*x,g;
}
bool DiscreteRoots(ll z)
{
    if(k%pw[z])
    {
        ll kp=k%pw[z],pk=1,mk=0;
        while(kp%p[z]==0)
            kp/=p[z],pk*=p[z],mk++;
        if(mk%n)
            return 1;
        ll pwk=pw[z]/pk,h=pwk/p[z]*(p[z]-1),g,le=0,e[21];
        if(Pow(kp,h/gcd(n,h),pwk)!=1)
            return 1;
        for(ll i=2;i*i<=h;i++)
            if(h%i==0)
            {
                while(h%i==0)
                    h/=i;
                e[le++]=i;
            }
        if(h>1)
            e[le++]=h;
        for(g=1;;g++)
        {
            bool f=1;
            for(int i=0;i<le&&f;i++)
                if(Pow(g,pwk/p[z]*(p[z]-1)/e[i],pwk)==1) 
                    f=0;
            if(f)
                break;
        }
        ll x=BSGS(g,kp,pwk),ka,kb,y=exgcd(n,pwk/p[z]*(p[z]-1),ka,kb);
        if(x%y)
            return 1;
        h=pwk/p[z]*(p[z]-1)/y,kb=Pow(p[z],mk/n,pw[z]),r[z].resize(y*(pk/kb)),r[z][0]=Pow(g,(ka*(x/y)%h+h)%h,pwk),ka=kb,kb=Pow(g,h,pwk);
        for(ll i=1;i<y;i++)
            r[z][i]=r[z][i-1]*kb%pwk;
        for(ll i=y;i<y*(pk/ka);i++)
            r[z][i]=(r[z][i-y]+pwk)%(pw[z]/ka);
        for(ll i=0;i<y*(pk/ka);i++)
            r[z][i]=r[z][i]*ka;
    }
    else
    {
        ll x=0,y=1,o=Pow(p[z],n,pw[z]);
        while(y)
            x++,y=y*o%pw[z];
        r[z].resize(o=pw[z]/(y=Pow(p[z],x,pw[z]+1)));
        for(ll i=0;i<o;i++)
            r[z][i]=i*y;
    }
    return 0;
}
void dfs(ll st,ll x)
{
    if(st==lp)
    {
        e[le++]=x;
        return;
    }
    for(ll i:r[st])
        dfs(st+1,(x+i*c[st])%m);
    return;
}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld",&n,&m,&k),f=1,lp=le=0;
        if(n==1)
        {
            printf("1\n%lld\n",k);
            continue;
        }
        if(!(m&7))
        {
            ll mp=3;
            p[0]=2,pw[0]=8,m>>=3;
            while(!(m&1))
                m>>=1,pw[0]<<=1,mp++;
            if(k&pw[0]-1)
            {
                ll k2=k&pw[0]-1,r2=0;
                while(!(k2&1))
                    k2>>=1,r2++;
                if(r2%n)
                {
                    puts("0");
                    continue;
                }
                if(pw[0]>>r2>4)
                {
                    ll w_1=0,w5=BSGS(5,k2,pw[0]>>r2);
                    if(!~w5)
                        w_1=1,w5=BSGS(5,(pw[0]>>r2)-k2,pw[0]>>r2);
                    if(n&1)
                    {
                        ll x,y,g=exgcd(n,pw[0]>>r2+2,x,y);
                        if(w5%g)
                        {
                            puts("0");
                            continue;
                        }
                        w5/=g,x=x*w5&(pw[0]>>r2+2)/g-1,y=Pow(5,(pw[0]>>r2+2)/g,pw[0]>>r2),r[0].resize(g),r[0][0]=w_1?(pw[0]>>r2)-Pow(5,x,pw[0]>>r2):Pow(5,x,pw[0]>>r2);
                        for(ll i=1;i<g;i++)
                            r[0][i]=r[0][i-1]*y&(pw[0]>>r2)-1;
                    }
                    else if(w_1)
                    {
                        puts("0");
                        continue;
                    }
                    else
                    {
                        ll x,y,g=exgcd(n,pw[0]>>r2+2,x,y);
                        if(w5%g)
                        {
                            puts("0");
                            continue;
                        }
                        w5/=g,x=x*w5&(pw[0]>>r2+2)/g-1,y=Pow(5,(pw[0]>>r2+2)/g,pw[0]>>r2),r[0].resize(g<<1),r[0][g]=(pw[0]>>r2)-(r[0][0]=Pow(5,x,pw[0]>>r2));
                        for(ll i=1;i<g;i++)
                            r[0][i]=r[0][i-1]*y&(pw[0]>>r2)-1,r[0][i+g]=r[0][i+g-1]*y&(pw[0]>>r2)-1;
                    }
                }
                else
                    for(ll i=0;i<pw[0]>>r2;i++)
                        if(Pow(i,n,pw[0]>>r2)==k&(pw[0]>>r2)-1)
                            r[0].emplace_back(i);
                ll s=r[0].size();
                r[0].resize(s<<r2-r2/n);
                for(ll i=s;i<s<<r2-r2/n;i++)
                    r[0][i]=r[0][i-s]+(pw[0]>>r2)&(pw[0]>>r2/n)-1;
                for(ll i=0;i<s<<r2-r2/n;i++)
                    r[0][i]=r[0][i]<<r2/n;
            }
            else
            {
                ll l=1ll<<mp-(mp-1)/n-1;
                r[0].resize(l);
                for(ll i=0;i<l;i++)
                    r[0][i]=i<<(mp-1)/n+1;
            }
            lp++;
        }
        else if(!(m&3))
        {
            p[0]=2,pw[0]=4,m>>=2;
            for(ll i=0;i<4;i++)
                if(Pow(i,n,4)==(k&3))
                    r[0].emplace_back(i);
            if(r[0].empty())
            {
                puts("0");
                continue;
            }
            lp++;
        }
        else if(!(m&1))
            p[0]=2,pw[0]=2,m>>=1,r[0].emplace_back(k&1),lp++;
        for(ll i=3;i*i<=m;i+=2)
            if(m%i==0)
            {
                p[lp]=i,pw[lp]=1;
                while(m%i==0)
                    m/=i,pw[lp]*=i;
                if(DiscreteRoots(lp))
                {
                    f=0;
                    break;
                }
                lp++;
            }
        if(!f)
        {
            puts("0");
            while(lp--)
                r[lp].clear();
            continue;
        }
        if(m>1)
        {
            p[lp]=pw[lp]=m,m=1;
            if(DiscreteRoots(lp))
            {
                puts("0");
                while(lp--)
                    r[lp].clear();
                continue;
            }
            lp++;
        }
        for(ll i=0;i<lp;i++)
            m*=pw[i];
        for(ll i=0;i<lp;i++)
            c[i]=m/pw[i]*Pow(m/pw[i],pw[i]/p[i]*(p[i]-1)-1,pw[i])%m;
        dfs(0,0),std::sort(e,e+le),printf("%lld\n",le);
        for(ll i=0;i<le;++i,i<le?putchar(' '):putchar('\n'))
            printf("%lld",e[i]);
        while(lp--)
            r[lp].clear();
    }
    return 0;
}