【数论】POJ1845 Sumdiv

最近本渣渣做了一道快搞死我的题，就是这个！

下面隆重给出题目以及链接：

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 29696		Accepted: 7312

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

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

题目一眼就能看懂对吧？？不就是求A^B所有因子的和嘛！我第一眼看到题目这么短就心知不好，绝对有坑啊！！搜了一下，竟然要用到数论的知识，学离散的时候数论就学的不好，为了做这道题，把离散初等数论的章节又从头到尾复习了一遍。下面罗列一下这道题要用到的知识点：

埃氏筛法

这是里面最简单的一个知识点了！关于埃氏筛法前几天刚好写了篇随笔，可以戳一下看看。在这里是用来筛素的，提供一个素数表，方便后面的计算。

这是模板：

 1 //埃氏筛法筛素
 2 const int MAX_N = 10005;
 3 int prime[MAX_N];
 4 bool is_prime[MAX_N];
 5 void getPrime(){
 6     int p = 0;
 7     for(int i = 2; i < MAX_N; i++) is_prime[i] = true;
 8     for(int i = 2; i < MAX_N; i++){
 9         if(is_prime[i]){
10             prime[p++] = i;
11             for(int j = 2*i; j < MAX_N ; j+=i) is_prime[j] = false;
12         }
13     }
14 }

算术基本定理

这也是数论里的一个知识点：任何大于1的整数要么是素数，要么可以分解成素数的乘积，这样的分解是唯一的。

也就是说，题目中的A，我们可以把它分解成这样的形式：A = p₁^k₁* p₂^k₂... * p_n^k_n(p_i是素数且k_i>0)

比如36，我们可以把它分解2² * 3²，其中的2和3都是素数；再比如44，我们可以把它分解成2² * 11¹，其中的2和11都是素数。

把这个弄懂之后，显然可以得到：A^B = p₁^k₁B* p₂^k₂B... * p_n^k_nB(p_i是素数且k_i>0）

下面是代码，我们要得到 p_i以及 k_i的数组，以及A的素因子的个数，方便后面的计算：

 1 //素因子分解 
 2 int p[MAX_N],k[MAX_N];
 3 int cnt;
 4 void getFactors(long long x){
 5     memset(k,0,sizeof(k));
 6     cnt = 0;
 7     for(int i = 0; prime[i]*prime[i] <=x ; i++){
 8         if(x % prime[i] == 0){ 
 9             p[cnt] = prime[i];
10             while(x % prime[i] == 0){
11                 k[cnt]++;
12                 x /= prime[i];
13             }
14             cnt++;
15         }
16     }
17     if(x != 1){   
18         p[cnt] = x;
19         k[cnt++] = 1;
20     }
21 }

可能有人（本渣渣）会觉得奇怪，最后那步x！=1是怎么回事呢？注意，在前面for循环的部分，循环条件写的是prime[i]*prime[i]<=x，也就是说，我们只考虑了小于等于√x的素因子。而x最多只可能有一个大于 $\sqrt{x}$ $\sqrt{x}$ $\sqrt{x}$

约数和定理

对于已经分解的整数 A = p₁^k₁* p₂^k₂... * p_n^k_n(p_i是素数且k_i>0)

有A的所有因子之和为 S = (p₁⁰+p₁¹+p₁²+...p₁^k₁) * (p₂⁰+p₂¹+p₂²+...p₂^k₂) ... * (p_n⁰+p_n¹+p_n²+...p_n^k_n)

很容易可以想到，每一项 p_i^k_i都有 k_i+1 个因子，即 p_i⁰、p_i¹、p_i²... p_i^k_i，所以 p_i^k_i的因子和就是 p₁⁰+p₁¹+p₁²+...p₁^k₁。（不难看出，这是一个等比数列。）现在要求A的因子和，因为A是 p_i^k_i的乘积，所以只要把每一项 p_i^k_i的因子和都相乘就可以了。（这是我的理解，想看更严密的证明可以百度一下）

快速幂运算

快速幂也是很基础的一个点了，刚好前段时间也写了随笔，可以戳这里！

代码在这里：

 1 //快速幂运算
 2 long long quickPow(long long a,long long n){
 3     long long res = 1;
 4     while(n){
 5         if(n & 1)
 6             res = res * a % MOD;
 7         a =  a * a % MOD;
 8         n >>= 1;
 9     }
10     return res;
11 }

等比数列二分求和

这个大概是这道题的核心了！从上面给出的约数和定理，我们发现A的所有因子之和是等比数列和的乘积。因为一般的等比数列求和公式中有除法，所以求和时不能直接用公式。据说本题有三种做法，二分，逆元，变换取模。本渣渣现在只会第一种，所以先介绍一下二分求和的方法~

观察等比数列的式子 a + a² + ... + aⁿ ，我们发现：

（1）时，

（2） n 为偶数时，

（3） n 为奇数时，

这边我们用到的式子还要再加上1，公式同样也可以像上面这样推出来：

n为奇数，S_n = （1 + a^n/2+1) * S_n/2；

n为偶数，Sn = （1 + a^n/2+1)S_n/2-1 + a^n/2；

下面给出代码：

 1 //等比数列二分求和
 2 long long sum(long long p,long long n){  //计算1+p+p^2+...+p^n
 3     if(p == 0) return 0;
 4     if(n == 0) return 1;
 5     if(n & 1)
 6         return ((1 + quickPow(p,n/2+1)) * sum(p,n/2)) % MOD;
 7     else
 8         return ((1+quickPow(p,n/2+1)) * sum(p,n/2-1) + quickPow(p,n/2)) % MOD;
 9 
10 }

最后！

下面就可以开始写main函数啦，注意注意，一定要小心取模，本渣渣在这边wa了n次，最后才找到问题。话不多说，上代码：

 1 int main(){
 2     long long A,B;
 3     getPrime();   //得到素数表
 4     while(cin>>A>>B){
 5         getFactors(A);  //分解质因子
 6         long long ans = 1;
 7         for(int i = 0; i < cnt; i++){
 8             ans = ans * sum(p[i],k[i]*B) % MOD; 
 9         }
10         cout<<ans<<endl;
11     }
12 
13     return 0;
14 }

筒子萌！！有没有看到那个for循环！！有没有看到里面的 ans = ans * sum(p[i],k[i]*B) % MOD; 这句！

我之前写的是：ans *= sum(p[i],k[i]*B) % MOD; 然后在for循环外面又写了个 ans %= MOD; 自以为万无一失了，然后。。我真的查了很久的错QAQ。。疯狂wa，内心是崩溃的！

注意在取模运算啥的时候能不写这种 *= 、+= 咱就别写啊，真的挺坑的（大佬请自动忽视我这句话），我前面随处都在留意着取模，没想到这边疏忽了，功亏一篑啊。

下面是完整的代码，有兴趣的可以参考一下，另附我参考的几个题解（题解一题解二）和POJ讨论区的一些测试数据：

 1 #include <iostream>
 2 #include <cstring>
 3 using namespace std;
 4 
 5 #define MOD 9901
 6 
 7 // 求A^B的所有约数之和%9901
 8 
 9 //埃氏筛法筛素
10 const int MAX_N = 10005;
11 int prime[MAX_N];
12 bool is_prime[MAX_N];
13 void getPrime(){
14     int p = 0;
15     for(int i = 2; i < MAX_N; i++) is_prime[i] = true;
16     for(int i = 2; i < MAX_N; i++){
17         if(is_prime[i]){
18             prime[p++] = i;
19             for(int j = 2*i; j < MAX_N ; j+=i) is_prime[j] = false;
20         }
21     }
22 }
23 
24 
25 //素因子分解 A = p1^k1 * p2^k2 * ... * pn^kn
26 int p[MAX_N],k[MAX_N];
27 int cnt;
28 void getFactors(long long x){
29     memset(k,0,sizeof(k));
30     cnt = 0;
31     for(int i = 0; prime[i]*prime[i] <=x ; i++){
32         if(x % prime[i] == 0){  //如果x能被当前素数整除
33             p[cnt] = prime[i];
34             while(x % prime[i] == 0){
35                 k[cnt]++;
36                 x /= prime[i];
37             }
38             cnt++;
39         }
40     }
41     if(x != 1){   //如果x不为1，剩下的一定是大于等于√x的一个素因子
42         p[cnt] = x;
43         k[cnt++] = 1;
44     }
45 }
46 
47 //带模快速幂运算
48 long long quickPow(long long a,long long n){
49     long long res = 1;
50     while(n){
51         if(n & 1)
52             res = res * a % MOD;
53         a =  a * a % MOD;
54         n >>= 1;
55     }
56     return res;
57 }
58 
59 //等比数列二分求和
60 long long sum(long long p,long long n){  //计算1+p+p^2+...+p^n
61     if(p == 0) return 0;
62     if(n == 0) return 1;
63     if(n & 1)
64         return ((1 + quickPow(p,n/2+1)) * sum(p,n/2)) % MOD;
65     else
66         return ((1+quickPow(p,n/2+1)) * sum(p,n/2-1) + quickPow(p,n/2)) % MOD;
67 
68 }
69 
70 int main(){
71     long long A,B;
72     getPrime();   //得到素数表
73     while(cin>>A>>B){
74         getFactors(A);  //分解质因子
75         long long ans = 1;
76         for(int i = 0; i < cnt; i++){
77             ans = ans * sum(p[i],k[i]*B) % MOD;
78         }
79         cout<<ans<<endl;
80     }
81 
82     return 0;
83 }

View Code