【OJ】POJ1845 （数学, a*b%mod, a^b%mod, a/b%mod, 分治, 逆元）

题目链接

Language:Default Sumdiv Time Limit: 1000MS Memory Limit: 30000K Total Submissions: 28200 Accepted: 6939 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=p_1^{a_1}\times p_2^{a_2}\times p_3^{a_3}\times ...\times p_n^{a_n}$$
$$A^B=p_1^{B\times a_1}\times p_2^{B\times a_2}\times p_3^{B\times a_3}\times ...\times p_n^{B\times a_n}$$
$$\begin{gathered} A^B=(1+p_1+p_1^2+p_1^3+...+p_1^{B\times a_1})\times (1+p_2+p_2^2+p_2^3+...+p_2^{B\times a_2}) \\ \times ...\times (1+p_n+p_n^2+p_n^3+...+p_n^{B\times a_n}) \end{gathered}$$
$$设sum(p,n)=1+p+p^2+p^3+...+p^{n-1}$$
$$\begin{gathered} （法一：分治法） \\ sum(p,n)=(1+p^{\lfloor n/2\rfloor })\times sum(p,\lfloor n/2\rfloor )+(n\mathrm{m}\mathrm{o}\mathrm{d}2)\times p^{n-1} \end{gathered}$$

LL sum(LL q, LL n, LL mod)
{
    if(n == 0)return 0ull;
    if(n == 1)return 1ull;
    LL ans = sum(q, n >> 1, mod), tmp = (1ull + pow_mod(q, n >> 1, mod)) % mod;
    if(1u & n)
        return ans = (mul_mod(ans, tmp, mod) + pow_mod(q, n - 1, mod)) % mod;
    else
        return ans = mul_mod(ans, tmp, mod);
}

$$\begin{gathered} (法二：除法取模，逆元) \\ M=m\times (p-1) \\ sum(p,n)\mathrm{m}\mathrm{o}\mathrm{d}m=((p^n\mathrm{m}\mathrm{o}\mathrm{d}M-1)\mathrm{m}\mathrm{o}\mathrm{d}M)/(p-1) \end{gathered}$$

LL sum(LL q, LL n, LL mod)
{
    LL new_mod = (q - 1) * mod;
    LL numerator = pow_mod(q, n, new_mod) - 1;
    numerator = (numerator % new_mod + new_mod) % new_mod;
    return numerator / (q - 1);
}

主要公式在上面已经列举。快速幂，加减乘取模这种模板就不用多说了。另外还有一个质因数分解的模板，这个可以借助stl的vector+pair或者map。

代码

#include <iostream>
#include <vector>
#include <utility>

typedef long long LL;
using namespace std;

const LL mod = 9901;

LL mul_mod(LL a, LL b, LL mod)
{
    LL ans = 0;
    for(; b; b >>= 1)
    {
        if(1 & b)ans = (ans + a) % mod;
        a = (a << 1) % mod;
    }
    return ans;
}

//快速幂a^b,对mod取模,这是位运算版，且只适合非负数
LL pow_mod(LL a, LL b, LL mod)
{
    LL ans = 1 % mod;
    for(; b; b >>= 1)
    {
        if(1 & b)ans = mul_mod(ans, a, mod);
        a = mul_mod(a, a, mod);
    }
    return ans;
}

//质因数分解
//x>=2
vector<pair<LL,LL> > get_prime_factors(LL x)
{
    vector<pair<LL, LL> > prime_factors;
    for(LL i = 2; i*i<=x; i++)
    {
        if(x % i)continue;
        LL cnt = 0;
        while(x % i == 0)
        {
            x /= i;
            cnt++;
        }
        prime_factors.push_back(make_pair(i, cnt));
    }
    //注意，最后可能会留一个大于n^(1/2)的质数，n为x原本的大小
    //例如14=2*7
    if(x!=1)prime_factors.push_back(make_pair(x,1));
    return prime_factors;
}

/*
//一种map版的写法
//
//map使用[]时，复杂度为O(logn)
//如果没有元素，会先自动生成一个二元组
//自动生成的value为zero（0，""等）

map<LL,LL> prime_factor(LL n)
{
    map<LL,LL> res;
    for(LL i=2;i*i<=n;i++)
        while(n%i==0) {++res[i];n/=i;}
    if(n!=1) res[n]=1;
    return res;
}
*/

等比数列求和，对mod取模

//除法取模版（或者说时逆元版，虽然没有逆元）
//1+q+...+q^(n-1)=(q^n-1)/(q-1)
//(a/b)%m = (a%(b*m))/b
LL sum(LL q, LL n, LL mod)
{
    LL new_mod = (q - 1) * mod;
    LL numerator = pow_mod(q, n, new_mod) - 1;
    numerator = (numerator % new_mod + new_mod) % new_mod;
    return numerator / (q - 1);
}

/*
//二分版
//(由于mod只对加、减、乘有分配率，所以涉及到取模的，不能直接用等比数列求和公式)
//首项a0，公比q，项数n，模mod，不过这里a0就是1，所以不用管
LL sum(LL q, LL n, LL mod)
{
    if(n == 0)return 0ull;
    if(n == 1)return 1ull;
    LL ans = sum(q, n >> 1, mod), tmp = (1ull + pow_mod(q, n >> 1, mod)) % mod;
    if(1u & n)
        return ans = (mul_mod(ans, tmp, mod) + pow_mod(q, n - 1, mod)) % mod;
    else
        return ans = mul_mod(ans, tmp, mod);
}
*/

LL get_ans(LL A, LL B, LL mod)
{
    if(A == 0)return 0;
    else if(A == 1 or B == 0)return 1 % mod;
    vector<pair<LL,LL> > prime_factors = get_prime_factors(A);
    LL ans = 1ll;
    for(vector<pair<LL,LL> >::iterator it = prime_factors.begin(); it != prime_factors.end(); it++)
    {
        LL q = it->first, c = it->second;
        LL tmp = sum(q, c * B + 1, mod);
        ans = mul_mod(ans, tmp, mod);
    }
    return ans;
}

int main()
{
    LL A, B;
    cin >> A >> B;
    LL ans = get_ans(A, B, mod);
    cout << ans << endl;
    return 0;
}