Sumdiv POJ - 1845 [唯一分解+约数和定理+二分求公式]

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 25921		Accepted: 6405

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

题目：求A^B的约数和S%9901。

分析：首先题目所需知识：唯一分解定理、约数和定理、二分思想。

唯一分解定理：

任何一个数都能由p[0]^e[0] * p[1]^e[1] * p[2]^e[2] *...*p[n]^e[n]唯一表示（p[0],p[1],...,p[n]都是素数），即A = p[0]^e[0] * p[1]^e[1] * p[2]^e[2] *...*p[n]^e[n];

约数和定理：

A的所有约数和就是（1 + p[0]^1 + ... + p[0]^e[0]）*（1 + p[1]^1 + ... + p[1]^e[1]）*...* （1 + p[n]^1 + ... + p[n]^e[n]）;

约数和定理也是很显然，大胆猜想！！！！

那么我们S就可以求啦？？？不对9901不是素数啊，不能用等比数列公式，这下怎么办？？？

下面我们用二分的思想把公式求出来。。。

1 + p^1 + ... + p^n :

如果n是奇数那么直接(1 + p^(n/2 + 1)) *(1 + ... + p^(n/2));

如果n是偶数那么就等于(1 + p^((n-1)/2 + 1))*(1 + ... + p^((n-1)/2)) + p^n;

这样子递归就能求解啦！！！

附上ac代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define ll long long

using namespace std;
const ll mod = 9901;
const int maxn = 5e4+7;
ll a, b;
int prime[maxn], tot = 0;
bool p[maxn*100];

void getprime()
{
    memset(p, 0, sizeof(p));
    for(int i = 2; i < maxn; i++) {
        if(!p[i]) prime[tot++] = i;
        for(int j = 0; j < tot&&i * prime[j] < maxn; j++)
        {
            p[i*prime[j]] = true;
            if(i%prime[j] == 0) break;
        }
    }
}

ll pow(ll x, ll y)
{
    ll ans = 1;
    while(y)
    {
        if(y&1) ans = ans*x%mod;
        x = x*x%mod;
        y >>= 1;
    }
    return ans;
}



ll fd(ll x, ll y){
    if(y == 0) return 1;
    if(y == 1) return (1+x)%mod;
    if(y%2 == 0) return (((pow(x, y)%mod)) + (fd(x, y-1)%mod))%mod;
    else return (((1+pow(x, y/2+1))%mod)*(fd(x, y/2)%mod))%mod;
}


ll wf(ll x, ll y)
{
    ll ans = 1, k = x;
    for(int i = 0; i < tot&&k>1&&prime[i]*prime[i] <= k; i++)
    {
        if(k%prime[i] == 0)
        {
            int c = 0;
            while(k%prime[i] == 0) {
                k /= prime[i];
                c++;
            }
            c*=y;
            ans = ans*fd(prime[i], c)%mod;
        }
    }
    if(k > 1) ans = ans*fd(k, y)%mod;
    return ans;
}


int main()
{
    getprime();
    //freopen("in.txt", "r", stdin);
    while(~scanf("%lld%lld",&a, &b))
    {
        printf("%lld\n", wf(a, b));
    }
    return 0;
}