AtCoder が提供しているABC（AtCoder Beginner Contest）293 のE問題をC++とPythonで解いてみました。ABC293は、2023年3月11日21:00に実施されました。
AtCoder の紹介はこちらに、プログラミングの方針はこちらに記事があります。
E問題 Geometric Progression（Difficulty : 1453）
問題はリンク先をご覧ください。
ABC293 E問題 Geometric Progression
ダブリングを使って解くことができました。AtCoder Problems による Difficulty は 1453 でした。
考察
この記事は、公式解説で紹介されている twitter解説の「E問題はダブリングなのだ！」を参考にしました。
X が 2 のべき乗の場合
$X$ が $2^i$ の場合は、以下の手順で値 dp[i] を求めることができます。なお、適当なタイミングで M で割った余りを求めることにして、M で割る式は省略します。
- i = 0: dp[0] = $A^0 = 1$
- i = 1: dp[1] = $1 + A^1 \times 1$ = dp[0] + $A^1 \times$ dp[0]
- i = 2: dp[2] = $1 + A^1 + A^2(1 + A^1)$ = dp[1] + $A^2 \times$ dp[1]
- i = 3: dp[3] = $1 + A^1 + A^2 + A^3 + A^4(1 + A^1 + A^2 + A^3)$
= dp[2] + $A^4 \times$ dp[2]
まとめると、$X = 2^i$ の場合、以下の漸化式で求めることができます。
dp[i] = dp[i – 1] + $A^{2^{i-1}}\times$ dp[i – 1]
X が一般的な数の場合
次に一般的な $X$ について、値を求めます。
例えば、$X = 11 = 2^3 + 3$ について、次のように求めることができます。
$$\sum^{10}_{i = 0} A^i = 1 + A + \cdots + A^{10}$$
$= 1 + A + A^2 + A^3(1 + A + \cdots + A^7) = 1 + A + A^2 + A^3 \times$ dp[3]
前半は、以下のように書き換えることができます。
$= 1 + A + A^2 = 1 + A(1 + A) = 1 + A \times$ dp[1]
つまり、$X$ が $2^i$ の場合に値を求めておけば、$X = 11$ の場合は、$X = 8 + 2 + 1$ と分解して値を求めることができます。このように 2 のべき乗の場合を先に計算しておいて、値を求める手法をダブリングと呼ぶようです。
解答案
C++ プログラム例（ABC293E）
剰余付きの$N$乗を求めるために、繰り返し二乗法と呼ばれる手法を使いました。具体的な実装は、関数 power を参照してください。
上で考察したように2のべき乗に対して、配列 d[i]（i = 0, …, 60）を求めておきます（23から28行目）。$X$ に対して、一番大きなビットから処理をしていきます。処理しているビットより小さな数は、変数 temp_p として求めて、残りの処理を行います。
以下が、C＋＋プログラムとなります。
#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long int ull;
int power(int n, ull p, int mod) {
if (p == 0) {
return 1;
}
if (p % 2 == 0) {
long long int t = power(n, p / 2, mod);
return (int)(t * t % mod);
}
return (int)(((long long int)n * power(n, p - 1, mod))% mod);
}
int main()
{
ull a, x, m;
cin >> a >> x >> m;
vector<ull> dp(61);
ull p = a;
dp[0] = 1;
for (int i = 1; i <= 60; ++i) {
dp[i] = (dp[i - 1] + dp[i - 1] * p)% m;
p = (p * p)% m;
}
ull result = 0;
for (int i = 60; i >= 0; --i) {
if (x & (1ULL << i)) {
ull temp_p = x & ((1ULL << i) - 1);
result += (dp[i] * power(a, temp_p, m))% m;
result %= m;
}
}
cout << result << endl;
return 0;
}AC（Accepted＝正しいプログラム）と判定されました。
Python プログラム例（ABC293E）
Python では、剰余付きのべき乗を求める関数 pow が提供されています。基本的な考え方は、C++ 版を移植しました。
やはり、Python ではすっきりと書ける印象です。
"""AtCoder Beginner Contest 293 E"""
a, x, m = map(int, input().split())
dp = [0] * 61
dp[0] = 1
p = a
for i in range(1, 61):
dp[i] = (dp[i - 1] + dp[i - 1] * p) % m
p = p * p % m
result = 0
for i in range(60, -1, -1):
if x & (1 << i):
temp_p = x & ((1 << i) - 1)
result += dp[i] * pow(a, temp_p, m)
result %= m
print(result)こちらも「AC」と判定されました。
最後に
問題を読んだときに、M が素数であれば、$(A^x-1)/(A-1)$ を求めれば良いことが分かりました。しかし、M が素数ではない場合に計算する工夫が分かりませんでした。公式解説では、このような場合も適用できる計算を紹介していました。この方法では非常にすっきりと書けています。
ただし、この方法は自分で思いつくことは難しく、また数論を学び直す気力もわかないため、「タブリング」というヒントの言葉から、プログラムを書いてみました。
引き続き ABC の問題を紹介していきます。
