hdu1845(a^b的因子和%p)
题目链接：http://poj.org/problem?id=1845
思路：
１．整数唯一分解定理：
任意正整数都有且只有一种方式写出其素因子的乘积表达式。
a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn) 其中pi均为素数
２．约数和公式：
对于已经分解的整数a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)
有a的所有因子之和为
S` = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)
那么 a^b 的所有因子和为
S = (1+p1+p1^2+p1^3+...p1^(k1*b)) * (1+p2+p2^2+p2^3+….p2^(k2*b)) * (1+p3+ p3^3+…+ p3^(k3*b)) * .... * (1+pn+pn^2+pn^3+...pn^(kn*b))
对于数列 1, p, p^2, p^3 ... p^n % mod，其中 mod 为质数，打个表可以发现该数列是一个循环数列，其中存在一个循环节为 1, p, p^1, ... p^(mod-2)．其实这点在费马小定理中是有体现的，a^(mod-1) = 1 (% mod)．
那么对于求 cnt = 1 + p + p^2 + ... + p^n % mod，可以令
cc1 = (n + 1) / (mod - 1)
cc2 = (n + 1) % (mod - 1)
cnt1 = 1 + p + p^2 + ... + p^(mod - 2)
cnt2 = 1 + p + p^2 + ... + p^(cc2 - 1)
那么 cnt = cc1 * cnt1 + cnt2 % mod
对于求 a^b，可以先将 a 质因分解，得到 S 的表达式，对于 S 表达式，只需要按照上面的方法求出其中每个乘数项即可．时间复杂度为 O(loga * mod), 本题中 mod = 9901, 时间上是允许的．
代码：
1 #include <iostream> 2 #include <stdio.h> 3 #include <map> 4 #define ll long long 5 using namespace std; 6 7 const int mod = 9901; 8 const int MAXN = 1e5 + 10; 9 ll prime[MAXN], indx = 0; 10 map<int, int> num; 11 12 int get(ll a, ll b){//计算sigma(a^i),其中0<=i<=b 13 ll sol1 = 0, sol2 = -1, cnt = 1; 14 ll cc1 = (b + 1) / (mod - 1); 15 ll cc2 = (b + 1) % (mod - 1); 16 for(int i = 0; i < mod - 1; i++){ 17 sol1 = (sol1 + cnt) % mod; 18 if(sol2 == -1 && i + 1 == cc2) sol2 = sol1; 19 cnt = (cnt * a) % mod; 20 } 21 return (sol1 * cc1 % mod + sol2) % mod; 22 } 23 24 int main(void){ 25 ll a, b, sol = 1; 26 cin >> a >> b; 27 for(int i = 2; i * i <= a; i++){ 28 if(a % i == 0){ 29 prime[indx] = i; 30 while(a % i == 0){ 31 num[i]++; 32 a /= i; 33 } 34 indx++; 35 } 36 } 37 if(a > 1) prime[indx++] = a, num[a]++; 38 for(int i = 0; i < indx; i++){ 39 sol = (sol * get(prime[i], b * num[prime[i]])) % mod; 40 } 41 cout << sol << endl; 42 return 0; 43 }
但是当 mod 比较大时这个方法就不行了，mod 比较大时可以用下面这个代码，贴的 Kuangbin 模板
代码：
1 #include <stdio.h> 2 #include <math.h> 3 #include <iostream> 4 #include <algorithm> 5 #include <string.h> 6 #define ll long long 7 using namespace std; 8 9 //****************************************** 10 //素数筛选和合数分解 11 const int MOD = 9901; 12 const int MAXN=10000; 13 int prime[MAXN+1]; 14 15 void getPrime(void){ 16 memset(prime, 0, sizeof(prime)); 17 for(int i = 2; i <= MAXN; i++) 18 { 19 if(!prime[i]) prime[++prime[0]] = i; 20 for(int j=1; j<= prime[0]&&prime[j] <= MAXN/i; j++) 21 { 22 prime[prime[j]*i] = 1; 23 if(i % prime[j] == 0) break; 24 } 25 } 26 } 27 28 ll factor[100][2]; 29 int fatCnt; 30 31 int getFactors(ll x){ 32 fatCnt = 0; 33 ll tmp = x; 34 for(int i=1; prime[i] <= tmp/prime[i]; i++){ 35 factor[fatCnt][1] = 0; 36 if(tmp % prime[i] == 0){ 37 factor[fatCnt][0] = prime[i]; 38 while(tmp % prime[i] == 0){ 39 factor[fatCnt][1]++; 40 tmp /= prime[i]; 41 } 42 fatCnt++; 43 } 44 } 45 if(tmp!=1) 46 { 47 factor[fatCnt][0]=tmp; 48 factor[fatCnt++][1]=1; 49 } 50 return fatCnt; 51 } 52 53 //****************************************** 54 ll pow_m(ll a, ll n)//快速模幂运算 55 { 56 ll res = 1; 57 ll tmp = a % MOD; 58 while(n){ 59 if(n & 1){ 60 res *= tmp; 61 res%=MOD; 62 } 63 n >>= 1; 64 tmp *= tmp; 65 tmp %= MOD; 66 } 67 return res; 68 } 69 70 ll sum(ll p, ll n){//计算1+p+p^2+...+p^n 71 if(p == 0)return 0; 72 if(n == 0)return 1; 73 if(n & 1){//奇数 74 return ((1 + pow_m(p, n/2 + 1)) % MOD * sum(p, n / 2) % MOD) % MOD; 75 }else return ((1 + pow_m(p, n / 2 + 1)) % MOD * sum(p, n / 2 - 1) + pow_m(p, n / 2) % MOD) % MOD; 76 77 } 78 79 int main(void){ 80 int A, B; 81 getPrime(); 82 while(scanf("%d%d", &A, &B) != EOF){ 83 getFactors(A); 84 ll ans = 1; 85 for(int i = 0; i < fatCnt; i++){ 86 ans *= (sum(factor[i][0], B * factor[i][1]) % MOD); 87 ans %= MOD; 88 } 89 printf("%lld\n",ans); 90 } 91 return 0; 92 }
