博客

Pohlig-Hellman algorithm　解决一类离散对数问题

  • Post author:
  • Post category:其他

原理：

https://blog.csdn.net/qq_26060285/article/details/90322054

适用范围大致为模数Ｐ为质数，且(P-1)的最大质因子 <= 1e12。Ｐ为合数的情况直接用EXBSGS就好了，因为P为合数的时候大部分情况下求不了原根, 若P有原根，则 P = 2,4,q^a,2q^a(q为质数)。模板如下：

洛谷3846

#include"iostream"
#include"cstdio"
#include"vector"
#include"unordered_map"
#include"map"
#include"cstring"
#include"algorithm"
#include"cmath"
using namespace std;

namespace DLP{
    typedef long long LL;
    typedef __int128 i128;
    const int MX = 1e6+7;
    LL qpow(LL a, LL n){
        LL ret = 1;
        while(n > 0){
            if(n&1) ret = ret*a;
            a = a*a;
            n >>= 1;
        }
        return ret;
    }

    LL qpow(LL a, LL n, LL p){
        LL ret = 1;
        a %= p;
        n %= p;
        while(n > 0){
            if(n&1) ret = (i128)ret*a%p;
            a = (i128)a*a%p;
            n >>= 1;
        }
        return ret;
    }

    int prime[MX],tt;
    bool is_prime[MX];
    struct node{
        LL p;
        int c; 
    };
    void prime_init()
    {
        memset(is_prime,1,sizeof(is_prime));
        int n = MX - 7;
        for(int i = 2; i <= n; i++){
            if(is_prime[i]) prime[++tt] = i;
            for(int j = 1; j <= tt && prime[j] <= n/i; j++){
                is_prime[i*prime[j]] = 0;
                if(!i%prime[j]) break;
            }
        }
    }

    void find_fra(vector<node> &v, LL num){
        if(tt == 0) prime_init();
        for(int i = 1; i <= tt && prime[i] <= num; i++){
            if(num%prime[i] == 0){
                int cnt = 0;
                while(num%prime[i] == 0){
                    ++cnt;
                    num /= prime[i];        
                }
                v.push_back(node{prime[i],cnt});
            }
        }
        if(num > 1) v.push_back(node{num,1});
    }

    int getroot(LL p, LL phi, vector<node> &v) {
        for(int k=2;; k++) {
            int flag=1;
            for(int i=0; i<(int)v.size(); i++)
                if(qpow(k,phi/v[i].p,p)==1) {
                    flag=0;
                    break;
                }
            if(flag)return k;
        }
    }


    LL BSGS(LL a, LL b, LL p, LL mod) {
        a %= mod, b %= mod;
        //特殊情况，具体看题目要怎么做
        if(b == 1) return 0;
        if(a==0){
            if(b==0) return 1;
            else return -1;
        }

        LL t = 1;       
        int m = int(sqrt(1.0 * p) + 1);
        LL base = b;
        unordered_map<LL,LL> has; 
        for(int i = 0; i < m; ++i) {
            has[base] = i;
            base = (i128)base * a % mod;
        }
        base = qpow(a, m, mod);
        LL now = t;
        //对于A^x，从前往后，保证答案是从小到大的
        for(int i = 1; i <= m + 1; ++i) {
            now = (i128)now * base % mod;
            if(has.count(now)) return i * m - has[now];//找到最小后，找到最小的x使得a^x=1 mod p,通解为 最小+k*x
        }
        return -1;
    }

    LL get_xi(LL g, LL h, LL p, LL c, LL N, LL mod){
        vector<LL> pi;
        LL tp = 1;
        for(int i = 0; i <= c; i++){
            pi.push_back(tp);
            tp *= p;
        }
        LL gq = qpow(g,pi[c-1],mod);
        LL inv = 0;
        tp = 1;
        for(int i = c-1; i >= 0; i--){
            LL tx = tp*BSGS(gq, qpow( (i128)h*qpow(g,pi[c]-inv,mod)%mod, pi[i], mod), p, mod);
            inv += tx;
            tp *= p;
        }
        return inv;
    }
    LL exgcd(LL a,LL b,LL &x,LL &y) {
        if(b==0) {
            x=1,y=0;
            return a;
        }
        LL d=exgcd(b,a%b,y,x);
        y-=a/b*x;
        return d;
    }
    LL Inv(LL a,LL P){//求a在膜P下的逆
        if(a==0)return 1;
        LL x,y;
        exgcd(a,P,x,y);
        return (x%P+P)%P;
    }
    //ax+by=c
    LL solve_exgcd(LL a,LL b,LL c,LL &x,LL &y) {
        LL d;
        if(c%(d=__gcd(a,b)))return -1;//返回-1是无解
        a/=d;
        b/=d;
        c/=d;
        exgcd(a,b,x,y);
        x=(i128)x*c%b;
        while(x<=0) x+=b;
        return x;
    }

    LL CRT(vector<LL> xi, vector<node> v){
        LL sz = v.size();
        LL M = 1, ans = 0;
        vector<LL> m;

        for(int i = 0; i<sz; i++){
            m.push_back(qpow(v[i].p,v[i].c));
            M *= m[i];
        }
        for(int i=0; i<sz; i++){
            LL Mi = M/m[i];
            ans = ((i128)ans + (i128)Mi *Inv(Mi,m[i]) *xi[i])%M;
        }
        if(ans < 0) ans += M;
        return ans;
    }

    LL get_x(LL h, LL g, LL N, LL mod, vector<node> &v){
        vector<LL> xi;
        for(node tp : v){
            LL tg = qpow(g,N/qpow(tp.p,tp.c),mod);
            LL th = qpow(h,N/qpow(tp.p,tp.c),mod);
            xi.push_back(get_xi(tg, th, tp.p, tp.c, N, mod));
        }
        return CRT(xi,v);
    }

    //p为质数且(p-1)分解质因子较大的情况时，算法比较优秀，如果ｐ=2*q+1,q为质数, q<=1e12皆可做，则算法与BSGS差不多，不为质数的话可以考虑EXBSGS
    LL solve(LL g, LL h, LL p){
        if(h == 1) return 0;
        LL phi_p = p-1;     
        vector<node> v;
        find_fra(v,phi_p);
        int rt = getroot(p,phi_p,v);
        LL x = get_x(g,rt,phi_p,p,v);
        LL y = get_x(h,rt,phi_p,p,v);
        LL a = 0, b = 0;
        if(x==0){
            if(y==0) return 1;
            else if(y == 1) return 0;
            else return -1; //不满足返回-1
        }
        else return solve_exgcd(x,phi_p,y,a,b);
    }
}

typedef long long LL;
int main(){
#ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
#endif
    LL g,h,p;
    while(~scanf("%lld%lld%lld",&p,&g,&h)){
        LL ans = DLP::solve(g,h,p);
        if(ans != -1) printf("%lld\n",ans);
        else printf("no solution\n");
    }
    return 0;
}

hdu6632

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

namespace DLP{
    typedef long long LL;
    typedef __int128 i128;
    const int MX = 1e6+7;
    LL qpow(LL a, LL n){
        LL ret = 1;
        while(n > 0){
            if(n&1) ret = ret*a;
            a = a*a;
            n >>= 1;
        }
        return ret;
    }

    LL qpow(LL a, LL n, LL p){
        LL ret = 1;
        a %= p;
        n %= p;
        while(n > 0){
            if(n&1) ret = (i128)ret*a%p;
            a = (i128)a*a%p;
            n >>= 1;
        }
        return ret;
    }

    int prime[MX],tt;
    bool is_prime[MX];
    struct node{
        LL p;
        int c; 
    };
    void prime_init()
    {
        memset(is_prime,1,sizeof(is_prime));
        int n = MX - 7;
        for(int i = 2; i <= n; i++){
            if(is_prime[i]) prime[++tt] = i;
            for(int j = 1; j <= tt && prime[j] <= n/i; j++){
                is_prime[i*prime[j]] = 0;
                if(!i%prime[j]) break;
            }
        }
    }

    void find_fra(vector<node> &v, LL num){
        if(tt == 0) prime_init();
        for(int i = 1; i <= tt && prime[i] <= num; i++){
            if(num%prime[i] == 0){
                int cnt = 0;
                while(num%prime[i] == 0){
                    ++cnt;
                    num /= prime[i];        
                }
                v.push_back(node{prime[i],cnt});
            }
        }
        if(num > 1) v.push_back(node{num,1});
    }

    int getroot(LL p, LL phi, vector<node> &v) {
        for(int k=2;; k++) {
            int flag=1;
            for(int i=0; i<(int)v.size(); i++)
                if(qpow(k,phi/v[i].p,p)==1) {
                    flag=0;
                    break;
                }
            if(flag)return k;
        }
    }


    LL BSGS(LL a, LL b, LL p, LL mod) {
        a %= mod, b %= mod;
        //特殊情况，具体看题目要怎么做
        if(b == 1) return 0;
        if(a==0){
            if(b==0) return 1;
            else return -1;
        }

        LL t = 1;       
        int m = int(sqrt(1.0 * p) + 1);
        LL base = b;
        unordered_map<LL,LL> has; 
        for(int i = 0; i < m; ++i) {
            has[base] = i;
            base = (i128)base * a % mod;
        }
        base = qpow(a, m, mod);
        LL now = t;
        //对于A^x，从前往后，保证答案是从小到大的
        for(int i = 1; i <= m + 1; ++i) {
            now = (i128)now * base % mod;
            if(has.count(now)) return i * m - has[now];//找到最小后，找到最小的x使得a^x=1 mod p,通解为 最小+k*x
        }
        return -1;
    }

    LL get_xi(LL g, LL h, LL p, LL c, LL N, LL mod){
        vector<LL> pi;
        LL tp = 1;
        for(int i = 0; i <= c; i++){
            pi.push_back(tp);
            tp *= p;
        }
        LL gq = qpow(g,pi[c-1],mod);
        LL inv = 0;
        tp = 1;
        for(int i = c-1; i >= 0; i--){
            LL tx = tp*BSGS(gq, qpow( (i128)h*qpow(g,pi[c]-inv,mod)%mod, pi[i], mod), p, mod);
            inv += tx;
            tp *= p;
        }
        return inv;
    }
    LL exgcd(LL a,LL b,LL &x,LL &y) {
        if(b==0) {
            x=1,y=0;
            return a;
        }
        LL d=exgcd(b,a%b,y,x);
        y-=a/b*x;
        return d;
    }
    LL Inv(LL a,LL P){//求a在膜P下的逆
        if(a==0)return 1;
        LL x,y;
        exgcd(a,P,x,y);
        return (x%P+P)%P;
    }
    //ax+by=c
    LL solve_exgcd(LL a,LL b,LL c,LL &x,LL &y) {
        LL d;
        if(c%(d=__gcd(a,b)))return -1;//返回-1是无解
        a/=d;
        b/=d;
        c/=d;
        exgcd(a,b,x,y);
        x=(i128)x*c%b;
        while(x<=0) x+=b;
        return x;
    }

    LL CRT(vector<LL> xi, vector<node> v){
        LL sz = v.size();
        LL M = 1, ans = 0;
        vector<LL> m;

        for(int i = 0; i<sz; i++){
            m.push_back(qpow(v[i].p,v[i].c));
            M *= m[i];
        }
        for(int i=0; i<sz; i++){
            LL Mi = M/m[i];
            ans = ((i128)ans + (i128)Mi *Inv(Mi,m[i]) *xi[i])%M;
        }
        if(ans < 0) ans += M;
        return ans;
    }

    LL get_x(LL h, LL g, LL N, LL mod, vector<node> &v){
        vector<LL> xi;
        for(node tp : v){
            LL tg = qpow(g,N/qpow(tp.p,tp.c),mod);
            LL th = qpow(h,N/qpow(tp.p,tp.c),mod);
            xi.push_back(get_xi(tg, th, tp.p, tp.c, N, mod));
        }
        return CRT(xi,v);
    }

    //p为质数且(p-1)分解质因子较大的情况时，算法比较优秀，如果ｐ=2*q+1,q为质数, q<=1e12皆可做，则算法与BSGS差不多，不为质数的话可以考虑EXBSGS
    LL solve(LL g, LL h, LL p){
        if(h == 1) return 0;
        LL phi_p = p-1;     
        vector<node> v;
        find_fra(v,phi_p);
        int rt = getroot(p,phi_p,v);
        LL x = get_x(g,rt,phi_p,p,v);
        LL y = get_x(h,rt,phi_p,p,v);
        LL a = 0, b = 0;
        if(x==0){
            if(y==0) return 1;
            else if(y == 1) return 0;
            else return -1; //不满足返回-1
        }
        else return solve_exgcd(x,phi_p,y,a,b);
    }
}

typedef long long LL;
int main(){
#ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
#endif
    int T;
    scanf("%d",&T);
    while(T--){
        LL g,h,p;
        scanf("%lld%lld%lld",&p,&g,&h);
        printf("%lld\n",DLP::solve(g,h,p));
    }
    return 0;
}

poj2417

#include"iostream"
#include"cstdio"
#include"vector"
#include"map"
#include"map"
#include"cstring"
#include"algorithm"
#include"cmath"
using namespace std;

namespace DLP{
	typedef long long LL;
	typedef long long i128;
	const int MX = 1e6+7;
	LL qpow(LL a, LL n){
		LL ret = 1;
		while(n > 0){
			if(n&1) ret = ret*a;
			a = a*a;
			n >>= 1;
		}
		return ret;
	}

	LL qpow(LL a, LL n, LL p){
		LL ret = 1;
		a %= p;
		n %= p;
		while(n > 0){
			if(n&1) ret = (i128)ret*a%p;
			a = (i128)a*a%p;
			n >>= 1;
		}
		return ret;
	}

	int prime[MX],tt;
	bool is_prime[MX];
	struct node{
		LL p;
		int c; 
	};
	void prime_init()
	{
	    memset(is_prime,1,sizeof(is_prime));
	    int n = MX - 7;
	    for(int i = 2; i <= n; i++){
	        if(is_prime[i]) prime[++tt] = i;
	        for(int j = 1; j <= tt && prime[j] <= n/i; j++){
	            is_prime[i*prime[j]] = 0;
	            if(!i%prime[j]) break;
	        }
	    }
	}

	void find_fra(vector<node> &v, LL num){
		if(tt == 0) prime_init();
		for(int i = 1; i <= tt && prime[i] <= num; i++){
			if(num%prime[i] == 0){
				int cnt = 0;
				while(num%prime[i] == 0){
					++cnt;
					num /= prime[i];		
				}
				v.push_back(node{prime[i],cnt});
			}
		}
		if(num > 1) v.push_back(node{num,1});
	}

	int getroot(LL p, LL phi, vector<node> &v) {
	    for(int k=2;; k++) {
	        int flag=1;
	        for(int i=0; i<(int)v.size(); i++)
	            if(qpow(k,phi/v[i].p,p)==1) {
	                flag=0;
	                break;
	            }
	        if(flag)return k;
	    }
	}


	LL BSGS(LL a, LL b, LL p, LL mod) {
		a %= mod, b %= mod;
		//特殊情况，具体看题目要怎么做
		if(b == 1) return 0;
	    if(a==0){
	    	if(b==0) return 1;
	    	else return -1;
	    }

	    LL t = 1;	    
	    int m = int(sqrt(1.0 * p) + 1);
	    LL base = b;
	    map<LL,LL> has; 
	    for(int i = 0; i < m; ++i) {
	        has[base] = i;
	        base = (i128)base * a % mod;
	    }
	    base = qpow(a, m, mod);
		LL now = t;
		//对于A^x，从前往后，保证答案是从小到大的
	    for(int i = 1; i <= m + 1; ++i) {
	        now = (i128)now * base % mod;
	        if(has.count(now)) return i * m - has[now];//找到最小后，找到最小的x使得a^x=1 mod p,通解为 最小+k*x
	    }
	    return -1;
	}

	LL get_xi(LL g, LL h, LL p, LL c, LL N, LL mod){
		vector<LL> pi;
		LL tp = 1;
		for(int i = 0; i <= c; i++){
			pi.push_back(tp);
			tp *= p;
		}
		LL gq = qpow(g,pi[c-1],mod);
		LL inv = 0;
		tp = 1;
		for(int i = c-1; i >= 0; i--){
			LL tx = tp*BSGS(gq, qpow( (i128)h*qpow(g,pi[c]-inv,mod)%mod, pi[i], mod), p, mod);
			inv += tx;
			tp *= p;
		}
		return inv;
	}
	LL exgcd(LL a,LL b,LL &x,LL &y) {
	    if(b==0) {
	        x=1,y=0;
	        return a;
	    }
	    LL d=exgcd(b,a%b,y,x);
	    y-=a/b*x;
	    return d;
	}
	LL Inv(LL a,LL P){//求a在膜P下的逆
	    if(a==0)return 1;
	    LL x,y;
	    exgcd(a,P,x,y);
	    return (x%P+P)%P;
	}
	//ax+by=c
	LL solve_exgcd(LL a,LL b,LL c,LL &x,LL &y) {
	    LL d;
	    if(c%(d=__gcd(a,b)))return -1;//返回-1是无解
	    a/=d;
	    b/=d;
	    c/=d;
	    exgcd(a,b,x,y);
	    x=(i128)x*c%b;
	    while(x<=0) x+=b;
	    return x;
	}

	LL CRT(vector<LL> xi, vector<node> v){
		LL sz = v.size();
		LL M = 1, ans = 0;
		vector<LL> m;

		for(int i = 0; i<sz; i++){
			m.push_back(qpow(v[i].p,v[i].c));
			M *= m[i];
		}
		for(int i=0; i<sz; i++){
			LL Mi = M/m[i];
			ans = ((i128)ans + (i128)Mi *Inv(Mi,m[i]) *xi[i])%M;
		}
		if(ans < 0) ans += M;
		return ans;
	}

	LL get_x(LL h, LL g, LL N, LL mod, vector<node> &v){
		vector<LL> xi;
		for(int i = 0; i < v.size(); i++){
			node tp = v[i];
			LL tg = qpow(g,N/qpow(tp.p,tp.c),mod);
			LL th = qpow(h,N/qpow(tp.p,tp.c),mod);
			xi.push_back(get_xi(tg, th, tp.p, tp.c, N, mod));
		}
		return CRT(xi,v);
	}

	//p为质数且(p-1)分解质因子较大的情况时，算法比较优秀，如果ｐ=2*q+1,q为质数, q<=1e12皆可做，则算法与BSGS差不多，不为质数的话可以考虑EXBSGS
	LL solve(LL g, LL h, LL p){
		if(h == 1) return 0;
		LL phi_p = p-1;		
		vector<node> v;
		find_fra(v,phi_p);
		int rt = getroot(p,phi_p,v);
		LL x = get_x(g,rt,phi_p,p,v);
		LL y = get_x(h,rt,phi_p,p,v);
		LL a = 0, b = 0;
		if(x==0){
	    	if(y==0) return 1;
	    	else if(y == 1) return 0;
	    	else return -1; //不满足返回-1
	    }
	    else return solve_exgcd(x,phi_p,y,a,b);
	}
}

typedef long long LL;
int main(){
#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
#endif
	LL g,h,p;
	while(~scanf("%lld%lld%lld",&p,&g,&h)){
		LL ans = DLP::solve(g,h,p);
		if(ans != -1) printf("%lld\n",ans);
		else printf("no solution\n");
	}
	return 0;
}

 

版权声明：本文为qq_18869763原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。