poj 1845 Sumdiv(求逆元)

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 17590		Accepted: 4421

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

Source

Romania OI 2002

题意：给定两个正整数

和

，求

的所有因子和对9901取余后的值。

分析：很容易知道，先把

分解得到

，那么得到

，那么


的所有因子和的表达式如下



所以我们有两种做法。第一种做法是二分求等比数列之和。
代码：

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 10005;
const int MOD = 9901;

bool prime
;
int p
;
int cnt;

void isprime()
{
cnt = 0;
memset(prime,true,sizeof(prime));
for(int i=2; i<N; i++)
{
if(prime[i])
{
p[cnt++] = i;
for(int j=i+i; j<N; j+=i)
prime[j] = false;
}
}
}

LL power(LL a,LL b)
{
LL ans = 1;
a %= MOD;
while(b)
{
if(b & 1)
{
ans = ans * a % MOD;
b--;
}
b >>= 1;
a = a * a % MOD;
}
return ans;
}

LL sum(LL a,LL n)
{
if(n == 0) return 1;
LL t = sum(a,(n-1)/2);
if(n & 1)
{
LL cur = power(a,(n+1)/2);
t = (t + t % MOD * cur % MOD) % MOD;
}
else
{
LL cur = power(a,(n+1)/2);
t = (t + t % MOD * cur % MOD) % MOD;
t = (t + power(a,n)) % MOD;
}
return t;
}

void Solve(LL A,LL B)
{
LL ans = 1;
for(int i=0; p[i]*p[i] <= A; i++)
{
if(A % p[i] == 0)
{
int num = 0;
while(A % p[i] == 0)
{
num++;
A /= p[i];
}
ans *= sum(p[i],num*B) % MOD;
ans %= MOD;
}
}
if(A > 1)
{
ans *= sum(A,B) % MOD;
ans %= MOD;
}
cout<<ans<<endl;
}

int main()
{
LL A,B;
isprime();
while(cin>>A>>B)
Solve(A,B);
return 0;
}

第二种方法就是用等比数列求和公式，但是要用逆元。用如下公式即可



因为

可能会很大，超过int范围，所以在快速幂时要二分乘法。

代码：

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 10005;
const int MOD = 9901;

bool prime
;
int p
;
int cnt;

void isprime()
{
cnt = 0;
memset(prime,true,sizeof(prime));
for(int i=2; i<N; i++)
{
if(prime[i])
{
p[cnt++] = i;
for(int j=i+i; j<N; j+=i)
prime[j] = false;
}
}
}

LL multi(LL a,LL b,LL m)
{
LL ans = 0;
a %= m;
while(b)
{
if(b & 1)
{
ans = (ans + a) % m;
b--;
}
b >>= 1;
a = (a + a) % m;
}
return ans;
}

LL quick_mod(LL a,LL b,LL m)
{
LL ans = 1;
a %= m;
while(b)
{
if(b & 1)
{
ans = multi(ans,a,m);
b--;
}
b >>= 1;
a = multi(a,a,m);
}
return ans;
}

void Solve(LL A,LL B)
{
LL ans = 1;
for(int i=0; p[i]*p[i] <= A; i++)
{
if(A % p[i] == 0)
{
int num = 0;
while(A % p[i] == 0)
{
num++;
A /= p[i];
}
LL M = (p[i] - 1) * MOD;
ans *= (quick_mod(p[i],num*B+1,M) + M - 1) / (p[i] - 1);
ans %= MOD;
}
}
if(A > 1)
{
LL M = MOD * (A - 1);
ans *= (quick_mod(A,B+1,M) + M - 1) / (A - 1);
ans %= MOD;
}
cout<<ans<<endl;
}

int main()
{
LL A,B;
isprime();
while(cin>>A>>B)
Solve(A,B);
return 0;
}

其实有些题需要用到

模

的所有逆元，这里

为奇质数。那么如果用快速幂求时间复杂度为

，

如果对于一个1000000级别的素数

，这样做的时间复杂度是很高了。实际上有

的算法，有一个递推式如下



它的推导过程如下，设

，那么



对上式两边同时除

，进一步得到



再把

和

替换掉，最终得到



初始化

，这样就可以通过递推法求出

模奇素数

的所有逆元了。

另外

模

的所有逆元值对应

中所有的数，比如

，那么

对应的逆元是

。