AtCoder が提供しているABC（AtCoder Beginner Contest）300 のD問題をC++とPythonで解いてみました。ABC300は、2023年4月29日21:00に実施されました。
AtCoder の紹介はこちらに、プログラミングの方針はこちらに記事があります。
D問題 AABCC（Difficulty : 908）
問題はリンク先をご覧ください。
素数表を計算して、条件を満たす素数の組み合わせをカウントします。AtCoder Problems による Difficulty は 908 でした。
解答案
3つの素数 $a < b < c$ に対して、$a^2 \times b \times c^2 \leqq N$ となる $a, b, c$ の組み合わせを計算します。$N$ の上限は、$10^{12}$ です。
$N = 10^{12}$ の場合を考えます。
$a = 2, b = 3$ のときに、$c$ の範囲が一番大きくなります。具体的には、$c$ の上限は以下となります。
$$c = \sqrt{\frac{10^{12}}{2^2 \times 3}} = 288675.1345 \ldots$$
キリが良いので、30万以下の素数表を作ります。$b$ と $a$ の上限も求めてみます。$b$ の上限は、$a = 2$ のときに以下となります。
$$b = \sqrt[3]{\frac{10^{12}}{2^2}} = 6299.6052 \ldots$$
$a$ の上限は、$10^{12}$ の5乗根となります。
$$a = \sqrt[5]{10^{12}} = 251.1866 \ldots$$
ループ回数を確認します。
- a の範囲の素数は 55 個
- b の範囲の素数は 820 個
- c の範囲（30万以下の素数の数）の素数は 25997 個
55 × 820 × 25997 = 1172464700 ですが、c は条件を満たさなくなればループを抜けるため、この回数のループは実行されません。事実、入力例2が上限の計算結果を示しており、ループは、2817785 回しか繰り返されていません。
泥臭い方法ですが、30万までの素数表を計算して、上記で求めた範囲で $a^2 \times b \times c^2 \leqq N$ となる組み合わせの数を求めます。
C++ プログラム例（ABC300D）
エラトステネスの篩を使って30万以下の素数表 primes を計算しました（13-44行目）。エラトステネスの篩については、こちらの記事で解説しています。素数表を作りながら、$a$ の上限 limit_a と $b$ の上限 limit_b も計算しておきます（34、38行目）。
$a$ と $b$ の範囲を絞っても、$a^2 \times b \times c^2$ が64ビット整数の範囲を超える可能性がある場合は、$c$ のループを抜けます（61-63行目）。
$a^2 \times b \times c^2 \leqq N$ を満たす場合、変数 result をインクリメントします。最後に result の値を出力します。
以下が、C＋＋プログラムとなります。
#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long int ull;
#define N 300000
// \sqrt{(10^12)/(2*2*3))} = 288675.1345...
bool is_prime_table[N + 1];
vector<int> primes;
int limit_a = -1;
int limit_b = -1;
void generate_primes(int n)
{
is_prime_table[0] = false;
is_prime_table[1] = false;
for (int i = 2; i <= n; ++i) {
is_prime_table[i] = true;
}
for (int i = 2; i * i <= n; ++i) {
if (is_prime_table[i]) {
for (int j = i + i; j <= n; j += i) {
is_prime_table[j] = false;
}
}
}
for (int i = 2; i <= n; ++i) {
if (is_prime_table[i]) {
primes.push_back(i);
if ((limit_a == -1)&&(i >= 252)) {
// (10^12)^(1/5) = 251.1886...
limit_a = primes.size();
}
if ((limit_b == -1)&&(i >= 6300)) {
// ((10^12)/4)^(1/3) = 6299.6052...
limit_b = primes.size();
}
}
}
return;
}
int main()
{
generate_primes(N);
ull n;
cin >> n;
ull a, b, c;
int result = 0;
for (int i = 0; i < limit_a; ++i) {
a = (ull)primes[i];
for (int j = i + 1; j < limit_b; ++j) {
b = (ull)primes[j];
for (int k = j + 1; k < primes.size(); ++k) {
c = (ull)primes[k];
if (c * c * b > n) {
break;
}
if (c * c * b * a * a <= n) {
++result;
} else {
break;
}
}
}
}
cout << result << endl;
return 0;
}ほぼ上記のプログラムをコンテストに提出しました。
コンテスト後、解説を読み、2点改良しました。
- limit_a や limit_b を計算するのは止めて、ループを break するように変更した（48-49、53-54行目）。
- c についての3重ループではなく、二分探索で c の上限を求めて、b のインデックスとの差をとり、result を更新した（56-61行目）。
改良したプログラムが以下になります。プログラムが見やすくなり、計算速度も向上しました。
#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long int ull;
#define N 300000
// \sqrt{(10^12)/(2*2*3))} = 288675.1345...
bool is_prime_table[N + 1];
vector<int> primes;
void generate_primes(int n)
{
is_prime_table[0] = false;
is_prime_table[1] = false;
for (int i = 2; i <= n; ++i) {
is_prime_table[i] = true;
}
for (int i = 2; i * i <= n; ++i) {
if (is_prime_table[i]) {
for (int j = i + i; j <= n; j += i) {
is_prime_table[j] = false;
}
}
}
for (int i = 2; i <= n; ++i) {
if (is_prime_table[i]) {
primes.push_back(i);
}
}
return;
}
int main()
{
generate_primes(N);
ull n;
cin >> n;
ull a, b, c;
int s = primes.size();
int result = 0;
for (int i = 0; i < s; ++i) {
a = (ull)primes[i];
if (a * a * a * a * a > n) {
break;
}
for (int j = i + 1; j < s; ++j) {
b = (ull)primes[j];
if (a * a * b * b * b > n) {
break;
}
c = (ull)sqrt(n / (a * a * b));
if (c <= b) {
break;
}
int upper = upper_bound(primes.begin(), primes.end(), c) - primes.begin();
result += upper - j - 1;
}
}
cout << result << endl;
return 0;
}どちらも、AC（Accepted＝正しいプログラム）と判定されました。
Python プログラム例（ABC300D）
Python 版は、C++ の後者のプログラムを移植しました。
- 平方を求める sqrt を使うために math モジュールをインポートしました（2,37行目）。
- 二分探索を行うため、bisect モジュールをインポートしました（3,40行目）。
以下が、Python のプログラムです。
"""AtCoder Beginner Contest 300 D"""
import math
import bisect
def generate_primes(n):
is_prime_table = [True] * (n + 1)
is_prime_table[0] = False
is_prime_table[1] = False
i = 2
while i * i <= n:
if is_prime_table[i]:
for j in range(i + i, n + 1, i):
is_prime_table[j] = False
i += 1
for i in range(2, n + 1):
if is_prime_table[i]:
primes.append(i)
primes = []
generate_primes(300000)
n = int(input())
s = len(primes)
result = 0
for i in range(s):
a = primes[i]
if a ** 5 > n:
break
for j in range(i + 1, s):
b = primes[j]
if a ** 3 * b ** 2 > n:
break
c = int(math.sqrt(n / (a * a * b)))
if c <= b:
break
upper = bisect.bisect_right(primes, c)
result += upper - j - 1
print(result)こちらも「AC」と判定されました。
最後に
ABC298とABC299がサイバー攻撃により unrated になりました。今回は無事に開催できました。rate が算出されると励みになります。関係者のみなさまに感謝いたします。
引き続き ABC の問題を紹介していきます。
