Sumdiv[POJ1845]

欢迎大家访问我的老师的OJ———caioj.cn

题面描述

传送门

思路

一句话题意:求$A^B$的所有约数之和$\mathrm{mod}9901(1\le A,B\le 5\ast 10^7)$

提要:

等比数列求和公式:$S_n=a_1\ast \frac{1-q^n}{1-q}$，其中$q$为公比(与公差类似)

乘法逆元:

若整数$b,m$互质，且$b\mid a$，则存在一个整数x，使得$a/b\equiv a\ast x(\mathrm{mod}m)$，
则称$x$为$b$的模$m$乘法逆元，记为$b^{-1}(\mathrm{mod}m)$。

因为$a/b\equiv a\ast b^{-1}\equiv a/b\ast b\ast b^{-1}(\mathrm{mod}m)$，所以$b\ast b^{-1}\equiv 1(\mathrm{mod}m)$

如果$m$是质数的话,并且$b<m$,根据费马小定理，可知，$b^{m-1}\equiv 1(\mathrm{mod}m)，即b\ast b^{m-2}\equiv 1(\mathrm{mod}m)$。

因此，当模数$m$为质数时，$b^{m-2}$即为$b$的乘法逆元。

由于我们这里仅涉及到$m$为$9901$，即为质数的情况，所以不做拓展。

由于A分解质因数，可表示为:$p_1^{c_1}\ast p_2^{c_2}\ast ......\ast p_m^{c_m}$，根据算术基本定理的“约数和”推论，

$A^B$的所有约数之和为:$(1+p_1+p_1^2+......p_1^{B\ast c_1})\ast ......\ast (1+p_m+p_m^2+......+p_m^{B\ast c_m})=\prod _{i=1}^m(\sum _{j=0}^{B\ast c_i}(p_i{)}^j)$，

仔细观察可以发现，$\sum _{j=0}^{B\ast c_i}(p_i{)}^j$其实就是一个等比数列求和,

即$(p_i^{B\ast c_i+1}-1)/(p_i-1)$，这是我们只要对分子$p_i^{B\ast c_i+1}-1$快速幂求出，
分母通过快速幂，求出乘法逆元得出即可。

特别地，对于$p_i-1$是$9901$的倍数,此时乘法逆元并不存在，但是$p_1\mathrm{mod}9901=1$，所以$1+p_1+p_1^2+......+p_1^{b\ast c_1}\equiv 1+1+1^2+......+1^{B\ast c_1}\equiv B\ast c_1+1(\mathrm{mod}9901).$

AC code

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<cstdlib>
#define ll long long
#define gc getchar()
using namespace std;
const int N=15;
const int mod=9901;
ll p[N];int c[N];
inline void qr(ll &x)
{
	x=0;int f=1;char c=gc;
	while(c<'0'||c>'9'){if(c=='-')f=-1;c=gc;}
	while(c>='0'&&c<='9'){x=x*10+(c^48);c=gc;}
	x*=f;
}
inline void mul(ll &a,ll b)
{
	ll c=a*b/mod;
	a=a*b-c*mod;
}
inline ll pow_mod(ll a,ll b)
{
	ll ans=1;a%=mod;
	while(b)
	{
		if(b&1)mul(ans,a);
		b>>=1;mul(a,a);
	}
	return ans;
}
inline void div(ll x)
{
	p[0]=0;
	for(int i=2;i*i<=x;i++)
	{
		if(x%i==0)
		{
			p[++p[0]]=i;
			c[p[0]]=0;
			while(x%i==0)x/=i,++c[p[0]];
		}
	}
	if(x>1)p[++p[0]]=x,c[p[0]]=1;
}
int main()
{
	ll a,b;qr(a),qr(b);
	div(a);ll ans=1;
	for(int i=1;i<=p[0];i++)
	{
		if((p[i]-1)%mod==0)
		{
			ans=(b*c[i]+1)*ans%mod;
			continue;
		}
		ll x=pow_mod(p[i],b*c[i]+1);
		x=(x-1+mod)%mod;
		ll y=p[i]-1;
		y=pow_mod(y,mod-2);
		ans=ans*x*y%mod; 
	}
	printf("%lld\n",ans);
	return 0;
}