POJ - 1845 Sumdiv —— 因子和 等比数列求和

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 25199		Accepted: 6244

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

题意：给定a和b，求a的b次方的所有因子的和

思路：

由算术基本定理将a拆分成素数的幂积，即a=p1^k1+p2^k2+...pn^kn

则a^b=p1^k1b+p2^k2b+...pn^knb

那么其因子和 (1+p1+p1^2+..p1^k1b)*(1+p2+p2^2+...p2^k2b)*...(1+pn+pn^2+...pn^knb)

 ___ps:不理解的话可以将括号乘进去，运算一下就明白了

要求这个式子,可以把每个括号看成等比数列的前n项和

但是若用前n项和公式求，比较麻烦，会涉及到乘法逆元，我鼓捣了半天没弄出来....

推荐的办法是用递归求其前n项和，理论如下：

①若n为奇数,一共有偶数项，则：

1 + p + p^2 + p^3 +...+ p^n 

= (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

前半部分恰好就是原式的一半，那么只需要不断递归二分求和就可以了，后半部分为幂次式，用快速幂解。
②若n为偶数,一共有奇数项,则:

1 + p + p^2 + p^3 +...+ p^n

= (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

前半部分还是原式的一半

这样求出每个的和加上快速幂就能搞出来了

需要注意的几个坑点：

一是特判a==1的情况

二是分解a的时候，要遍历到根号a下，否则会T，最后再判断一下有没有取因子取尽，没有取尽则最后肯定剩的是个素数

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <cmath>
#include <vector>
#define max_ 50000100
#define inf 0x3f3f3f3f
#define ll long long
#define les 1e-8
#define mod 9901
using namespace std;
ll a,b;
int pl=0;
ll num[100100];
ll cnt[100100];
ll fpow(ll a,ll b)
{
	ll tmp=a%mod;
	ll ans=1;
	while(b)
	{
		if(b&1)
		ans=(ans*tmp)%mod;
		tmp=(tmp*tmp)%mod;
		b/=2;
	}
	return (ans%mod);
}
ll f(ll p,ll n)//求等比数列的前n项和
{
	if(n==0)
	return 1;
	if(n&1)
	return (f(p,n/2)*(1+fpow(p,n/2+1)))%mod;
	else
	return (f(p,n/2-1)*(1+fpow(p,n/2+1))%mod+fpow(p,n/2)%mod)%mod;
}
int main(int argc, char const *argv[]) {
	scanf("%lld%lld",&a,&b);
	if(a==1LL)
	{
		printf("1\n" );
		return 0;
	}
	ll tmp=a;
	int flag=0;
	for(ll i=2;i*i<=a;i++)
	{
		if(tmp%i==0)
		{
			num[++pl]=i;
			flag=1;
		}
		while(tmp%i==0)
		{
			tmp/=i;
			cnt[pl]++;
		}
		if(tmp==1)
		break;
	}
	if(tmp!=1)
	{
		num[++pl]=tmp;
		cnt[pl]=1;
	}
	ll ans=1;
	for(ll i=1;i<=pl;i++)
	{
		ans=(ans*(f(num[i],cnt[i]*b)%mod))%mod;
	}
	printf("%lld\n",ans%mod );
	return 0;
}