题解：
各种算法的证明及结论大家自行问度娘吧
本标程包含所有部分分
代码中有简单解说
ahahaha······
1 #include<map> 2 #include<cmath> 3 #include<cstdio> 4 #include<cstdlib> 5 #include<cstring> 6 #include<iostream> 7 #include<algorithm> 8 using namespace std; 9 typedef long long lo; 10 inline void Scan(int &x) 11 { 12 char c; 13 while((c = getchar()) < '0' || c > '9'); 14 x = c - '0'; 15 while((c = getchar()) >= '0' && c <= '9') 16 x = x * 10 + c - '0'; 17 } 18 const int maxn = 1233; 19 const int maxp = 1000233; 20 int T; 21 int type, y, z, p; 22 inline bool Prime(int p) 23 { 24 int s = sqrt(p); 25 for(int i = 2; i <= s; ++i) 26 if(!(p % i)) 27 return false; 28 return true; 29 } 30 inline int Pow(int x, int n, int p) //快速幂 31 { 32 int sum = 1; 33 while(n) 34 { 35 if(n & 1) sum = (lo) sum * x % p; 36 x = (lo) x * x % p; 37 n >>= 1; 38 } 39 return sum; 40 } 41 struct violence //暴力 42 { 43 int c[maxn][maxn]; 44 int fac[maxp], inv[maxp]; 45 inline void Inv() //阶乘逆元 46 { 47 fac[0] = 1; 48 for(int i = 1; i < p; ++i) fac[i] = (lo) fac[i - 1] * i % p; 49 inv[p - 1] = Pow(fac[p - 1], p - 2, p); 50 for(int i = p - 2; i >= 0; --i) inv[i] = (lo) inv[i + 1] * (i + 1) % p; 51 } 52 inline int Bsgs(int a, int b, int p) //求a^x=b(mod p) 53 { 54 int sum = 1; 55 for(int i = 0; i <= p; ++i) 56 { 57 if(sum == b) return i; 58 sum = (lo) sum * a % p; 59 } 60 return -1; 61 } 62 inline int Com(int n, int m, int p) //求C(n,m)(mod p) 63 { 64 for(int i = 0; i <= n; ++i) 65 for(int j = 0; j <= i; ++j) 66 { 67 if(!j) c[i][j] = 1; 68 else c[i][j] = ((lo) c[i - 1][j] + c[i - 1][j - 1]) % p; 69 } 70 return c[n][m]; 71 } 72 int Lucas(int n, int m, int p) //Lucas求C(n,m)(mod p) 73 { 74 if(!m) return 1; 75 if(n < m) return 0; 76 if(n < p) return (lo) fac[n] * inv[m] % p * inv[n - m] % p; 77 return (lo) Lucas(n / p, m / p, p) * Lucas(n % p, m % p, p) % p; 78 } 79 }; 80 violence vio; 81 struct divisor //Gcd 82 { 83 inline int Nor(int m, int n) //普通Gcd 84 { 85 if(m < n) swap(m, n); 86 int r; 87 while(r = m % n) 88 { 89 m = n; 90 n = r; 91 } 92 return n; 93 } 94 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的值) 95 { 96 if(!b) g = a, x = 1, y = 0; 97 else 98 { 99 Ex(b, a % b, g, y, x); 100 y -= x * (a / b); 101 } 102 } 103 }gcd; 104 inline int Inv(int a, int p) //用ExGcd求逆元 105 { 106 int g; 107 lo x, y; 108 gcd.Ex(a, p, g, x, y); 109 return ((x % p) + p) % p; 110 } 111 struct thoery 112 { 113 int n; 114 int fac[maxp]; 115 inline int Fac(int n, int p, int k, lo &e) //计算n!(mod k)，k=p^a，p为质数，e为答案中质因数p的个数 116 { 117 if(n < p) return fac[n]; 118 e += n / p; 119 return (lo) Pow(fac[k - 1], n / k, k) * fac[n % k] % k * Fac(n / p, p, k, e) % k; 120 } 121 /* 122 举个例子： n = 19, p = 3, k = 9 123 19!=1*2*3*4*5*6*7*8*9*10*11*12*13*14*16*16*17*18*19(mod 9) 124 =(1*2*3*4*5*6*7*8*9)*(10*11*12*13*14*15*16*17*18)*19(mod 9) 125 =(1*2*4*5*7*8)*(10*11*13*14*16*17)*(3^6)*(1*2*3*4*5*6)*19(mod 9) 126 =(1*2*4*5*7*8)*(1*2*4*5*7*8)*(3^6)*(1*2*3*4*5*6)*19(mod 9) 127 可以发现左边的两段都是一样的，段数有n/k个 128 中间的幂加入e 129 右边的还是一个阶乘，递归处理 130 最后面的剩下的部分必定小于k个，取用之前求出的阶乘即可 131 */ 132 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)!) 133 { 134 fac[0] = 1; 135 for(int i = 1; i < k; ++i) 136 { 137 if(i % p) fac[i] = (lo) fac[i - 1] * i % k; 138 else fac[i] = fac[i - 1]; 139 } 140 lo t = 0, e = 0; 141 lo ansn = Fac(n, p, k, t); 142 lo ansm = Fac(m, p, k, e); 143 lo ansc = Fac(n - m, p, k, e); 144 t -= e; 145 return ansn * Inv((lo) ansm * ansc % k, k) % k * Pow(p, t, k) % k; 146 } 147 inline int Ask(int n, int m, int p, int k, int mo) //求中国剩余定理中模k的答案，a=C(n,m)(mod k)，k=p^a，p为质数 148 { 149 int g = mo / k; 150 int h = Inv(g, k); 151 int a = Com(n, m, p, k); 152 return (lo) a * g % mo * h % mo; 153 } 154 inline int Ans(int n, int m, int p) //分解p的质因数，分成多个质因数的幂相乘，做中国剩余定理 155 { 156 int k; 157 int c = p; 158 int s = sqrt(p); 159 int ans = 0; 160 for(int i = 2; i <= s; ++i) 161 { 162 if(!(c % i)) 163 { 164 k = 1; 165 while(!(c % i)) k *= i, c /= i; 166 ans = (ans + Ask(n, m, i, k, p)) % p; 167 } 168 if(!(c - 1)) return ans; 169 } 170 if(c != 1) ans = (ans + Ask(n, m, c, c, p)) % p; 171 return ans; 172 } 173 }; 174 thoery crt; 175 struct algorith 176 { 177 map <int, int> vis; 178 inline int Nor(int y, int z, int p) //Bsgs，求y^x=z(mod p)最小非负整数的x，p为质数 179 { 180 int s = ceil(sqrt(p)); 181 vis.clear(); 182 int x = z; 183 for(int i = 0; i <= s; ++i) 184 { 185 vis[x] = i; 186 x = (lo) x * y % p; 187 } 188 x = 1; 189 y = Pow(y, s, p); 190 for(int i = 1; i <= s; ++i) 191 { 192 x = (lo) x * y % p; 193 if(vis[x]) return i * s - vis[x]; 194 } 195 return -1; 196 } 197 /* 198 解释一下： 199 a^x=b(mod c) 200 设m=ceil(sqrt(c)) 201 设x=i*m-j; 202 a^(i*m-j)=b(mod c) 203 a^(i*m)=b*(a^J)(mod c) 204 a^m^i=b*(a^j)(mod c) 205 枚举i将所有a^m^i加入map 206 在枚举j查找是否存在b*(a^j) 207 */ 208 inline int Ex(int a, int b, int c) //普通Bsgs只适用于a与p互质的情况，Exbsgs即将其消至a与p互质，再用普通Bsgs求答案 209 { 210 if(!(z - 1)) return 0; 211 int g, n = 0, d = 1; 212 while(true) 213 { 214 g = gcd.Nor(a, c); 215 if(!(g - 1)) break; 216 if(b % g) return -1; 217 ++n; 218 c /= g, b /= g; 219 d = (lo) d * (a / g) % c; 220 if(!(b - d)) return n; 221 } 222 d = (lo) b * Inv(d, c) % c; 223 int v = Nor(a, d, c); 224 return v < 0 ? -1 : v + n; 225 } 226 /* 227 再解释一下： 228 a^x=b(mod c) 229 转换为a^x+cy=b 230 g=gcd(a,c) 231 将整个式子除去g 232 如果b不整除g，那么根据扩展欧几里得，方程无解 233 不断除去g(对于每个式子求一次),假设除了n次 234 设d为a^n除去g的乘积 235 那么d*a^(x-n)=b'(mod c') 236 即a^(x-n)=b'/d(mod c') (要是某一时刻b'等于d了，那么b'/d就是1，答案即为n) 237 此时c'与a互质了 238 我们就可以用普通Bsgs求答案了 239 答案加上n就完美啦 240 */ 241 }; 242 algorith bsgs; 243 int main() 244 { 245 Scan(T); 246 while(T--) 247 { 248 Scan(type), Scan(y), Scan(z), Scan(p); 249 bool flag = Prime(p); 250 switch(type) 251 { 252 case 1: 253 { 254 printf("%d ", Pow(y, z, p)); 255 break; 256 } 257 case 2: 258 { 259 int ans = (p <= maxn) ? vio.Bsgs(y, z, p) : (flag) ? bsgs.Nor(y, z, p) : bsgs.Ex(y, z, p); 260 if(ans != -1) printf("%d ", ans); 261 else printf("Math Error "); 262 break; 263 } 264 case 3: 265 { 266 if(z <= maxn && y <= maxn) 267 { 268 printf("%d ", vio.Com(z, y, p)); 269 break; 270 } 271 if(flag) 272 { 273 vio.Inv(); 274 printf("%d ", vio.Lucas(z, y, p)); 275 break; 276 } 277 else 278 { 279 printf("%d ", crt.Ans(z, y, p)); 280 break; 281 } 282 } 283 } 284 } 285 }