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

原理：

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;