无聊的计算器【数论多合一】

---

题解：

各种算法的证明及结论大家自行问度娘吧

本标程包含所有部分分

代码中有简单解说

ahahaha······

#include<map>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long lo;
inline void Scan(int &x)
{
    char c;
    while((c = getchar()) < '0' || c > '9');
    x = c - '0';
    while((c = getchar()) >= '0' && c <= '9')
        x = x * 10 + c - '0';
}
const int maxn = 1233;
const int maxp = 1000233;
int T;
int type, y, z, p;
inline bool Prime(int p)
{
    int s = sqrt(p);
    for(int i = 2; i <= s; ++i)
        if(!(p % i))
            return false;
    return true;
}
inline int Pow(int x, int n, int p) //快速幂 
{
    int sum = 1;
    while(n)
    {
        if(n & 1) sum = (lo) sum * x % p;
        x = (lo) x * x % p;
        n >>= 1;
    }
    return sum;
}
struct violence //暴力 
{
    int c[maxn][maxn];
    int fac[maxp], inv[maxp];
    inline void Inv() //阶乘逆元 
    {
        fac[0] = 1;
        for(int i = 1; i < p; ++i) fac[i] = (lo) fac[i - 1] * i % p;
        inv[p - 1] = Pow(fac[p - 1], p - 2, p);
        for(int i = p - 2; i >= 0; --i) inv[i] = (lo) inv[i + 1] * (i + 1) % p;
    }
    inline int Bsgs(int a, int b, int p) //求a^x=b(mod p)
    {
        int sum = 1;
        for(int i = 0; i <= p; ++i)
        {
            if(sum == b) return i;
            sum = (lo) sum * a % p;
        }
        return -1;
    }
    inline int Com(int n, int m, int p) //求C(n,m)(mod p) 
    {
        for(int i = 0; i <= n; ++i)
            for(int j = 0; j <= i; ++j)
            {
                if(!j) c[i][j] = 1;
                else c[i][j] = ((lo) c[i - 1][j] + c[i - 1][j - 1]) % p;
            }
        return c[n][m];
    }
    int Lucas(int n, int m, int p) //Lucas求C(n,m)(mod p) 
    {
        if(!m) return 1;
        if(n < m) return 0;
        if(n < p) return (lo) fac[n] * inv[m] % p * inv[n - m] % p;
        return (lo) Lucas(n / p, m / p, p) * Lucas(n % p, m % p, p) % p;
    }
};
violence vio;
struct divisor //Gcd 
{
    inline int Nor(int m, int n) //普通Gcd 
    {
        if(m < n) swap(m, n);
        int r;
        while(r = m % n)
        {
            m = n;
            n = r;
        }
        return n;
    }
    inline int Ex(int a, int b, int &g, lo &x, lo &y) //ExGcd，求ax+by=Gcd(a,b)的解,当为Gcd(a,b)为1时，相当于求a在模p意义下的逆元(答案为x的值) 
    {
        if(!b) g = a, x = 1, y = 0;
        else
        {
            Ex(b, a % b, g, y, x);
            y -= x * (a / b);
        }
    }
}gcd;
inline int Inv(int a, int p) //用ExGcd求逆元 
{
    int g;
    lo x, y;
    gcd.Ex(a, p, g, x, y);
    return ((x % p) + p) % p;
}
struct thoery
{
    int n;
    int fac[maxp]; 
    inline int Fac(int n, int p, int k, lo &e) //计算n!(mod k)，k=p^a，p为质数，e为答案中质因数p的个数 
    {
        if(n < p) return fac[n];
        e += n / p;
        return (lo) Pow(fac[k - 1], n / k, k) * fac[n % k] % k * Fac(n / p, p, k, e) % k;
    }
    /*
        举个例子： n = 19, p = 3, k = 9
        19!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*16*16*17*18*19(mod 9)
           =(1*2*3*4*5*6*7*8*9)*(10*11*12*13*14*15*16*17*18)*19(mod 9)
           =(1*2*4*5*7*8)*(10*11*13*14*16*17)*(3^6)*(1*2*3*4*5*6)*19(mod 9)
           =(1*2*4*5*7*8)*(1*2*4*5*7*8)*(3^6)*(1*2*3*4*5*6)*19(mod 9)
           可以发现左边的两段都是一样的，段数有n/k个
           中间的幂加入e
           右边的还是一个阶乘，递归处理
           最后面的剩下的部分必定小于k个，取用之前求出的阶乘即可 
    */
    inline int Com(int n, int m, int p, int k) //求C(n,m)(mod k), k=p^a, p为质数, C(n,m)=n!/(m!(n-m)!) 
    {
        fac[0] = 1;
        for(int i = 1; i < k; ++i)
        {
            if(i % p) fac[i] = (lo) fac[i - 1] * i % k;
            else fac[i] = fac[i - 1];
        }
        lo t = 0, e = 0;
        lo ansn = Fac(n, p, k, t);
        lo ansm = Fac(m, p, k, e);
        lo ansc = Fac(n - m, p, k, e);
        t -= e;
        return ansn * Inv((lo) ansm * ansc % k, k) % k * Pow(p, t, k) % k;
    }
    inline int Ask(int n, int m, int p, int k, int mo) //求中国剩余定理中模k的答案，a=C(n,m)(mod k)，k=p^a，p为质数 
    {
        int g = mo / k;
        int h = Inv(g, k);
        int a = Com(n, m, p, k);
        return (lo) a * g % mo * h % mo;
    }
    inline int Ans(int n, int m, int p) //分解p的质因数，分成多个质因数的幂相乘，做中国剩余定理 
    {
        int k;
        int c = p;
        int s = sqrt(p);
        int ans = 0;
        for(int i = 2; i <= s; ++i)
        {
            if(!(c % i))
            {
                k = 1;
                while(!(c % i)) k *= i, c /= i;
                ans = (ans + Ask(n, m, i, k, p)) % p;
            }
            if(!(c - 1)) return ans;
        }
        if(c != 1) ans = (ans + Ask(n, m, c, c, p)) % p;
        return ans;
    }
};
thoery crt;
struct algorith
{
    map <int, int> vis;
    inline int Nor(int y, int z, int p) //Bsgs，求y^x=z(mod p)最小非负整数的x，p为质数 
    {
        int s = ceil(sqrt(p));
        vis.clear();
        int x = z;
        for(int i = 0; i <= s; ++i)
        {
            vis[x] = i;
            x = (lo) x * y % p;
        }
        x = 1;
        y = Pow(y, s, p);
        for(int i = 1; i <= s; ++i)
        {
            x = (lo) x * y % p;
            if(vis[x]) return i * s - vis[x];
        }
        return -1;
    }
    /*
        解释一下： 
        a^x=b(mod c)
        设m=ceil(sqrt(c)) 
        设x=i*m-j;
        a^(i*m-j)=b(mod c)
        a^(i*m)=b*(a^J)(mod c)
        a^m^i=b*(a^j)(mod c)
        枚举i将所有a^m^i加入map
        在枚举j查找是否存在b*(a^j) 
    */
    inline int Ex(int a, int b, int c) //普通Bsgs只适用于a与p互质的情况，Exbsgs即将其消至a与p互质，再用普通Bsgs求答案 
    {
        if(!(z - 1)) return 0;
        int g, n = 0, d = 1;
        while(true)
        {
            g = gcd.Nor(a, c);
            if(!(g - 1)) break;
            if(b % g) return -1;
            ++n;
            c /= g, b /= g;
            d = (lo) d * (a / g) % c;
            if(!(b - d)) return n;
        }
        d = (lo) b * Inv(d, c) % c;
        int v = Nor(a, d, c);
        return v < 0 ? -1 : v + n;
    }
    /*
        再解释一下：
        a^x=b(mod c)
        转换为a^x+cy=b
        g=gcd(a,c)
        将整个式子除去g
        如果b不整除g，那么根据扩展欧几里得，方程无解
        不断除去g(对于每个式子求一次),假设除了n次 
        设d为a^n除去g的乘积
        那么d*a^(x-n)=b'(mod c')
        即a^(x-n)=b'/d(mod c') (要是某一时刻b'等于d了，那么b'/d就是1，答案即为n)
        此时c'与a互质了
        我们就可以用普通Bsgs求答案了
        答案加上n就完美啦 
    */
};
algorith bsgs;
int main()
{
    Scan(T);
    while(T--)
    {
        Scan(type), Scan(y), Scan(z), Scan(p);
        bool flag = Prime(p);
        switch(type)
        {
            case 1:
            {
                printf("%d\n", Pow(y, z, p));
                break;
            }
            case 2:
            {
                int ans = (p <= maxn) ? vio.Bsgs(y, z, p) : (flag) ? bsgs.Nor(y, z, p) : bsgs.Ex(y, z, p);
                if(ans != -1) printf("%d\n", ans);
                else printf("Math Error\n");
                break;
            }
            case 3:
            {
                if(z <= maxn && y <= maxn)
                {
                    printf("%d\n", vio.Com(z, y, p));
                    break;
                }
                if(flag)
                {
                    vio.Inv();
                    printf("%d\n", vio.Lucas(z, y, p));
                    break;
                }
                else
                {
                    printf("%d\n", crt.Ans(z, y, p));
                    break;
                }
            }
        }
    }
}