POJ 1845 Sumdiv 求某数的幂取模

Sumdiv

Time Limit: 1000MS Memory Limit: 30000KTotal Submissions: 11049 Accepted: 2617

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).

 

思路：因为任何数都可以表示为 n = p1^e1*p2^e2*p3^e3...pn^en（p1,p2,..,pn为素数),且其因子的和为sum = (1+p1+p1^2+...+p1^e1)*(1+p2+p2^2+...+p2^e2)*...*(1+pn+pn^2+...+pn^en),所以这题的目标就是求sum.假设A = p1^e1*p2^e2*p3^e3...pn^en,则 A^B = (p1^e1*p2^e2*p3^e3...pn^en)^B;故目标实际是求(1+p1+p1^2+...+p1^(e1*B))*....*(1+pn+pn^2+...+pn^(en*B))，令k属于[1,n],很明显目标其实是求pk的前ek*B项和的乘积。又由于A和B都很大，所以要参考以下知识点：

 

1.素数筛选+分解出某整数的素数因子(分解出A的素数因子,也就是求出p1,e1,p2,e2...pn,en）：这个好办。

2.反复平方法求某数的整数幂的取模，也就是a^b%c。(用来求数列前n项和的模用到)

View Code

 1 __int64 Mod(__int64 p,__int64 n,int mod) 2 { 3     __int64 temp = 1; 4     while(n > 0) 5     { 6         if(n % 2) 7             temp = (temp * p) % mod; 8         n /= 2; 9         p = p * p % mod;10     }11     return temp;12 }

 

3.二分法求等比数列的前n项和。

对于等比数列前n项为 1 + q + q^2 + .... + q^n;

（1）若n为奇数，则前n项有偶数项，则Sn = (1+q^(n/2+1)+q*(1+q^(n/2+1)+...+q^(n/2)*(1+q^(n/2+1),提公因式(1+q^(n^2+1)得(1+q+q^2+...+q^(n/2))*(1+q^(n/2+1));

（2）同理，若n为偶数，则前n项有奇数项，则Sn = (1+q+q^2+...+q^(n/2-1))*(1+q^(n/2+1))+q^(n/2);

4.同余模定理  (a+b)%c = (a%c + b%c)%c ;  (a*b)%c = (a%c*b%c)%c;

 

View Code

  1 #include <iostream>  2 #include <cstdio>  3 #include <cstring>  4 #define MAX 10001  5   6 using namespace std;  7   8 struct Node  9 { 10     __int64 p,num; 11 }; 12  13 Node node[MAX]; 14 int prime[MAX]; 15 bool vis[MAX]; 16 int idx; 17 __int64 A,B; 18 int index; 19  20 //筛选素数 21 void findPrime() 22 { 23     memset(vis,0,sizeof(vis)); 24     idx = 0; 25     for(int i=2; i<=100; i++) 26         if(!vis[i]) 27             for(int j=i*i; j<=10000; j+=i) 28                 vis[j] = true; 29     for(int i=2; i<=10000; i++) 30         if(!vis[i]) 31             prime[idx++] = i; 32 } 33 //找出素数因子和它的指数 34 void findSum(int a) 35 { 36     index = 0; 37     for(int i=0; i<idx && prime[i]<=a; i++) 38     { 39         __int64 time = 0; 40         while(a % prime[i] == 0) 41         { 42             a /= prime[i]; 43             time ++; 44         } 45         if(time == 0) 46             continue; 47         node[index].p = prime[i]; 48         node[index].num = time; 49         index ++; 50     } 51     if(a > 1) 52     { 53         node[index].p = a; 54         node[index].num = 1; 55         index ++; 56     } 57 } 58 //反复平方法求幂的模 59 __int64 Mod(__int64 p,__int64 n,int mod) 60 { 61     __int64 temp = 1; 62     while(n > 0) 63     { 64         if(n % 2) 65             temp = (temp * p) % mod; 66         n /= 2; 67         p = p * p % mod; 68     } 69     return temp; 70 } 71  72 //二分法求等比数列前n项和 73 __int64 sum(__int64 p,__int64 n) 74 { 75     if(n == 0) 76         return 1; 77     if(n % 2) 78         return (sum(p,n/2)*(1+Mod(p,n/2+1,9901)))%9901; 79     else 80         return (sum(p,n/2-1)*(1+Mod(p,n/2+1,9901))+Mod(p,n/2,9901))%9901; 81 } 82  83 __int64 solve() 84 { 85     __int64 ret = 1; 86     for(int i=0; i<index; i++) 87     { 88         ret = ((ret%9901)*(sum(node[i].p,node[i].num*B)%9901))%9901; 89     } 90     return ret ; 91 } 92  93 int main() 94 { 95     findPrime(); 96     while(~scanf("%I64d%I64d",&A,&B)) 97     { 98         findSum(A); 99         __int64 ans = solve();100         printf("%I64d\n",ans);101     }102     return 0;103 }