N次剩余

Determinant · 2024-09-12 20:29:18 · 题解

本文将介绍一种不求离散对数求任意模数 n 次剩余的方法。具体地，设质因数分解 m 的复杂度为 \operatorname{factorize}(m)，答案个数为 s，且假设值域范围内整数四则运算复杂度为 O(1)，则可以在 O(\operatorname{factorize}(m)+s\log^2m) 的时间内得到答案。

初步转化

根据维基百科，求任意模数二次剩余等价于对模数进行质因数分解。因此本题不弱于质因数分解，不妨先做个质因数分解，转化为求 k^{1/n}\bmod p^c 这样的问题，然后用中国剩余定理合并。显然合并可以在 O(s\log m) 时间内完成。

这里为了避免一些问题，我们先将 n 对 \varphi(p^c) 取模，k 对 p^c 取模，然后讨论一些 corner case（若符合某一条，则执行该条并忽略后面的所有内容）：

(1)若 k=0，答案为 \{0\}。

现在 k\neq0。

(2)若 n=0，k=1，答案为 [1,p-1]\cap\Z。

现在 n\neq0,k\neq0。

(3)记 v_p(n)=\max\{k|n\bmod p^k=0\}，令 h=v_p(k)。当 h\bmod n\neq0 时无解；否则先求 (k/p^h)^{1/n}\bmod p^{c-h} 的答案，然后在 p 进制意义下在每个值的右边补上 h/n 位，全部填 0，然后在左边补上 h(1-1/n) 位，随便填数，将所有值的集合并起来即可。

注记：不难注意到 v_p(k^{1/n})=v_p(k)/n，因此判定无解和在右边补 0 是合理的；展开 (a_0+a_1p+a_2p^2+...)^n，发现最高的 h(1-1/n) 位对结果完全没有贡献，于是可以随便填。至于中间的 c-h 位，可以证明两个问题的解是一一对应的。

现在 n\neq0,k\bmod p\neq0。

(4)令 q=\varphi(p^c)/(1+[p=2][c\geq3])，则有解的充要条件为 k^{q/\gcd(n,q)}\bmod p^c=1。

注记：现在 k 一定是模 p^c 乘法群中的元素。当 p\neq2 或 c\leq2，模 p^c 乘法群同构于循环群 \mathbf{C_{\varphi(p^c)}}，于是 k 可以表示为 g^a 的形式，g 是原根；由于 (n/\gcd(n,q))^{-1}\bmod q 一定存在，这部分不用考虑了；为使 a/\gcd(n,q)\bmod q 存在，a\bmod\gcd(n,q) 必须等于 0，于是 k^{q/\gcd(n,q)}\bmod p^c=1。反过来，满足这个条件，也可以推出一个根是 g^a，于是一定有解。当 p=2 且 c\geq3，乘法群同构于 \mathbf{C_2}\times\mathbf{C_{2^{n-2}}}，可以设 k=g_1^{a_1}g_2^{a_2}，进行类似的分析。

现在 n\neq0,k\bmod p\neq0，且一定有解。

模质数的情形

模 p^c 的情形还是太过复杂，我们先考虑模 p 意义下的根。

不妨将 n 对 p-1 取余，k 对 p 取余。现在有 k\neq0，但可能出现 n=0 的情况。为了方便，我们仍然先讨论一些 corner case。

(5)当 n=0，答案为 [1,p-1]\cap\Z。

(6)当 n=1，答案为 \{k\}。

现在 n\geq2,k\neq0。继续规约：看起来，令 g=\gcd(n,p-1)，我们可以先求出 q=(n/g)^{-1}\bmod(p-1)，然后转化为求 k^q 的 g 次方根，这样就规约到了 n|(p-1) 的情况。但还有一个问题：n/g 的逆元一定存在吗？取 n=16,p-1=18，我们发现实际上并不一定。

这里的解决方法是把 n 加上 p-1 的随机倍数。令 p-1 任一质因数为 q，则要使 \gcd(n+c(p-1),p-1)=1，若 v_q(n)\neq v_q(p-1)，则只需 c\bmod q\neq0；否则，若 v_q(n)=v_q(p-1)，也只要 (n/q^v)+c((p-1)/q^v)\bmod q\neq0，即存在一个 d 使得 c\bmod q\neq d 即可。于是这样的 c 是大量存在的，直接随机就可以快速找到。实际上，本题数据中直接取 1 就能过。

现在 n|(p-1)。求出最大的 t 使得 n^t|(p-1)，令 d=(p-1)/n^t。对 n 做质因数分解，随机一个数求 d 次方，期望随机 o(\log m) 次就能找到一个 n^t 次单位根 g。现在求出 g^0,g^{n^{t-1}},g^{2n^{t-1}},...,g^{(n-1)n^{t-1}} 放入哈希表。

令 D=d^{-1}\bmod n，N=(Dd-1)/n，K=k^D。由于一定有解，必然有 K^{(p-1)/n}\bmod p=1。两边同时开 n 次方，得到一个 n 次单位根，查表，设其为 g^{a_0n^{t-1}}，得 K^{(p-1)/n^2}g^{-a_0n^{t-1}}\bmod p=1。

继续这样两边同时开根然后除过去的操作，最终 K^dg^{-(a_0n+a_1n^2+...+a_{t-2}n^{t-1})}\bmod p=1，也即 k^{Nn}g^{-(a_0n+a_1n^2+...+a_{t-2}n^{t-1})}\bmod p=k^{-1}。

两边开根，k^{1/n}=k^Ng^{a_0+a_1n+...+a_{t-2}n^{t-2}}，于是右边就是一个合法的根；而 n 次单位根已经全部在上面的表中了，可以直接得到所有的根。分析复杂度，每次开根需要计算 O(\log m) 次快速幂，因此需要 O(\log^2 m) 的时间。还有一个小问题：d^{-1}\bmod n 不存在怎么办？当 n 为质数时一定不会出现这种情况，因此把 n 分解质因数之后多次开质数次方根即可，这不影响复杂度为 O(s\log^2m)。

牛顿迭代

现在我们已经求出了模 p 意义下的根，我们要做的就是求出模 p^c 意义下的根。

正如经典的求模 2^k 逆元一样，在这里可以使用牛顿迭代。

定义 \mathbb Q 上的 p-adic 赋值：|\frac ab|=p^{v_p(b)-v_p(a)},a,b\in\Z-\{0\}，特别地，|0|=0。

可以验证 p-adic 赋值满足如下三条性质：

（1）正定性：|x|\geq0,|x|=0\iff x=0

（2）积性：|xy|=|x||y|

（3）强三角不等式：|x+y|\leq\max\{|x|,|y|\}

对于解析函数 f(z)，设 f(\zeta)=0,f'(\zeta)\neq0，定义 \gamma=\sup\limits_{k=2}^\infty\left |\dfrac{f^{(k)}(\zeta)}{k!f'(\zeta)}\right |^{\frac1{k-1}}.

引理 若 \gamma|z-\zeta|<1，则 |f'(z)/f'(\zeta)|=1.

证明 在 \zeta 处展开 f：

而

于是 |\sum\limits_{k=1}^\infty\dfrac{f^{(k+1)}(\zeta)(z-\zeta)^k}{k!f'(\zeta)}|<1，|f'(z)/f'(\zeta)|=1.

定理 若 \gamma|z-\zeta|<1，则 |z-f(z)/f'(z)-\zeta|\leq\gamma|z-\zeta|^2.

证明 在 \zeta 处展开 f：

由引理知 |f'(\zeta)/f'(z)|=1，而 \left|\dfrac{(k-1)f^{(k)}(\zeta)(z-\zeta)^k}{k!f'(\zeta)}\right|\leq|k-1|\gamma^{k-1}|z-\zeta|^k\leq\gamma|z-\zeta|^2.

于是 |z-f(z)/f'(z)-\zeta|\leq\gamma|z-\zeta|^2.

本题中对 f(z)=z^n-k，有 \gamma\geq|1/f'(\zeta)|\geq|1/n|。也就是说，只要求出了模 p^{v_p(n)+1} 意义下的根，就能通过 O(\log\log m) 次牛顿迭代提升到模 p 的任意幂次。那么模 p^{v_p(n)+1} 意义下的根怎么求呢？不妨直接枚举。

由二项式定理可知，k 在 p 进制下最高的 v_p(n) 位对 k^n\bmod p^c 无任何贡献；而另一方面，由前面对乘法群结构的分析，根的个数恰有 p^{v_p(n)}\gcd(n,p-1) 个。于是对每个模 p 意义下的根，我们枚举其 p 进制下第 2 至 v_p(n)+1 位怎么填，做牛顿迭代。应当恰好有一种填法能成功迭代得到一个特解，我们在它 p 进制下最高的 v_p(n) 位随意填数，便得到了所有解。

另外，对于 p=2 且 c\geq3 且 2|n 的情况，需要特殊处理。由于乘法群结构不同，此时根的个数为 2^{\min\{v_2(n)+1,c-1\}}，意味着如果按照刚才的方法，会漏掉一半的根。这里我们注意到，若 x 是方程的根，则 -x 也一定是方程的根，于是在解集中加入所有已知根的相反数即可。

最终，我们得到了 O(\operatorname{factorize}(m)+s\log^2m) 求解任意模数 n 次剩余的方法。实际效率上，对于我的实现，在洛谷模板题比前最优解快了三分之一；理论上，在限制答案个数不太多的前提下，依据 LOJ 上离散对数和分解质因数模板题的最快记录，原先大概能做 10^{20}；通过这个做法规避离散对数之后，大概可以做到 10^{30}。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;
mt19937 rnd(114514);
struct FactorList{
    map<int,int>mp;
    FactorList(){}
    FactorList(int x){mp[x]=1;}
    void operator +(const int &h){mp[h]++;}
    void operator +=(const FactorList &h){
        for(auto it:h.mp)mp[it.first]+=it.second;
    };
    void print(){int fl=0;for(auto it:mp){if(fl)putchar('*');fl=1;printf("%d^%d",it.first,it.second);}puts("");}
};
namespace Factorize{
    const int N=1000;
    int m,c[31],pr[3]={2,7,61};FactorList f[N];
    void init(){
        for(int i=2;i<N;++i){
            int fl=1;
            for(int j=2;j*j<=i;++j)if(i%j==0){fl=0;f[i]=f[i/j];f[i]+=j;break;}
            if(fl)f[i]=i;
        }
    }
    bool chk(int n,int x){
        int v=pr[x],p=n;n--;m=1;while(!(n&1))n>>=1,++m;c[m]=1;
        while(n){if(n&1)c[m]=1ll*c[m]*v%p;v=1ll*v*v%p;n>>=1;}
        for(int i=m-1;i;--i)c[i]=1ll*c[i+1]*c[i+1]%p;
        for(int i=1;i<=m;++i){
            if(c[i]!=1&&c[i]!=p-1)return 0;
            if(i==1&&c[i]==p-1)return 0;
            if(i==m||c[i]==p-1)return 1;
        }
    }
    bool chk(int n){
        for(int i=0;i<3;++i)if(!chk(n,i))return 0;return 1;
    }
    int getd(int n){
        int s=0,t=0,c=rnd()&((1<<30)-1);
        while(1){
            int val=1;for(int i=1;;i<<=1){
                for(int j=1;j<=i;++j){
                    t=(1ll*t*t+c)%n;val=1ll*val*abs(t-s)%n;
                    if(!(i%31)){ll q=__gcd(n,val);if(q>1)return q;}
                }
                int q=__gcd(n,val);if(q==n)return getd(n);if(q>1)return q;s=t;
            }
        }
    }
    FactorList sol(int n){
        if(n<N)return f[n];
        if(chk(n))return n;
        int d=getd(n);FactorList h=sol(d);h+=sol(n/d);return h;
    }
};
void exgcd(int a,int b,int &x,int &y){
    if(!b){x=1;y=0;return;}
    exgcd(b,a%b,y,x);y-=a/b*x;
}
int inv(int x,int p){int a,b;if(x<0)x+=p;exgcd(x,p,a,b);if(a<0)a+=p;return a;}
ll fpowll(ll x,int y,ll p){
    ll s=1;while(y){if(y&1)s=(__int128)s*x%p;x=(__int128)x*x%p;y>>=1;}return s;
}
int fpow(int x,int y,int p){
    int s=1;while(y){if(y&1)s=1ll*s*x%p;x=1ll*x*x%p;y>>=1;}return s;
}
int fpow(int x,int y){
    int s=1;while(y){if(y&1)s=s*x;x=x*x;y>>=1;}return s;
}
vector<int> enumerate(int l,int r,int d){
    vector<int>v;for(int i=l;i<r;i+=d)v.push_back(i);
    return v;
}
vector<int> getroot0(int n,int k,int p){
    int t=0,d=p-1,w,u,npow[32]={1};while(d%n==0)d/=n,t++;for(int i=1;i<t;++i)npow[i]=npow[i-1]*n;
    FactorList N=Factorize::sol(n);unordered_map<int,int>mp;
    if(__gcd(n,d)>1){
        vector<int>w0,w1={k},w2;
        for(auto i:N.mp){
            for(int j=1;j<=i.second;++j){
                w2.clear();
                for(int k:w1){
                    w0=getroot0(i.first,k,p);
                    for(int l:w0)w2.push_back(l);
                }
                swap(w1,w2);
            }
        }
        return w1;
    }
    while(1){
        w=fpow(u=rnd(),npow[t-1]*d,p);if(w<0)w+=p;if(w==0||w==1)continue;
        int fl=1;for(auto it:N.mp)if(fpow(w,n/it.first,p)==1){fl=0;break;}if(fl)break;
    }
    for(int i=0,c=1;i<n;++i)mp[c]=i,c=1ll*c*w%p;u=fpow(u,d,p);
    int ans=0,d_=max(inv(n,d),1),k_=fpow(k,d_*n-1,p);
    for(int i=t-2;i>=0;--i){
        int h=1ll*fpow(k_,npow[i],p)*fpow(u,ans,p)%p;
        ans+=npow[t-1]*(n-mp[h]);ans/=n;
    }
    k=1ll*fpow(k,d_,p)*fpow(u,ans,p)%p;
    vector<int>v;for(int i=0;i<n;++i)v.push_back(k),k=1ll*k*w%p;return v;
}
vector<int> getroot(int n,int k,int p){
    if(!k)return {0};
    if(!n){if(k==1)return enumerate(1,p,1);return {};}
    if(n==1)return {k};
    if(!k)return {0};n+=p-1;
    int g=__gcd(n,p-1);
    k=fpow(k,inv(n,p-1),p);n=g;if(n==1)return {k};
    if(fpow(k,(p-1)/n,p)!=1)return {};
    return getroot0(n,k,p);
}
vector<int> getroot(int n,int k,int p,int pn){
    int q=fpow(p,pn),fac=1,fac1=1,fac2=q;k%=q;
    if(!k)return enumerate(0,q,fpow(p,(pn+n-1)/n));
    int h=0,u=p;while(k%p==0)k/=p,q/=p,h++;
    if(h%n)return {};
    fac=fpow(p,h/n);fac1=fpow(p,pn-h+h/n);pn-=h;
    vector<int>v=getroot(n%(p-1),k%p,p);vector<int>V;
    int n0=n;while(n%p==0)u*=p,n/=p,fac1/=p;u=min(u,q);fac1=max(fac1,fac*p);
    ll q0=1ll*u/p;
    for(int i:v)for(int j=i;j<u;j+=p){
        int x=j,x1,fl=0;
        for(int z=1;z<7;++z){
            ll t0=fpowll(x,n0-1,q0*q),t=(lll)t0*x%(q0*q);
            if((t-k)%q==0){fl=1;break;}
            x1=((k-t)/q0*inv(t0%q*n%q,q)+x)%q;if((x1-x)%u)break;x=x1;
        }
        if(!fl)continue;
        if(x<q)x+=q;x*=fac;x%=fac1;fl=0;for(int j=x;j<fac2;j+=fac1){V.push_back(j);if(j+x==fac2)fl=1;}
        if(p==2&&!fl&&q0>1)for(int j=x;j<fac2;j+=fac1)V.push_back(fac2-j);
        break;
    }
    return V;
}
bool chkroot(int n,int k,int p,int pn){
    int h=0,u=p,q=fpow(p,pn),ord;k%=q;if(!k)return 1;while(k%p==0)k/=p,q/=p,h++;if(h%n)return 0;if(q==1)return 1;
    ord=1ll*(p-1)*q/p;n%=ord;if(!n)return (k==1);
    if(p==2&&q<=4){for(int i=0;i<q;++i)if(fpow(i,n,q)==k)return 1;return 0;}
    if(p==2)ord>>=1;n=__gcd(n,ord);return (fpow(k,ord,p)==1);
}
void sol(){
    int n,m,k;FactorList M;vector<int>w0,w1,w2;w1.push_back(0);
    scanf("%d%d%d",&n,&m,&k);M=Factorize::sol(m);
    for(auto it:M.mp)if(!chkroot(n,k,it.first,it.second)){puts("0");return;}
    for(auto it:M.mp){
        int h=fpow(it.first,it.second);h=(1ll*inv(m/h,h)*m/h)%m;
        w0=getroot(n,k,it.first,it.second);
        if(w0.empty()){w1={};break;}
        w2.resize(w1.size()*w0.size());
        int k=0;for(int i:w1)for(int j:w0)w2[k++]=((1ll*j*h+i)%m);
        swap(w1,w2);
    }
    printf("%d\n",(int)w1.size());sort(w1.begin(),w1.end());for(int i=0;i<w1.size();++i){printf("%d",w1[i]);if(i<w1.size()-1)putchar(' ');else puts("");}
}
int main(){
    Factorize::init();
    int t;scanf("%d",&t);while(t--)sol();
    return 0;
}

闲话

写这篇文章，可以说是为了一瓶醋包了一碗饺子。当时在搜 p-adic 牛顿迭代有没有什么理论依据，找到了 On a p-adic newton method，看了看分析还蛮有意思的，然后就想着能不能有什么应用，就搞出来了这个东西。

在这篇题解中，我们只用了原论文中结论的一个特例。我们只用到 p-adic 赋值的正定性、积性和强三角不等式，实际上这正是非阿基米德赋值的三条性质。于是它可以推广到一般非阿基米德赋值域上。例如，在有理分式域上，类似地，定义 |f(x)/g(x)|=2^{-c}，c 为 f(x)/g(x) 的最低非零项幂次数，便立即可以得到 oi 范围内所有多项式牛顿迭代的合理性。对于阿基米德赋值域(如赋予通常绝对值的实数域)也可以有类似的分析，但由于没有强三角不等式，需要类似等比数列求和的计算(具体过程可参见 Complexity and Real Computation Ch.8)，最终可知当 \gamma|z-\zeta|<\frac{5-\sqrt{17}}4 时二阶收敛。