poj 1845 Sumdiv （数论）

最新推荐文章于 2019-03-10 11:34:00 发布

clover_hxy

最新推荐文章于 2019-03-10 11:34:00 发布

阅读量246

点赞数

分类专栏：数论

本文链接：https://blog.csdn.net/clover_hxy/article/details/55094581

版权

数论专栏收录该内容

158 篇文章 0 订阅

订阅专栏

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 20394		Accepted: 5144

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

Source

Romania OI 2002

[Submit] [Go Back] [Status] [Discuss]

题解：数论

这道题要求A^B的约数和。

约束和有计算的公式F(n)=∏(1<=i<=k)[pi^(qi+1)-1)/(pi-1)]

证明如下：若a,b互质，则f[ab]=f[a]*f[b] 这个应该比较好理解，因为a,b互质，所以a,b的因数也互质，那么a,b中的所有因数两两自由组合形成的因数一定是ab的因数。

f[a^k]=1+a+a^2+...+a^k ，这个用等比数列求和公式可以直接计算=(a^(k+1)-1)/(a-1)

那么我们将n质因数分解，然后再合并就能得到F(n)=∏(1<=i<=k)[pi^(qi+1)-1)/(pi-1)]

但是这道题中不能保证pi-1与mod互质，即pi-1%mod=0的情况可能存在，那么我们就无法用快速幂或者扩欧来求解了。

有一种比较简单的解决方法就是a/b%c=a*inv(b)%c=(a%(b*c)/b)%c ,但是如果这里的模数变成(pi-1)*mod的话，在进行快速幂运算的时候容易爆LL。

所有这里采用二分递归的方式来求等比数列的和。

1+a+a^2+...+a^k

如果k为奇数，那么数列总共有偶数项（1+a+a^2+...+a^(n/2))*(1+a^*(n/2+1)

如果k为偶数，那么数列总共有奇数项 (1+a+a^2+...+a^(n/2-1)*(1+a^(n/2+1))+a^(n/2)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 500003
#define LL long long 
#define p 9901
using namespace std;
LL a,b,pi[N],qi[N],ans;
LL quickpow(LL num,LL x)
{
	LL base=num%p; 
	LL ans=1;
	while (x) {
		if (x&1) ans=(ans*base)%p;
		x>>=1;
		base=(base*base)%p;
	}
	return ans%p;
}
LL calc(LL a,LL n)
{
	if (n==0) return 1;
	if (n&1) return calc(a,n/2)*(1+quickpow(a,(n/2)+1))%p;
	else return (calc(a,n/2-1)*(1+quickpow(a,(n/2)+1))%p+quickpow(a,n/2))%p;
}
int main()
{
	freopen("a.in","r",stdin);
	freopen("my.out","w",stdout);
	scanf("%I64d%I64d",&a,&b);
	if (a==1) {
		printf("1\n");
		return 0;
	}
	LL x=a; int cnt=0;
	for (int i=2;i*i<=x;i++)
	 if (x%i==0) {
	 	pi[++cnt]=i;
	 	while (x%i==0) {
	 		x/=i;
	 		qi[cnt]++;
		 }
	 }
	if (x>1) pi[++cnt]=x,qi[cnt]=1; 
	ans=1;
	for (int i=1;i<=cnt;i++) {
		LL m=(LL)(pi[i]-1)*(LL)p;
		//cout<<ans<<endl;
		ans=ans*calc(pi[i],qi[i]*b)%p;                               
	}
	printf("%I64d\n",ans);
}