AtCoder

ABC300 D問題(AABCC)を解く

AtCoder_ABC300_D

AtCoder が提供しているABC(AtCoder Beginner Contest)300 のD問題をC++とPythonで解いてみました。ABC300は、2023年4月29日21:00に実施されました。

AtCoder の紹介はこちらに、プログラミングの方針はこちらに記事があります。

D問題 AABCC(Difficulty : 908)

問題はリンク先をご覧ください。

ABC300 D問題 AABCC

素数表を計算して、条件を満たす素数の組み合わせをカウントします。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 に挑戦して、解説記事を書いていくつもりです。

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

CAPTCHA