Skip to content

算法刷题【洛谷P1593】因子和(附等比数列求和公式推导)

信息学竞赛

2022-10-25

题目描述

输入两个整数 aabb,求 aba^b 的因子和。

由于结果太大,只要输出它对 99019901 取模的结果。

输入格式

仅一行,为两个整数 aabb

输出格式

输出一行一个整数表示答案对 99019901 取模的结果。

输入输出样例

In 1:

2 3

Out 1:

15

数据范围

对于全部的测试点,保证 1a5×1071 \leq a \leq 5 \times 10^70b5×1070 \leq b \leq 5 \times 10^7

题解

本文视频讲解:

[video(video-Dw9oUPTm-1643561388309)(type-csdn)(url-https://live.csdn.net/v/embed/185217)(image-https://live-file.csdnimg.cn/release/live/file/1643560899980.png?x-oss-process=image/resize,l_300)(title-【洛谷P1593】因子和视频讲解)]

为了完成洛谷P1593这道题也是拼了……用到了三个不会的知识:因子和公式等比数列求和公式费马小定理求分数

这里另摆一份很好的题解:洛谷 [P1593 因子和] {快速幂+费马小定理求逆元+求解质因子} 奋斗的珂珂~ - 代码先锋网


先摆出因子和公式的参考资料:


等比数列定义很简单,设有数列 [a1,a2,...,an][a_1, a_2,...,a_n] ,存在一个定值 qq 使得任意 1in11 ≤ i ≤ n-1 都有 ai+1=ai×qa_{i+1} = a_i \times q 。由此定义可知 an=a1×qn1a_n = a_1 \times q^{n-1}

这个数列的求和公式为 Sn=a1+a2+...+an=a1qn1q1S_n = a_1 + a_2 + ... + a_n = a_1 \frac{q^n-1}{q-1} ,下面我们来证明这个公式:

qSn=qa1+qa2+...+qan=a2+a3+...+an+1∵ qS_n = qa_1 + qa_2 + ... + qa_n = a_2 + a_3 + ... + a_{n+1}

qSnSn=(q1)Sn=an+1a1=(a1qn)a1=a1(qn1)∴ qS_n - S_n= (q-1) S_n = a_{n+1} - a_1 = (a_1 * q^n) - a_1 = a_1(q^n-1)

Sn=a1qn1q1∴ S_n = a_1 \frac{q^n-1}{q-1}

(参考资料:等比数列公式及推导_高三网


费马小定理求分数快速幂的公式及证明:

(来源:分数取模(费马小定理)_give it a try-CSDN博客_分数取模运算


根据整数的唯一分解定理,整数a进行质因数分解对应的式子唯一,有:

a=p1k1p2k2p3k3pnkna = p_1^{k_1} * p_2^{k_2} *p_3^{k_3}* … * p_n^{k_n}

又因为本题要分解的是aba^b,所以上面的式子又可以写成这样:

ab=p1k1bp2k2bp3k3bpnknba^b= p_1^{k_1*b} * p_2^{k_2*b} *p_3^{k_3*b}* … * p_n^{k_n*b}

证明很简单,就是把上面第一个式子乘上 bb 次即可得第二个式子

接下来我们要求的是因子和,所以就有:

ans=(1+p11+p12+p13++p1k1b)(1+p21+p22+p23++p1k2b)...(1+pn1+pn2+pn3++pnknb)ans= (1+p_1^1 + p_1^2 +p_1^3+ … + p_1^{k_1*b})*(1+p_2^1 + p_2^2 +p_2^3+ … + p_1^{k_2*b})*...*(1+p_n^1 + p_n^2 +p_n^3+ … + p_n^{k_n*b})

于是求和转化成了等比数列的和的乘积

而对于每一个等比数列,可根据公式求和:sum=pikib+11pi1sum=\frac{p_i^{k_i*b+1}-1}{p_i-1} (等比数列求和公式中代入 a1=1a_1 = 1

注意这里的 pp 与上面截图中的 pp 含义不一样!这里指的是公比(类比等差数列的公差),上面指的是 mod (快速幂中取余的那个数)

有了这些知识来看代码:

#include <bits/stdc++.h>
using namespace std;
const int mod = 9901;
int ys[50000000];

// 一个不那么简单的快速幂模板
long long ksm(long long x, long long y) {
    if (y == 0) return 1;
    if (y % 2 == 0) return ((long long)pow(ksm(x, y / 2) % mod, 2) % mod);
    return (x * (long long)pow(ksm(x, y / 2) % mod, 2) % mod);
}

int main() {
    int a, b;
    cin >> a >> b;
    int aa = a;  // 由于aa的值会被修改,所以需要一个变量保存

    ys[2] = 0;
    while (a % 2 == 0) {
        a /= 2;
        ys[2]++;
    }
    ys[2] *= b;  // 因为是b次幂因此要乘以b
    for (int i = 3; i <= a / i; i += 2) {
        ys[i] = 0;
        while (a % i == 0) {
            a /= i;
            ys[i]++;
        }
        ys[i] *= b;
    }
    if(a != 1) {
        // 分解质因数,若有质因数超过根号a,则只能是a本身
        ys[a] = b;
    }

    long long ans = 1;
    for (int i = 2; i <= aa; i++) {
        if (!ys[i]) continue;  // 不存在该质因数
        int temp;
        if (i % mod == 1) {
            temp = (ys[i] + 1) % mod;  // 详见博文最后的解释
        } else
            temp = ((ksm(i, ys[i] + 1) - 1 + mod) % mod) * ksm(i - 1, mod - 2) % mod;  // 使用费马小定理求解
        if (temp != 0) ans = ans * temp % mod;
    }
    cout << ans << endl;

    return 0;
}

关于代码第42行特判 i % mod == 1 (也就是判断 i-1 是否为 mod 的倍数;同样适用于45行 temp != 0 特判):

i-1mod 的倍数时,等比数列中每一项取余 99019901 的结果都为 11 ,等比数列的和便是 ki×b(上方证明公式中的表示)=ys[i](代码中的表示)k_i \times b\text{(上方证明公式中的表示)} = ys[i]\text{(代码中的表示)} ,因此此处特判改为 temp = ys[i] + 1 ;同理,若 imod 的倍数,则 temp 值为 00 恒成立(洛谷测试点中不存在此极端情况因此不考虑也能AC(又更新:目前洛谷已加强))

Ukraine 在俄罗斯对乌克兰发动的野蛮的侵略战争中矢志不渝地支持乌克兰