Sumdiv
Time Limit: 1000MS   Memory Limit: 30000K
Total Submissions: 18877   Accepted: 4745

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

Source

Romania OI 2002

题意：给定两个正整数和，求的所有因子和对9901取余后的值。

 

分析：很容易知道，先把分解得到，那么得到，那么

     的所有因子和的表达式如下

 

    

 

所以我们有两种做法。第一种做法是二分求等比数列之和。

 

代码：

[cpp]  view plain  copy
 
  在CODE上查看代码片 派生到我的代码片
  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4.   
  5. using namespace std;  
  6. typedef long long LL;  
  7. const int N = 10005;  
  8. const int MOD = 9901;  
  9.   
  10. bool prime[N];  
  11. int p[N];  
  12. int cnt;  
  13.   
  14. void isprime()  
  15. {  
  16.     cnt = 0;  
  17.     memset(prime,true,sizeof(prime));  
  18.     for(int i=2; i<N; i++)  
  19.     {  
  20.         if(prime[i])  
  21.         {  
  22.             p[cnt++] = i;  
  23.             for(int j=i+i; j<N; j+=i)  
  24.                 prime[j] = false;  
  25.         }  
  26.     }  
  27. }  
  28.   
  29. LL power(LL a,LL b)  
  30. {  
  31.     LL ans = 1;  
  32.     a %= MOD;  
  33.     while(b)  
  34.     {  
  35.         if(b & 1)  
  36.         {  
  37.             ans = ans * a % MOD;  
  38.             b--;  
  39.         }  
  40.         b >>= 1;  
  41.         a = a * a % MOD;  
  42.     }  
  43.     return ans;  
  44. }  
  45.   
  46. LL sum(LL a,LL n)  
  47. {  
  48.     if(n == 0) return 1;  
  49.     LL t = sum(a,(n-1)/2);  
  50.     if(n & 1)  
  51.     {  
  52.         LL cur = power(a,(n+1)/2);  
  53.         t = (t + t % MOD * cur % MOD) % MOD;  
  54.     }  
  55.     else  
  56.     {  
  57.         LL cur = power(a,(n+1)/2);  
  58.         t = (t + t % MOD * cur % MOD) % MOD;  
  59.         t = (t + power(a,n)) % MOD;  
  60.     }  
  61.     return t;  
  62. }  
  63.   
  64. void Solve(LL A,LL B)  
  65. {  
  66.     LL ans = 1;  
  67.     for(int i=0; p[i]*p[i] <= A; i++)  
  68.     {  
  69.         if(A % p[i] == 0)  
  70.         {  
  71.             int num = 0;  
  72.             while(A % p[i] == 0)  
  73.             {  
  74.                 num++;  
  75.                 A /= p[i];  
  76.             }  
  77.             ans *= sum(p[i],num*B) % MOD;  
  78.             ans %= MOD;  
  79.         }  
  80.     }  
  81.     if(A > 1)  
  82.     {  
  83.         ans *= sum(A,B) % MOD;  
  84.         ans %= MOD;  
  85.     }  
  86.     cout<<ans<<endl;  
  87. }  
  88.   
  89. int main()  
  90. {  
  91.     LL A,B;  
  92.     isprime();  
  93.     while(cin>>A>>B)  
  94.         Solve(A,B);  
  95.     return 0;  
  96. }  

 

第二种方法就是用等比数列求和公式，但是要用逆元。用如下公式即可

 

                     

 

因为可能会很大，超过int范围，所以在快速幂时要二分乘法。

 

代码：

[cpp]  view plain  copy
 
  在CODE上查看代码片 派生到我的代码片
  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4.   
  5. using namespace std;  
  6. typedef long long LL;  
  7. const int N = 10005;  
  8. const int MOD = 9901;  
  9.   
  10. bool prime[N];  
  11. int p[N];  
  12. int cnt;  
  13.   
  14. void isprime()  
  15. {  
  16.     cnt = 0;  
  17.     memset(prime,true,sizeof(prime));  
  18.     for(int i=2; i<N; i++)  
  19.     {  
  20.         if(prime[i])  
  21.         {  
  22.             p[cnt++] = i;  
  23.             for(int j=i+i; j<N; j+=i)  
  24.                 prime[j] = false;  
  25.         }  
  26.     }  
  27. }  
  28.   
  29. LL multi(LL a,LL b,LL m)  
  30. {  
  31.     LL ans = 0;  
  32.     a %= m;  
  33.     while(b)  
  34.     {  
  35.         if(b & 1)  
  36.         {  
  37.             ans = (ans + a) % m;  
  38.             b--;  
  39.         }  
  40.         b >>= 1;  
  41.         a = (a + a) % m;  
  42.     }  
  43.     return ans;  
  44. }  
  45.   
  46. LL quick_mod(LL a,LL b,LL m)  
  47. {  
  48.     LL ans = 1;  
  49.     a %= m;  
  50.     while(b)  
  51.     {  
  52.         if(b & 1)  
  53.         {  
  54.             ans = multi(ans,a,m);  
  55.             b--;  
  56.         }  
  57.         b >>= 1;  
  58.         a = multi(a,a,m);  
  59.     }  
  60.     return ans;  
  61. }  
  62.   
  63. void Solve(LL A,LL B)  
  64. {  
  65.     LL ans = 1;  
  66.     for(int i=0; p[i]*p[i] <= A; i++)  
  67.     {  
  68.         if(A % p[i] == 0)  
  69.         {  
  70.             int num = 0;  
  71.             while(A % p[i] == 0)  
  72.             {  
  73.                 num++;  
  74.                 A /= p[i];  
  75.             }  
  76.             LL M = (p[i] - 1) * MOD;  
  77.             ans *= (quick_mod(p[i],num*B+1,M) + M - 1) / (p[i] - 1);  
  78.             ans %= MOD;  
  79.         }  
  80.     }  
  81.     if(A > 1)  
  82.     {  
  83.         LL M = MOD * (A - 1);  
  84.         ans *= (quick_mod(A,B+1,M) + M - 1) / (A - 1);  
  85.         ans %= MOD;  
  86.     }  
  87.     cout<<ans<<endl;  
  88. }  
  89.   
  90. int main()  
  91. {  
  92.     LL A,B;  
  93.     isprime();  
  94.     while(cin>>A>>B)  
  95.         Solve(A,B);  
  96.     return 0;  
  97. }