POJ 1845 Sumdiv 【二分 || 逆元】
任意门:http://poj.org/problem?id=1845、
Sumdiv
Time Limit: 1000MS |
|
Memory Limit: 30000K |
Total Submissions: 30268 |
|
Accepted: 7447 |
Description
Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).
Input
The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.
Output
The only line of the output will contain S modulo 9901.
Sample Input
2 3
Sample Output
15
Hint
2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).
Source
题目概括:
给出 A 和 B ,求 AB 所有因子之和,答案 mod 9901;
解题思路:
可对 A 先进行素因子分解 A = p1a1 * p2a2 * p3a3 * ...... * pnan
即 AB = p1a1*B * p2a2*B * p3a3*B * ...... * pnan*B
那么 AB 所有因子之和 = (1 + P11 + P12 + P13 + ... + P1a1*B) * (1 + P21 + P22 + P23 + ... + P2a2*B) * ......* (1 + Pn1 + Pn2 + Pn3 + ... + Pnan*B) ;
对于 (1 + P1 + P2 + P3 + ... + Pa1*B) 我们可以
①二分求和 (47 ms)
AC code:
1 #include <cstdio> 2 #include <iostream> 3 #include <cstring> 4 #include <algorithm> 5 #include <cmath> 6 #define INF 0x3f3f3f3f 7 #define LL long long 8 using namespace std; 9 const LL MOD = 9901; 10 const int MAXN = 1e4+10; 11 LL p[MAXN]; 12 bool isprime[MAXN]; 13 int cnt; 14 LL A, B; 15 16 void init() //打表预处理素因子 17 { 18 cnt = 1; 19 memset(isprime, 1, sizeof(isprime)); 20 for(LL i = 2; i < MAXN; i++){ 21 if(isprime[i]){ 22 p[cnt++] = i; 23 for(int j = 1; j < cnt && p[j]*i < MAXN; j++) 24 isprime[p[j]*i] = false; 25 } 26 } 27 } 28 29 LL qpow(LL x, LL n) 30 { 31 LL res = 1LL; 32 while(n){ 33 if(n&1) res = (res%MOD*x%MOD)%MOD; 34 x = (x%MOD*x%MOD)%MOD; 35 n>>=1LL; 36 } 37 return res; 38 } 39 40 LL sum(LL x, LL n) 41 { 42 if(n == 0) return 1; 43 LL res = sum(x, (n-1)/2); 44 if(n&1){ 45 res = (res + res%MOD*qpow(x, (n+1)/2)%MOD)%MOD; 46 } 47 else{ 48 res = (res + res%MOD*qpow(x, (n+1)/2)%MOD)%MOD; 49 res = (res + qpow(x, n))%MOD; 50 } 51 return res; 52 } 53 54 int main() 55 { 56 LL ans = 1; 57 scanf("%lld %lld", &A, &B); 58 init(); 59 //printf("%d\n", cnt); 60 for(LL i = 1; p[i]*p[i] <= A && i < cnt; i++){ //素因子分解 61 //printf("%lld\n ", p[i]); 62 if(A%p[i] == 0){ 63 LL num = 0; 64 while(A%p[i] == 0){ 65 num++; 66 A/=p[i]; 67 } 68 ans = ((ans%MOD*sum(p[i], num*B)%MOD)+MOD)%MOD; 69 } 70 } 71 if(A > 1){ 72 // puts("zjy"); 73 ans = (ans%MOD*sum(A, B)%MOD+MOD)%MOD; 74 } 75 76 printf("%lld\n", ans); 77 78 return 0; 79 }
②应用等比数列求和公式 ( 0ms )
因为要保证求模时的精度,所以要求逆元。
这里用的方法是一般情况都适用的 : A/B mod p = A mod (p*B) / B;
但是考虑乘法会爆int,所以自定义一个二分乘法。
AC code:
1 // 逆元 + 二分乘法 2 #include <cstdio> 3 #include <iostream> 4 #include <cmath> 5 #include <algorithm> 6 #include <cstring> 7 #define INF 0x3f3f3f3f 8 #define LL long long 9 using namespace std; 10 const int MAXN = 1e4+10; 11 const LL mod = 9901; 12 int p[MAXN], cnt; 13 bool isp[MAXN]; 14 LL A, B; 15 16 void prime() 17 { 18 memset(isp, 1, sizeof(isp)); 19 isp[0] = isp[1] = false; 20 for(int i = 2; i < MAXN; i++){ 21 if(isp[i]){ 22 p[cnt++] = i; 23 for(int k = 0; k < cnt && i*p[k] < MAXN; k++){ 24 isp[i*p[k]] = false; 25 } 26 } 27 } 28 } 29 30 LL mutil(LL x, LL y, LL MOD) 31 { 32 LL res = 0; 33 while(y){ 34 if(y&1) res = (res+x)%MOD; 35 x = (x+x)%MOD; 36 y>>=1LL; 37 } 38 return res; 39 } 40 41 LL q_pow(LL x, LL n, LL MOD) 42 { 43 LL res = 1LL; 44 while(n){ 45 if(n&1) res = mutil(res, x, MOD)%MOD; 46 x = mutil(x, x, MOD)%MOD; 47 n>>=1LL; 48 } 49 return res; 50 } 51 52 int main() 53 { 54 LL num = 0; 55 LL ans = 1; 56 scanf("%lld %lld", &A, &B); 57 prime(); 58 for(int i = 0; p[i]*p[i] <= A; i++){ 59 if(A%p[i] == 0){ 60 num = 0; 61 while(A%p[i] == 0){ 62 A/=p[i]; 63 num++; 64 } 65 LL m = mod*(p[i]-1); 66 ans *= (q_pow(p[i], num*B+1, m) + m-1)/(p[i] - 1); 67 ans = ans%mod; 68 //ans += ((q_pow(p[i], num*B, mod*(p[i]-1))-p[i])/(p[i]-1))%mod; 69 } 70 } 71 72 if(A != 1){ 73 // ans += ((q_pow(A, B, mod*(A-1))-A)/(A-1))%mod; 74 LL m = mod*(A-1); 75 ans*=(q_pow(A, B+1, m)+m-1)/(A-1); 76 ans = ans%mod; 77 } 78 79 printf("%lld\n", ans); 80 81 return 0; 82 }