HDU5446
正确的拆乘防爆，超大素数取模，Lucas定理，CRT合并。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 |
#include <iostream> #include <algorithm> #include <queue> #include <vector> #include <cstdlib> #include <cstdio> #include <string> #include <cstring> #include <ctime> #include <iomanip> #include <cmath> #include <set> #include <stack> #include <cmath> #include <map> using namespace std; const int maxn = 100005; long long mod = 1; long long factor[12]; long long A[12]; long long fac[12][maxn]; inline long long mul(long long a,long long b,long long moder){ long long ret = 0; if (a > b) { swap(a, b); } long long fac = a; while (b) { if (b & 1) { ret += fac; if (ret >= moder) { ret -= moder; } } fac += fac; if (fac >= moder) { fac -= moder; } b>>= 1; } return ret; } long long pow_mod(long long n, long long m, long long p) { n %= p; long long ret = 1; for(; m; m>>=1, n = mul(n, n, p)) if(m & 1) ret = mul(ret, n, p); return ret; } long long inv(long long n, long long p) { return pow_mod(n, p-2, p); } long long C(long long n, long long m, int i,long long p) { if(n < m) return 0; return mul(fac[i][n], inv(mul(fac[i][m], fac[i][n-m], p), factor[i]),p); } long long lucas(long long n, long long m, int i) { if(m == 0) return 1; long long p = factor[i]; return mul(lucas(n/p, m/p, i), C(n%p, m%p, i,p), p); } long long CRT(int k) { long long m = mod, x = 0; for(int i = 0; i < k; i++) { long long w = m / factor[i]; x += mul(mul(inv(w, factor[i]), A[i], m), w, m); x %= m; } return (x + m) % m; } int main() { int caseCnt = 0; scanf(" %d",&caseCnt); for (int d = 1; d <= caseCnt; d++) { long long n,m,k; memset(fac, 0, sizeof(fac)); mod = 1; scanf(" %lld %lld %lld",&n,&m,&k); long long maxfac = 0; for (int i = 0; i < k; i++) { cin >> factor[i]; maxfac = max(factor[i],maxfac); mod *= factor[i]; } for(int i = 0; i < k; i++) { fac[i][0] = 1; for(int j = 1; j <= maxfac; j++) fac[i][j] = mul(fac[i][j-1], j, factor[i]); } memset(A, 0, sizeof(A)); for(int i = 0; i < k; i++) { A[i] = (A[i] + lucas(n, m, i)) % factor[i]; } printf("%lld\n",CRT(k)); } return 0; } |
(Extended) Baby Step Giant Step Algorithm
用于求出满足A^x=B(mod C)的最小的x。
- 参数配置：
- HASHMOD：Hashtable的Hash空间大小，爆内存要改小
- 调用方法：exBSGS(A,B,C);
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 |
#pragma mark - Extended Baby-Step-Gaint-Step Algorithm (exBSGS) struct Hashtable{//Mod Based Hash Table static const int HASHMOD = 1200611;//A Prime long long top, hash[HASHMOD+100], value[HASHMOD+100], stk[1 << 16]; int locate(const long long x){ int h = x % HASHMOD; while(hash[h] != -1 && hash[h] != x){ h++; } return h; } void insert(const long long key,const long long val){ int h = locate(key); if(hash[h] == -1){ hash[h] = key; value[h] = val; stk[++top] = h; } } long long find(const long long key){ int h = locate(key); return (hash[h] == key)?value[h]:-1; } void clear(){ while(top){ hash[stk[top--]] = -1; } } Hashtable(){ top = 0; memset(hash, -1, sizeof(hash)); } } hashtab; struct Triple{ long long x,y,z; Triple(const long long _x,const long long _y,const long long _z){ x = _x; y = _y; z = _z; } Triple(){ x = y = z = 0; } }; Triple exgcd(const long long a,const long long b){ //for inv use(A/B % C):let a = B,b = C,return.x = inv,return.z = gcd //for solve use:return.x = specx,return.y = specy,returnz = gcd if (b == 0) { return Triple(1,0,a); } Triple last = exgcd(b, a%b); return Triple(last.y, last.x - a / b * last.y, last.z); } long long exBSGS(long long A,long long B,long long C){ //A^x=B(mod C) //Ensure B is legal B %= C; if (B < 0) { B += C; } //Procceed A,B,C to normal BSGS if C is not prime long long tmp = 1 % C, cnt = 0, D = 1 % C; for (int i = 0; i < 64; i++) { if (tmp == B) { return i; } tmp = tmp * A % C; } for (Triple res; (res = exgcd(A, C)).z != 1; cnt++) { if (B % res.z) { return -1; } B /= res.z; C /= res.z; D = D * A / res.z % C; } //Normal BSGS long long sqrtn = (long long)(ceil(sqrt(C))); hashtab.clear(); long long base = 1 % C; for (int i = 0; i < sqrtn; i++) { hashtab.insert(base, i); base = base * A % C; } long long j = -1,i = 0; for (; i < sqrtn; i++) { Triple res = exgcd(D, C); long long c = C / res.z; res.x = (res.x * B/res.z % c + c) % c; j = hashtab.find(res.x); if (j != -1) { return i * sqrtn + j + cnt; } D = D * base % C; } return -1; } #pragma mark - |
指数循环节
若Ａ^start = B (mod P)求出最小的len使得Ａ^(k*len) = B (mod P)，原理是这种对于一个值的循环节长度一定是Phi(P)的约数，所以不断约掉因子并尝试。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
long long fastpow(long long a,long long p,long long m){ long long ret = 1; long long fac = a; while(p){ if(p&1){ ret *= fac; ret %= m; } p >>= 1; fac *= fac; fac %= m; } return ret; } long long euler_phi(long long x){ long long ret = 1; for(long long i = 2; i * i <= x; i++){ if(x % i == 0){ x /= i; ret *= i-1; while(x % i == 0){ x /= i; ret *= i; } } } if(x > 1){ ret *= x-1; } return ret; } long long loop_len(long long start,long long A,long long D,long long P){ long long phi = euler_phi(P),ret = phi; for(long long i = 2; i * i <= phi; i++){ while(phi % i == 0){ phi /= i; } while(ret % i == 0 && fastpow(A,start + ret/i,P) == D){ ret /= i; } } if (phi > 1) { while(ret % phi == 0 && fastpow(A,start + ret/phi,P) == D){ ret /= phi; } } return ret; } |
