任意模数的K次剩余 /HDU 3930 Broot/BZOJ 2219 数论之神/51nod 1123 X^A mod B

前言

难以置信, 遇到这个问题竟然是因为做一道CTF Reverse的题, 实际上我以前打ACM的时候都没有写过任意模数的K次剩余这种东西. 最简单的模质数的K次剩余的例题是
“`HDU 3930 Broot“`, 进阶版本(模任意奇数的K次剩余)的例题是“`BZOJ 2219 数论之神“`, 而终极版本(任意模数的K次剩余)似乎可以交的题只有“`51nod 1123“`, 截至本文写作时只有44个AC. 然而这个CTF题本身就有50多个通过了. ~~CTF>ACM 实锤~~ (暴论)

其实是因为这个Reverse的题里面幂指数是质数,所以会跟模数互质,然后就可以用RSA的方法来做了. 其次该题目更多的工作在于去花指令混淆上, 而非算法破解.

OK那么就来看一下这个问题从最简单的版本到终极版本怎么求解吧…

模质数的K次剩余

已知$a$, $b$, $p$, 求使得

$$
x^{a} \equiv b \quad(\bmod \ p)
$$

成立的所有 $x$. 其中 $p$ 是质数.

由于 $p$ 是质数, 所以 $p$ 存在原根 $g$ , 此时对于模 $p$ 意义下的任意数 $w$ ($0\le w \le p-1$) 存在唯一的 $i$ ($0\le i \le p-1$) 使得 $w\equiv g^i\quad(\bmod \ p)$.

由此可以将最终的答案用 $g$ 来表示: $x=g^{c}$, 带入上式转化为求解

$$
(g^{c})^{a}\equiv (g^{a})^c\equiv b \quad(\bmod \ p)
$$

此时 $g$, $a$, $b$, $p$, 已知, 只需要解出来 $c$. 此时相当于求解离散对数, 使用 Baby-Step-Giant-Step 可以在 $\mathcal{O}(\sqrt{p})$ 时间内得到一个特解 $x_0 \equiv g^{c} \quad(\bmod \ p)$.

在已知一个特解的情况下求出所有解是简单的, 由原根的性质可知 $g^{\varphi(p)}\equiv 1 \quad(\bmod \ p)$, 因此:

$$
\forall t \in \mathbb{Z}, x^{a} \equiv g^{c \cdot a+t \cdot \varphi(p)} \equiv b \quad(\bmod \ p)
$$

因此所有的解为:

$$
\forall t \in \mathbb{Z}, a \mid t \cdot \varphi(p), x \equiv g^{c+\frac{t \cdot \varphi(p)}{a}} \quad(\bmod p)
$$

上面幂次部分要能整除必须要有$\frac{a}{\operatorname{gcd}(a, \varphi(p))} \mid t$, 可以设 $t=\frac{a}{\operatorname{gcd}(a, \varphi(p))} \cdot i$, 于是得到所有的解为:

$$
\forall i \in \mathbb{Z}, x \equiv g^{c+\frac{\varphi(p)}{\operatorname{gcd}(a, \varphi(p))} \cdot i} \quad(\bmod p)
$$

HDU 3930 Broot 代码:

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 i128;

const int maxn=100;

/*
 *  Complexity:
 *  Get primitive root: O(sqrt(M)+0.25*M*log(M))
 *  Discrete logarithm: O(sqrt(M))
 *  Total: O(sqrt(M)+0.25*M*log(M))
 */

ll gcd(ll a, ll b){
    return b?gcd(b,a%b):a;
}

ll qpow(i128 a, i128 b, i128 p){
    ll res=1;
    while(b){
        if(b&1)
            res=(res*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return res;
}

int pi[maxn],pcnt;
ll get_g(ll P){
    pcnt=0;
    ll phi= P - 1;
    ll tphi=phi;
    if(!(tphi & 1ll)){
        pi[++pcnt]=2;
        while(!(tphi & 1)) tphi>>=1;
    }
    for(ll i=3;i*i <= tphi; i+=2){
        if(tphi % i == 0){
            pi[++pcnt]=i;
            while(tphi % i == 0) tphi/=i;
        }
    }
    if(tphi != 1) pi[++pcnt]=tphi;
    for(ll g=2; g < P; ++g){
        bool flag=true;
        for(int i=1;i<=pcnt;++i){
            if(qpow(g,phi/pi[i], P) == 1){
                flag=false;
                break;
            }
        }
        if(flag) return g;
    }
    return -1;
}

unordered_map<ll,int> Hash;
ll bsgs(ll a, ll b, ll M){
    if(b==1) return 0;
    Hash.clear();
    ll sq=sqrt(M*1.0);
    i128 x=1,p=1;
    for(int i=0;i<sq;++i,p=p*a%M) Hash[p*b%M]=i;
    for(ll i=sq,j;i<=M;i+=sq){
        if(Hash.find(x=x*p%M)!=Hash.end()){
            return i-Hash[x];
        }
    }
    return -1;
}


vector<ll> ans;
void k_residue(ll k, ll y, ll M){
    ans.clear();
    if(y==0){
        ans.push_back(0);return ;
    }
    ll g=get_g(M);
    ll c=bsgs(qpow(g,k,M),y,M);
    if(c==-1) {
        ans.push_back(-1);return;
    }
    ll delta=(M-1)/gcd(k,M-1);
    for(ll cur=c%delta;cur<M-1;cur+=delta)
        ans.push_back(qpow(g,cur,M));
    return ;
}

int main() {
    ll k,y,M;
    int cs=0;
    while(scanf("%lld %lld %lld",&k,&M,&y)==3){
        ++cs;printf("case%d:\n",cs);
        k_residue(k,y,M);
        sort(ans.begin(),ans.end());
        for(int i=0,len=ans.size();i<len;++i){
            printf("%lld\n",ans[i]);
        }
    }
    return 0;
}

模任意奇数的K次剩余

已知$a$, $b$, $p$, 求使得

$$
x^{a} \equiv b \quad(\bmod \ p)
$$

成立的所有 $x$. 其中 $p$ 是任意奇数.

此时套用原来的解法会存在以下两个问题:

  1. $p$ 不一定有原根
  2. 即使 $p$ 有原根, $b$ 跟 $p$ 不互质的话则无法用原根表示 $b$, 也就无法转化为离散对数求解.

对于问题1, 可以考虑将 $p$ 做质因子分解, 即$p=p_{1}^{e_1}…p_{j}^{e_j}$. 由原根的存在性定理可知 $p_{j}^{e_j}$ 一定有原根. 将每个 $p_{j}^{e_j}$ 作为模数分别求解得到 $n_j$ 个解, 求完了之后对于其中的每一个解, 令其与模数为 $p_{1}^{e_1}…p_{j-1}^{e_{j-1}}$ 时的所有解(假设有 $n_{j-1}$ 个)用中国剩余定理两两合并, 这可以用两个数组迭代地实现. 本次合并完后应该得到 $n_{j}\times n_{j-1}$ 个解, 注意每次合并完之后去重一下.

对于问题2, 我们在将问题分解到对每个 $p_{j}^{e_j}$ 上之后, 对 $b \ \bmod(\ p_{j}^{e_j})$ 的结果与 $p_{j}^{e_j}$ 的关系进行分类讨论.

  1. $gcd(b \ \bmod(\ p_{j}^{e_j}), p_{j}^{e_j})=1$

    此时 $b \ \bmod(\ p_{j}^{e_j})$ 可以表示为 $g^c$, 用数论的语言来说就是可以求指标. 此时直接套用模质数的K次剩余即可.

  2. $1\lt gcd(b \ \bmod(\ p_{j}^{e_j}), p_{j}^{e_j}) \lt p_{j}^{e_j}$

    假设 $x^a \ \bmod(\ p_{j}^{e_j})=b \ \bmod(\ p_{j}^{e_j})=u\times p_{j}^{k}$, 其中 $gcd(u,p)=1$. 令同余等式同时除以一个 $p_{j}^{k}$ 得到 $(\frac{x}{p_{j}^{\frac{k}{a}}})^a \ \bmod(\ p_{j}^{e_j-k})=u$ (如果不能除说明无解). 此时模数$p_{j}^{e_j-k}$ 与 $u$ 互质, 回到情况1, 直接套用模质数的K次剩余, 在求出来特解之后, 找所有解的过程中将每个解乘上 $p_{j}^{\frac{k}{a}}$得到原问题的解, 并且乘上之后要对$p_{j}^{k}$ 取模(之前的套用模质数的时候运算都需要对$p_{j}^{e_j-k}$取模, 包括枚举解过程中的增量).

  3. $gcd(b \ \bmod(\ p_{j}^{e_j}), p_{j}^{e_j}) = p_{j}^{e_j}$

    此时其实意味着$b \ \bmod(\ p_{j}^{e_j})=0$. 若$x^a\equiv 0 \ \bmod(\ p_{j}^{e_j})$, 则必有 $p^{\lceil\frac{e_j}{a}\rceil}|x$. 此时只需要枚举其倍数即可, 则解为 $\forall i \in \mathbb{Z}, x \equiv i\times p^{\lceil\frac{e_j}{a}\rceil} \ \bmod(\ p_{j}^{e_j})$

这样就解决了模任意奇数K次剩余的问题, 例题是
“`BZOJ 2219 数论之神“`, 我没有写这个题的代码(BZOJ已经无了), 而是直接写了更进阶的任意模数K次剩余的代码, 见下一小节.

模任意数的K次剩余

已知$a$, $b$, $p$, 求使得

$$
x^{a} \equiv b \quad(\bmod \ p)
$$

成立的所有 $x$. 其中 $p$ 是任意整数.

此时再套用模奇数的K次剩余的方法, 有一个问题是 $p$ 的因子中可能会包括 $2$ , 而我们知道对于 $\forall k \in \mathbb{Z},2^k$ , 只有 $2$ 和 $4$ 是有原根的:

原根存在定理:一个数 $m$ 存在原根当且仅当 $m=2,4, p^{\alpha}, 2 p^{\alpha}$, 其中 $p$ 为奇素数, $\alpha \in \mathbb{N}^{*}$

那么我们在将模数分割到每个 $p_j^{e_j}$ 上的时候就需要对 $p_j=2$ 这种情况做特殊处理. 此时不能用原根将问题转化为离散对数问题了, 前面的分析全部失效, 但当解的数量较小的时候, 可以用倍增的方法来求出所有模 $2^k$ 意义下的解.

具体地, 假设我们要求的最终答案为 $x$, 满足 $x^{a} \equiv b \quad(\bmod \ p)$, 且 $gcd(2^e,p)=2^e, e>0$, 那么如果上式在模 $2^{k-1}$ 意义下存在一个解为 $x_i$, 则它对应于 $2^k$ 的解只能是 $x_i$ 或者 $x_{i}+2^{k-1}$. 我们只需要对于每个$i$把这两个可能的解代进$x^{a} \equiv b \quad(\bmod \ 2^i)$里检验一下就可以了, 如果不能使得式子成立就丢掉. 当然直观来看这种做法似乎会得到最多 $2^{e}$ 个解, 但实际上一般解是很少的, 可以保证较快的算出来.

将模 $2^e$ 意义下的所有解求出之后, 其他的就和模任意质数的K次剩余一样了, 对不同的质因子求一遍然后用中国剩余定理合并即可. 这样的话就解决了任意模数K次剩余的问题.

时间复杂度:

  • 质因子分解 $\mathcal{O}(\sqrt{p})$
  • 找原根 $\mathcal{O}(\sqrt{p}+0.25\times p\times log(p))$
  • 离散对数BSGS $\mathcal{O}(\sqrt{p})$
  • CRT合并, $2^e$ 因子的处理 $\mathcal{O}(答案的个数*log(p))$

例题
“`51nod 1123 X^A mod B“`, 代码:

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 i128;

const int maxn=100;

ll upperInt(ll a, ll b){
    return (a+b-1)/b;
}

ll gcd(ll a, ll b){
    return b?gcd(b,a%b):a;
}

ll extEuclid(ll a,ll b,ll &x,ll &y){
    ll d;
    if(b==0){
        x=1;y=0;
        return a;
    }
    d=extEuclid(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

ll qpow(i128 a, i128 b, i128 p){
    ll res=1;
    while(b){
        if(b&1)
            res=(res*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return res;
}


void factor(ll P, ll pi[], ll ai[], ll &pcnt){
    pcnt=0;
    ll tP=P;
    if(!(tP & 1ll)){
        pi[++pcnt]=2;
        ai[pcnt]=0;
        while(!(tP & 1)) tP>>=1,ai[pcnt]++;
    }
    for(ll i=3;i*i <= tP; i+=2){
        if(tP % i == 0){
            pi[++pcnt]=i;
            ai[pcnt]=0;
            while(tP % i == 0) tP/=i,ai[pcnt]++;
        }
    }
    if(tP != 1) {
        pi[++pcnt]=tP;
        ai[pcnt]=1;//error2
    }
}

ll piG[55],aiG[55],pcntG;
ll getG(ll P, ll phi){
    factor(phi,piG,aiG,pcntG);
    for(ll g=2; g < P; ++g){
        bool flag=true;
        for(int i=1;i<=pcntG;++i){
            if(qpow(g,phi/piG[i], P) == 1){
                flag=false;
                break;
            }
        }
        if(flag) return g;
    }
    return -1;
}

unordered_map<ll,int> Hash;
ll bsgs(ll a, ll b, ll M){
    if(b==1) return 0;
    Hash.clear();
    ll sq=sqrt(M*1.0);
    i128 x=1,p=1;
    for(int i=0;i<sq;++i,p=p*a%M) Hash[p*b%M]=i;
    for(ll i=sq;;i+=sq){
        if(Hash.find(x=x*p%M)!=Hash.end()){
            return i-Hash[x];
        }
        if(i>M) break;//error3
    }
    return -1;
}

ll as[55],ps[55];
ll crt(ll as[], ll ps[], ll len)
{
    int i;
    long long d,x,y,m,n,ret;
    ret=0;
    n=1;
    for (i=0;i<len;i++)
        n*=ps[i];
    for (i=0;i<len;i++){
        m= n / ps[i];
        d=extEuclid(ps[i], m, x, y);
        ret=(ret+y*m*as[i])%n;
    }
    return (n+ret%n)%n;
}

ll gs[55],Ps[55],ss[55],piR[55],aiR[55],pcntR;
ll ans[100005],nans,tans[100005],tnans;// error1
void merge(ll sol,ll prevPi, ll curPi){
    as[0] = sol, ps[0] = curPi;
    for (int j = 1; j <= nans; ++j) {
        as[1] = ans[j], ps[1] = prevPi;
        tans[++tnans] = crt(as,ps,2);
    }
}


void kResidue(ll k, ll y, ll P){
    nans=0;

    ll pi2_ai=1;
    if(!(P & 1ll)){
        pi2_ai=2;
        P>>=1;
        if(y&1) ans[++nans]=1;
        else ans[++nans]=0;

        while(!(P&1ll)){
            pi2_ai*=2;
            tnans=0;
            for(int i=1;i<=nans;++i){
                if(qpow(ans[i], k, pi2_ai) == y % pi2_ai) tans[++tnans]=ans[i];
                if(qpow(ans[i]+(pi2_ai >> 1), k, pi2_ai) == y % pi2_ai) tans[++tnans]= ans[i] + (pi2_ai >> 1);
            }
            for(int i=1;i<=tnans;++i)
                ans[i]=tans[i];
            if(tnans==0) {
                ans[nans = 1] = -1;return; //error5
            }
            nans=tnans;
            P>>=1;
        }

    }
    factor(P, piR, aiR, pcntR);

    ll privPi=pi2_ai;//error
    for(int i=1; i <= pcntR; ++i){
        tnans = 0;
        ll pi_ai=pow(piR[i],aiR[i]);
        ll divisor=gcd(y%pi_ai,pi_ai);
        if(divisor==pi_ai){
            ll t = upperInt(aiR[i],k);
            ll base=pow(piR[i],t);
            for(ll sol=base;sol<=pi_ai;sol+=base){// error6 base==pi_ai, error7 +=base
                if(i==1 && pi2_ai == 1){
                    tans[++tnans]=sol;
                }else {
                    merge(sol,privPi,pi_ai);
                }
            }
        }else{
            ll phi=pi_ai/piR[i]*(piR[i]-1);
            ll upperphi=phi;
            ll yy=y%pi_ai;
            ll tpi_ai=pi_ai;
            ll pek=1;
            if(divisor!=1){
                int tmp=0;
                while((divisor/pek)%piR[i]==0){
                    pek*=piR[i];
                    tmp++;
                }
                if(tmp/k*k!=tmp){
                    ans[nans=1]=-1;return ;
                }
                tpi_ai/=pek;
                phi=tpi_ai/piR[i]*(piR[i]-1);
                yy/=pek;
                pek=qpow(piR[i],tmp/k,pi_ai);
            }
            ll g=getG(tpi_ai,phi);
            ll c=bsgs(qpow(g,k,tpi_ai),yy,tpi_ai);
            if(c==-1){
                ans[nans=1]=-1;return ;
            }
            ll delta= phi / gcd(k, phi);
            for(ll cur=c%delta,sol; cur < upperphi; cur+=delta){
                sol=(qpow(g, cur, pi_ai)*pek)%pi_ai;
                if(i==1 && pi2_ai == 1){
                    tans[++tnans]=sol;
                }else{
                    merge(sol,privPi,pi_ai);
                }
            }
        }
        sort(tans+1,tans+1+tnans);
        tnans=unique(tans+1,tans+1+tnans)-tans-1;
        for(int j=1;j<=tnans;++j){
            ans[j]=tans[j];
        }
        nans=tnans;
        privPi*=pi_ai;//error4
    }

    sort(ans+1,ans+1+nans);
    nans=unique(ans+1,ans+1+nans)-ans-1;

    return ;
}

int main() {
//    int v7[4]; // [rsp+10h] [rbp-50h]
//    int v8[4]; // [rsp+20h] [rbp-40h]
//    int v9[4]; // [rsp+30h] [rbp-30h]
//    v7[0] = -855924449;
//    v7[1] = 2093860801;
//    v7[2] = 767425486;
//    v7[3] = 480046732;
//    v8[0] = 873063593;
//    v8[1] = 947071579;
//    v8[2] = 1040105093;
//    v8[3] = 780771613;
//    v9[0] = -785653175;
//    v9[1] = -1371619377;
//    v9[2] = -1173543495;
//    v9[3] = -1099902543;
//    for(int _=0;_<4;++_){
//        kResidue((unsigned int)v8[_],(unsigned int)v7[_],(unsigned int)v9[_]);
//        if(ans[1]==-1){
//            puts("No Solution");
//            continue;
//        }
//        for(int i=1;i<=nans;++i){
//            printf("%lld ",ans[i]);
//        }
//        puts("");
//    }


    ll k,y,M;
    int cs=0;
    scanf("%d",&cs);
    for(int _=0;_<cs;++_){
        scanf("%lld %lld %lld",&k,&M,&y);
        kResidue(k,y,M);
        if(ans[1]==-1){
            puts("No Solution");
            continue;
        }
        for(int i=1;i<=nans;++i){
            printf("%lld ",ans[i]);
        }
        puts("");
    }
    return 0;
}

/*
1
90 100 0

1
3 36 17

Input
1
3 531441 330750
Output
264 19947 39630 59313 78996 98679 118362 138045 157728 177411 197094 216777 236460 256143 275826 295509 315192 334875 354558 374241 393924 413607 433290 452973 472656 492339 512022
 */

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据