POJ-1845 Sumdiv（唯一分解定理，快速幂，约数和公式，同余模公式）

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 31351		Accepted: 7684

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

题目链接：

http://poj.org/problem?id=1845

题意：

计算A^B 所有约数和s。s%MOD。A和B都很大（0<=A,B<=50000000).

思路：

本题需要用到唯一分解定理，快速幂，约数和公式，同余模公式。

(1).整数的唯一分解定理

任何整数都可以且尽可以分解成若干个素数相乘的形式。

$$A=p{1}^{k1}\ast p{2}^{k2}\ast ......\ast p{n}^{kn}$$

其中的pi各表示一个素数。

(2).约数和公式

按照整数的分解定理，我们有一个整数的约数之和公式为：
$$S=(1+p1+p{1}^2+...+p{1}^{k1})\ast (1+p2+p{2}^2+...+p{2}^{k2})\ast .......\ast (1+pn+p{n}^2+...+p{n}^{kn})$$

还有一个常用的是
​ 因数的个数为(1+k1)(1+k2)(1+k3)…(1+kn)。

(3).同余模公式
$$①(a+b)\%mod=(a\%mod+b\%mod)\%mod;$$

$$②(a\ast b)\%m=(a\%m\ast b\%m)\%m;$$

一、质因数分解

for(int i=2;i*i<=A;) //根号法的体现
{
    int k=0;//p是存素因子的，n是存其指数
    if (A%i==0)
    {
        p[k]=i;
        n[k]=0;
        while(A%i) //找到一个素因子就用此法计算其个数
        {
            n[k]++;
            A/=i;
        }
        k++;
    }
    if(i==2) //除了2其余素数都是奇数哦
         i++;
    else
        i+=2;
}

if(A!=1)//常规来讲，这时候的A已经被分解剩下1了，除非A本身就是素数
{
    p[k]=A;
    n[k++]=1;
}

二、二分法求等比数列的和

解决了A的分解问题，自然需要来解决一下S的求解问题，很明显S是一系列以pi为公比的等比数列和 之积。只要能解决等比数列和的问题那这道题就迎刃而解了啊。最最直接的方法是利用求和公式，但是别忘了，我们的S可是还需要对9901取模的，[pi^(ki*B)-1]/(pi-1)这个结果中pi-1未必和9901互素！因而，解决这个问题，就只好用二分法：

对于一个等比数列求和　　S=1+q+q^2+……qn

如果n为奇数，那么一共就有偶数个项了，

S=[1+q+……+q^(n/2)] * [1+q^(n/2+1)]

如果n为偶数，那么一共就是有奇数个项了，

S=[1+q+……+q^(n/2-1] * [1+q^(n/2+1)]+q(n/2)

*如果记不住上面的公式的话，就举个例子自己算算，结果自然就很清楚了。*

ll sum(ll a, ll b)
{
    if(b == 0)
        return 1;
    if(b & 1) {
        return ((sum(a,b/2) % MOD) * (1+quick_Mod(a,b/2+1))%MOD)%MOD;
    }
    else
        return (((sum(a,b/2-1)%MOD) * (1+quick_Mod(a,b/2+1))%MOD)%MOD + quick_Mod(a,b/2) % MOD)%MOD;

}

代码：

#include<cstdio>
#include <iostream>
#include <cstring>
using namespace std;
typedef  long long ll;
const int maxn = 50000000;
const int MOD = 9901;

int tot;
int p[2000];
int n[2000];
ll quick_Mod(ll a, ll b)
{
    ll ans = 1;
    a %= MOD;
    while(b) {
        if(b & 1) {
            ans = (ans * a) % MOD;
        }
        a = (a * a) % MOD;
        b >>= 1;
    }
    return ans % MOD;
}
ll sum(ll a, ll b)
{
    if(b == 0)
        return 1;
    if(b & 1) {
        return ((sum(a,b/2) % MOD) * (1+quick_Mod(a,b/2+1))%MOD)%MOD;
    }
    else
        return (((sum(a,b/2-1)%MOD) * (1+quick_Mod(a,b/2+1))%MOD)%MOD + quick_Mod(a,b/2) % MOD)%MOD;

}
ll solve(ll a, ll b)
{
    tot = 0;

    memset(p, 0 , sizeof(p));
    for(int i = 2; i * i <= a; i++) {
        if(a % i == 0) {
            n[tot] = i;
            while(a % i == 0) {
                p[tot]++;
                a /= i;
            }
            tot++;
        }
    }
    if(a != 1) {
        n[tot] = a;
        p[tot++] = 1;
    }
    ll ans = 1;
    for(int i = 0; i < tot; i++) {
        ans = (ans * (sum(n[i], p[i]*b))%MOD) % MOD;
    }
    return ans % MOD;
}
int main()
{
    ll a, b;
    while(cin >> a >> b) {
        cout << solve(a,b)%MOD << endl;
    }


    return 0;
}