Poj - 1845 - Sumdiv

Sumdiv

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

 

 

 

 

分析：

1、

(1)、如题输入 2 3

        2^3 = 8.    

 8的因数有 1 2 4 8，它们的和为 15

   发现 2^3的因数是：2^0  2^1   2^2  2^3，这些数的和为15

 

（2）、再举例 3 4

  3^4 = 81

81的因数有1 3 9 27 81 ，它们的和为 148

  3^4的因数是： 3^0  3^1   3^2  3^3  3^4 ，这些数的和为 148，与结果符合。

发现共同点 2和3都是素数，若不是素数呢

 

（3）若A非素数如  4 2

      4^2 = 16

16的因数有 1 2 4 8 16 ，它们的和为  31

4^2的因数是：4^0  4^1  4^2 ，它们的和为21，与结果不符。

 

    若把4分解

（2*2）^2 = (2^2)*(2^2)

 即 (2^0 + 2^1 + 2^2)* (2^0 + 2^1 + 2^2) =49 也与结果不符。

 

若分解为（2*2）^2 =2^(2*2)=2^4

即2^0  2^1   2^2  2^3  2^4 ，它们的和为31 ，与结果相符合。

 

所以可总结而出 

     A^B  的所有因数和可转化为

A应分解成（p1^r1）*（p2^r2）*…pk^rk

其中 p1、p2……pk 为A的不重复分解因数

所以 

  故 A^B = p1^(r1*B) * p2^(r2*B) *...* pk^(rk*B);

 

 

举个例子 A=12  B=2时

题面算法 A^B=144

144的因数有 1 2 3 4 6 8 9 12 16 18 24 36 

48 72 144 ，这些数的和为403。

 

推导出各因数的和

A^B=(2*2*3)^2=( (2^2) *3 )^2

=(2^4) * 3^2

Ans=（2^0 + 2^1 + 2^2 + 2^3 + 2^4）*（3^0+3^1+3^2）

 =31*13

 =403      符合

 

 

2、

而现在面临的问题即求从0开始的等比数列和，若用等比求和公式，因为 (0 <= A,B <= 50000000)，所以结果必然会溢出，出现错误，所以求等比数列1+pi+pi^2+pi^3+...+pi^n可以由递归形式的二分求得：（模运算不能用等比数列和公式)

 

 若r为奇数,则数列一共有偶数项，则：
       p^0 + p^1 + p^2 + p^3 +...+ p^r

=     p^0    +    p^1 +    p^2 +……+ p^(r/2)

   + p^(r/2+1) +  p^(r/2+2)+p^(r/2+3)+…+p^r;

（1）前一半与后一半合并相加为

（p^0+ p^(r/2+1)）+（p^1+p^(r/2+2)）+…

…+（p^(r/2)+p^r）；

（2）发现可提取公因式（p^0+ p^(r/2+1)）；

即（p^0+ p^(r/2+1)）*（ p^0+p^1+p^2+……+p^(r/2) ）

如：1+ p + p^2 + p^3 + p^4 + p^5 = (1 + p^3)*(1 + p + p^2) ；

 

 

   若r为偶数,则数列一共有奇数项，会多出一项，则：

    p^0 + p^1 + p^2 + p^3 +...+ p^r

=     p^0   +   p^1  +……+ p^(r/2-1) + p^(r/2)

   + p^(r/2+1) + p^(r/2+2) + … + p^r;

（1）前一半与后一半合并相加为

（p^0+ p^(r/2+1)）+（p^1+p^(r/2+2)）+…

…+（p^(r/2-1)+p^r） +  p^(r/2)；

（2）发现可提取公因式（p^0+ p^(r/2+1)）；

即（p^0+ p^(r/2+1)）*（p^0+p^1+…+p^(r/2-1)）+p^(r/2)；

如：1 + p + p^2 + p^3 + p^4 = (1 + p) * (1 + p^3) + p^2；

 

 

举个例子   p=2,r=5 时

数列为      2^0 + 2^1 + 2^2 + 2^3 + 2^4 + 2^5

   

即 2^0 + 2^1 + 2^2

          + 2^3 + 2^4 + 2^5 

上下合并相加为

  （2^0+2^3）+（2^1+2^4）+ （2^2+2^5）；

提取公因式后为

   （2^0+2^3）* （2^0+2^1+2^2）

此时后面又可看成一个等比数列2^0+2^1+2^2

其为p=2，r=2 的数列

所以用递归的二分求和即可求出等比数列的和。

 

 

 

备注：

1、二分法递归求等比数列和时，多在纸上写几个例子推导即可容易写出。

 

2、用打表出9901内的素数来求A的分解因数仍会超时，所以要用别的方法来求A分解因数

 

3、基于两个最基本公式：

(A*B)%C = ((A%C)*(B%C))%C 

和 (A+B)%C = ((A%C)+(B%C))%C

在二分递归求和的过程中即可对9901取模，既能节省时间，也能避免数字太大。

 

而由于(0 <= A,B <= 50000000)，所以pow(a,b)函数的计算级别太大容易超时，所以要对此写函数优化

需注意如下代码的二分求幂也会超时

 

所以可用位运算来求幂 ：

 

举个例子：

3 ^ 999 = 3 * 3 * 3 * … * 3

直接乘要做998次乘法。

但事实上可以先求出2^k次幂：

3 ^ 2 = 3 * 3
3 ^ 4 = (3 ^ 2) * (3 ^ 2)
3 ^ 8 = (3 ^ 4) * (3 ^ 4)
3 ^ 16 = (3 ^ 8) * (3 ^ 8)
3 ^ 32 = (3 ^ 16) * (3 ^ 16)
3 ^ 64 = (3 ^ 32) * (3 ^ 32)
3 ^ 128 = (3 ^ 64) * (3 ^ 64)
3 ^ 256 = (3 ^ 128) * (3 ^ 128)
3 ^ 512 = (3 ^ 256) * (3 ^ 256)

再相乘：

3 ^ 999
= 3 ^ (512 + 256 + 128 + 64 + 32 + 4 + 2 + 1)
= (3 ^ 512) * (3 ^ 256) * (3 ^ 128) * (3 ^ 64) * (3 ^ 32) * (3 ^ 4) * (3 ^ 2) * 3

这样只要做16次乘法。即使加上一些辅助的存储和运算，也比直接乘高效得多（尤其如果这里底数是成百上千位的大数字的话）。

我们发现，把999转为2进制数：1111100111，其各位就是要乘的数。

 

如输入  2^5时

 5的二进制为 101

 

所以答案为  (1*2^4 ) * (0*2^2) * (1*2^1) =32

#include<iostream>
using namespace std;
const int mod=9901;
int ys[mod];
int cnt[mod];
__int64 pow(__int64 a, __int64 b){     // 位运算。
    __int64 sum = 1;
    while(b > 0){
        if(b & 1)
            sum = (sum * a) % mod;
        b >>= 1;
        a = a * a % mod;
    }
    return sum;
}
__int64 sum(__int64 p,__int64 r)
{
	if(r==0)
	{
		return 1;
	}
	else if(r%2!=0)   //奇数情况
	{
		return sum(p,r/2)%mod*(1+pow(p,r/2+1));	
	}else           //偶数情况
	{
		return sum(p,r/2-1)%mod*(1+pow(p,r/2+1))+pow(p,r/2);	
	}
	
}
int main()
{
	int a,b;
	int i,q;
	while(scanf("%d%d",&a,&b)!=EOF)
	{
		q=0;
		memset(ys,0,sizeof(ys));
		memset(cnt,0,sizeof(cnt));
		for(i=2;i*i<=a;i++)
		{
			if(a%i==0)
			{
				ys[q]=i;
				while(a % i == 0){
					cnt[q] ++;
					a /= i;
				}
				q++;
			}
		}
		if(a != 1){
			ys[q] = a;
			cnt[q] = 1;
			q++;
		}
		
		int res=1;
		for(i=0;i<q;i++)
		{
			res=res*(sum(ys[i],cnt[i] *b)%mod)%mod;
		}
		printf("%d\n",res);
	}
	return 0;
}