【PTA】7-2 国王的奖励——分数取模、快速幂、int64的乘法模运算【C/C++】

• 1 问题
• • 1.1 题面描述
  • 1.2 输入描述
  • 1.3 输出描述
  • 1.4 样例描述
  • 1.5 样例解释
• 2 分析
• • 2.1 数学抽象
  • 2.2 解决方法
  • • 2.2.1 分数取模——乘法逆元
    • 2.2.2 求解 q^n^ ——快速幂
    • 2.2.3 int64 的乘法模运算——位运算法
• 3 结果
• 4 参考文章

7-2 国王的奖励

1.1 题面描述

国际象棋是很久以前由一个印度人 Shashi 发明的，当他把该发明献给国王的时候，国王很高兴，就许诺可以给这个发明人任何他想要的奖赏。Shashi 要求以这种方式给他一些粮食：棋盘的第 1 个方格里只放 1 粒麦粒，第 2 格 q 个,第 3 格 q² 个,第 4 格 q³ 个……，直到 n 个格子全部放满。这个奖赏最终会是什么样子的呢？Shashi 已经有算法了，请你算一算吧。
当然 Shashi 并不关心具体的数有多少了，他有一个检验你的答案是否与他心意相通的办法:把你求出的答案对 100000007 取模看看和 Shashi 算的是否一样就行了。

1.2 输入描述

本题输入具有多组样例。
第一行为一个数 t (1 ≤ t ≤ 500)，代表测试数据的组数。
之后的 t 行，每行两个数 q, n (1 ≤ q, n ≤ 10¹⁸)，含义见题目描述。

1.3 输出描述

输出共 t 行，每行一个数，为你算的答案。
友情提醒：你算的答案应该比 100000007 小。

1.4 样例描述

样例输入：

3
2 1
2 2
2 3

样例输出：

1
3
7

1.5 样例解释

对于样例：

• q = 2, n = 1 时仅有1个麦粒
• q = 2, n = 2 时有 1 + 2 = 3 个麦粒
• q = 2, n = 3 时有 1 + 2 + 2² = 7 个麦粒

作者｜cugbacm
单位｜中国地质大学（北京）
代码长度限制｜16 KB
时间限制｜1000 ms
内存限制｜256 MB

该问题本质为等比数列求和问题，q 为公比，n 为项数。
需要注意：

• 每次输入的 q, n 可能不同；
• q, n 的范围超过 int32，需要用 int64 来存储；
• 结果要对 100000007 取模。

等比数列求和方法一般有直接求和、二分递归与等比求和通项公式，我们在本文中主要讨论等比求和公式法，二分递归可参见参考文章1。

等比数列求和通项公式如下：
S(n) = n (q = 1)
S(n)=q₀ * (1 - qⁿ) / (1 - q) (q ≠ 1)

加和部分时间复杂度只有 O(1)，主要取决于求解 qⁿ 的算法。
但模运算对除法不存在以下关系：
(a mod p) / (b mod p) = (a / b) mod p
不能直接使用通项公式。

对分数取模有两种方法：乘法逆元和变换模值，我们在此仅对逆元进行讨论，变换模值可参见参考文章1。

若在 mod p 意义下，对于一个整数 a，有 a * b ≡ 1 (mod p)，那么这个整数 b 与 a 互为乘法逆元。

求乘法逆元的方法主要有：将_a * b ≡ 1 (mod p)_ 转化为 a * x + p * y = 1，用拓展欧几里得算法解二元一次方程解出 x。
但本题中 p = 100000007 是质数，因此可以用小费马定理（详见参考文章2）直接求解。

小费马定理：
a^p-1 mod p = 1，p 为质数

因此有
(a / b) mod p = (b^p-1 mod p) * ((a / b) mod p) = a * b^p-2 mod p。
即
S(n) = n mod p (q = 1)
S(n)=q₀ * (1 - qⁿ) * (1 - q)^p-2 mod p (q ≠ 1)

注意若 b 是 p 的倍数则不能用逆元，因为此时 b mod p ≡ 0，无法求解逆元。

代码如下：

// q: 公比, n：项数, mod: 模, qpow: 求幂函数
if (q == 1) {
	printf("%lldn", n % mod);
}
else {
	printf("%lldn", (qpow(q, n) - 1) * qpow(q - 1, mod - 2) % mod);
}

解决了等比求和通项公式的问题，只剩求解 qⁿ 了，此处选择了位运算的快速幂算法。

代码如下：

// mod: 模, qmul: 求积函数
long long qpow(long long base, long long power) {
	base %= mod;			// 对指数的底先取模，减少运算次数
	long long result = 1;
	while (power > 0) {
		if (power & 1) {	// power 为奇数时
			result = qmul(result, base);
		}
		base = qmul(base, base);
		power >>= 1;		// 即 power /= 2
	}
	return result;
}

这部分优化直接借用参考文章3，将大数乘法的时间复杂度从反复翻倍法（按位相乘）的 O(log₂n) 降到了 O(1)！！

代码如下：

// mod: 模
long long qmul(long long a, long long b) {
	// 前两句处理了负数和数值过大的情况
	a = (a % mod + mod) % mod;
	b = (b % mod + mod) % mod;
	long long c = a * (long double) b / mod;
	long long ans = a * b - c * mod;
	// 若结果不在 [0, p) 区间则调整一下
	if (ans < 0) {
		ans += mod;
	}
	else if (ans >= mod) {
		ans -= mod;
	}
	return ans;
}

完整代码如下：

#include
#define mod 100000007

long long qmul(long long a, long long b) {
	a = (a % mod + mod) % mod;
	b = (b % mod + mod) % mod;
	long long c = a * (long double) b / mod;
	long long ans = a * b - c * mod;
	if (ans < 0) {
		ans += mod;
	}
	else if (ans >= mod) {
		ans -= mod;
	}
	return ans;
}

long long qpow(long long base, long long power) {
	base %= mod;
	long long result = 1;
	while (power) {
		if (power & 1) {
			result = qmul(result, base);
		}
		base = qmul(base, base);
		power >>= 1;
	}
	return result;
}

int main() {
	int t;
	scanf("%d", &t);
	long long q, n;
	for (int i = 0; i < t; i++) {
		scanf("%lld%lld", &q, &n);
		if (q == 1) {
			printf("%lldn", n % mod);
		}
		else {
			printf("%lldn", (qpow(q, n) - 1) * qpow(q - 1, mod - 2) % mod);
		}
	}
	return 0;
}

提交结果如下：
7ms，低到惊人的耗时！

1. 《poj 1845 Sumdiv 数论–等比数列和(逆元或者递归)》
2. 《费马小定理（易懂）》
3. 《小技巧1——长整型：64位整数的乘法模运算》