POJ 1845 Sumdiv （唯一分解定理+求等比数列前n项和）

最新推荐文章于 2021-11-03 14:46:45 发布

原创最新推荐文章于 2021-11-03 14:46:45 发布 · 913 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#数论

ACM 专栏收录该内容

59 篇文章

订阅专栏

本文介绍一种高效算法来解决大数幂次的因数之和问题，并给出具体实现代码。通过质因数分解及等比数列求和原理，结合模运算公式，实现了对大数A^B的因数之和求模运算。

题目大意：先求出A^B（0 <= A,B <= 50000000）的因数之和sum，然后求sum % 9901的值。

由于A，B的取值范围很大，所以直接运算肯定不行。

唯一分解定理

任何一个大于1的自然数 N,如果N不为质数，那么N可以唯一分解成有限个质数的乘积

所以 A = p1^q1 * p2^q2 * …… * pn^qn（pi为素数）
从而 A^B = (p1^q1 * p2^q2 * …… * pn^qn) ^ B
= p1^(q1*B) * p2^(q2*B) * …… * pn^(qn*B)
由于pi是素数，很容易推得pi^(qi*B)的因数之和为

    1 + pi + pi^2 + ...... + pi^(qi*B)

从而A^B的因数之和为

(1 + p1 + p1^2 + ...... + p1^(q1*B)) * (1 + p2 + p2^2 + ...... + p2^(q2*B)) * ...... * (1 + pn + pn^2 + ...... + pn^(qn*B))

如果先求出这整个的和值，肯定超long long，所以要边乘边取模。

模运算公式

(a + b) % p = (a % p + b % p) % p
(a - b) % p = (a % p - b % p) % p
(a * b) % p = (a % p * b % p) % p
(a ^ b) % p = ((a % p)^b) % p

另外由于qi * B的值很大，一步步算肯定会超时，注意到每一项都是一个等比数列的和。

求等比数列前N项和

求表达式等比数列的值

当时，
时，则有
时，有

由以上分析，用分治法求等比数列前N项和。
参考：http://blog.youkuaiyun.com/ACdreamers/article/details/7851144

代码

#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>

using namespace std;
typedef long long ll;
const int maxn = 10000;
ll p[maxn] , q[maxn];

ll pow_mod(ll a , ll n , ll m)
{
    if(n == 0) return 1;
    ll x = pow_mod(a , n >> 1 , m);
    ll ret = x * x % m;
    if(n & 1) ret = ret * a % m;
    return ret;
}

//求(1 + a + a^2 + ... + a^n) % m
ll sum(ll a , ll n , ll m)
{
    if(!n) return 1;
    if(n == 1) return a % m;
    ll ret = 0;
    if(n&1) {
        ret = (1 + pow_mod(a , (n - 1) / 2 + 1 , m)) % m;
        ret = (ret * sum(a , (n - 1) / 2 , m)) % m;
        ret = (ret + pow_mod(a , (n - 1) / 2 + 1 , m)) % m;
    }
    else {
        ret = (1 + pow_mod(a , n / 2 , m)) % m;
        ret = (ret * sum(a , n / 2 , m)) % m;
    }
    return ret;
}
int main()
{
    int a , b;
    while(cin >> a >> b)
    {
        if(!a) {cout << 0 << endl; continue;}
        if(!b) {cout << 1 << endl; continue;}
        int n = 0;
        // 根号法 + 递归法 分解整数A
        for(int i = 2; i * i <= a && a != 1; i++) if(a % i == 0)   
        {
            p[n] = i , q[n] = 0;
            while(a % i == 0) {q[n]++; a /= i;}
            n++;
            if(i != 2) i++; //奇偶法
        }
        if(a != 1) {p[n] = a; q[n++]=1;}

        ll ans = 1 , mod = 9901;
        for(int i = 0; i < n; i++) {
            ll tmp = (sum(p[i] , q[i] * b , mod) + 1) % mod;
            ans = ans * tmp % mod;
        }
        cout << ans << endl;
    }
    return 0;
}