POJ 1845.Sumdiv

最新推荐文章于 2022-02-26 10:03:02 发布

Pegasi_Tio

最新推荐文章于 2022-02-26 10:03:02 发布

阅读量144

点赞数

分类专栏： POJ ACM

本文链接：https://blog.youkuaiyun.com/PegasiTIO/article/details/89249456

版权

ACM 同时被 2 个专栏收录

115 篇文章

订阅专栏

POJ

17 篇文章

订阅专栏

本文介绍了一种算法，用于计算A^B的所有自然数约数之和并对9901取模。通过质因数分解和分治法，实现了高效求解。代码中包含了快速幂、质数筛选等关键步骤。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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

输入

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

输出

The only line of the output will contain S modulo 9901.

样例输入

2 3

样例输出

15

提示

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

求A^B的所有约数之和并mod9901

把A分解

根据唯一分解定理(质因数分解)：任何一个大于1的正整数都能唯一分解成为有限个质数的乘积，可写作：

N=p₁^c₁p₂^c₂…p_m^c_m

其中 $c_i$ 都是正整数， $p_i$ 都是质数，且满足 $p_1<p_2<···<p_m$

如果正整数 $N$ 是合数，则存在能整除 $N$ 的数 $T$ ，其中 $\le T \le \sqrt N$

表示 $A^B$

A^B=p₁^B*C₁p₂^B*C₂ ···p_n ^B*C_n

$A^B$ 的所有约数可以表示成集合{p₁^k₁p₂^k₂…p_m^k_n}，其中 $\le k_i \le B*C_i(1 \le i \le n)$

根乘分配律，所有约数之和为：

(1+p₁+p₁²+···+p₁^B*c₁)*(1+p₂+p₂²+···+p₂^B*c₂)*···*(1+p_n+p_n²+···+p_n^B*c_n)

不能利用等比数列求和，因为要对结果求余，而求余运算不对除法具有分配律

利用分治法求和： $sum(p,c)=1+p+p_2+···+p_c$

若c为奇数：

sum(p,c) =(1+p+···+p^(c-1)/2 )+(p^(c+1)/2+···+p^c)

=(1+p+···+p^(c-1)/2 )+p^(c+1)/2*(1+p+···+p^(c-1)/2 )

=(1+p^(c+1)/2)*sum(p,(c-1)/2)

若c为偶数,同理：

sum(p,c)=(1+p^c/2)*sum(p,c/2-1)+p^c

计算乘方时使用快速幂，这样这题就解完了

代码实现:

先计算sqrt(50000000) ，也就是能整除50000000的数最大不会超过多少，结果是7071.07我们只需要计算10000以内的素数即可，还有要非常注意的一点是，整除不代表能除完，也就是最后可能还会剩下一个数，这个数比10000以内最大的素数大，并且也是素数，所以要注意判断分解质因数后是否为1

例：

还想吐槽一下markdown里LaTeX的语法
$p_1^(c_1)$ ???我就不能在上标里打个带下标的嘛，还专门加了个括号，也可能用法不对吧，学了学html的语法打上去了

#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;
static const auto io_sync_off = []() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    return nullptr;
}();

using ll = long long;
using paii = pair<int, int>;

const int mod = 9901;
const int maxn = 10005;
bool is_not_prime[maxn];
int prime[maxn], tot = 0;

void getPrime() // 素数表
{
    is_not_prime[1] = 1;
    for (int i = 2; i <= maxn; ++i)
    {
        if (!is_not_prime[i])
            prime[++tot] = i;
        for (int j = 1; i <= tot; ++j)
        {
            is_not_prime[i * prime[j]] = true;
            if (i % prime[j] == 0)
                break;
        }
    }
}

ll power(ll a, ll b) // 快速幂
{
    int ans = 1;
    for (; b; b >>= 1)
    {
        if (b & 1)
            ans = ans * a % mod;
        a = a * a % mod;
    }
    return ans;
}

vector<paii> pfn;
void pf(int x) // Prime factorization分解质因数
{
    for (int i = 1; i <= tot; ++i)
    {
        int cnt = 0;
        while (x % prime[i] == 0)
        {
            x /= prime[i];
            ++cnt;
        }
        if (cnt)
            pfn.push_back(make_pair(prime[i], cnt));
        if (x == 1)
            return;
    }
    if (x != 1) // 重要！！！
        pfn.push_back(make_pair(x, 1));
}

ll sumdiv(ll p, ll c)
{
    if (c == 0) // 返回条件
        return 1;
    if (c & 1)
        return (1 + power(p, (c + 1) >> 1)) * sumdiv(p, (c - 1) >> 1) % mod;
    else
        return (1 + power(p, c >> 1)) * sumdiv(p, (c >> 1) - 1) % mod + power(p, c) % mod;
}

int main()
{
    // cout<<sqrt(50000000)<<endl;
    getPrime();
    int a, b;
    cin >> a >> b;
    pf(a);

    ll ans = 1;
    for (int i = 0; i < pfn.size(); ++i)
        ans = (ans * sumdiv(pfn[i].first, b * pfn[i].second)) % mod;
    cout << ans;
    return 0;
}