POJ 1845 Sumdiv 求某数的幂取模

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total 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 }