CodeIQ MAGAZINECodeIQ MAGAZINE

数学の問題をプログラミングで解こう!「スクエア・カルテット」問題解説

2016.01.18 Category:CodeIQ問題解説・リーダーボード Tag: ,

  • 10
  • このエントリーをはてなブックマークに追加

「数学の問題をプログラミングで解こう!」シリーズ。

平方数の和が等しくなるような整数の組を見つける問題です。
単純に考えると膨大な計算量を要しますが、うまく数式を変形するとぐっと計算量を抑えられます。
ぜひ本記事で解法をチェックして下さい!
by @riverplus Kawazoe

問題のおさらい

2つの自然数の組 (a, b) が与えられたとき、自然数 x, y に関する次の方程式を考えます。

    x2a2y2 + b2 … (※)

例えば、(a, b) = (3, 10) のとき、方程式(※)の解は (x, y) = (10, 3), (46, 45) の 2 組です。

自然数の組 (a, b) に対し、方程式(※)の全ての解の x + y の和を F(a, b) と定義します。

例えば F(3, 10) = 10 + 3 + 46 + 45 = 104 です。

同様に、F(10, 50) = 3500、F(20, 100) = 15022 となることが確かめられます。

標準入力から、半角空白区切りで 2 つの自然数 a, b(1 ≦ a < b ≦ 105)が与えられます。
標準出力に F(a, b) の値を出力するプログラムを書いてください。

方針

与えられた a, b に対し、等式を成立させる整数解を求める問題です。

単純な解法としては、y = √(x2 + a2b2) と変形するやり方が考えられます。x の値を一つ選び、x2 + a2b2 の値が平方数になるかどうかをチェックするというやり方です。

ただこの方法はあまり高速ではありません。x の値の取り得る範囲が広すぎるためです。yx – 1 でなくてはなりませんから、上の y の式に代入すると、√(x2 + a2b2) ≦ x – 1 となります。両辺を 2 乗して整理すると、x ≦ (b2a2 + 1) / 2 という条件が得られます。a < b ≦ 105 という条件では、調べるべき範囲が膨大になり、多くの計算量を要してしまいます。

そこで、本問のポイントは、次の式変形を行うことです。

    (xy) (xy) = b2a2

変数 x, y を左辺に集めて因数分解しました。このように、整数に関する問題では、与式を因数分解した形で表すと道筋が見えてくることが多いです。さて、この式をじっと見てみると、x + yxy の積が b2a2 に等しいと分かります。つまり、b2a2 の約数が分かれば、それがx + y または xy の候補になるということです。

b2a2(以下、K とします)の約数を求めるコードを書きましょう。割る数 d を 1 から始めて 1 ずつ増やしながら、K を割り切るかどうか確かめればOKです。なお、dK までループさせる必要はありません。dK の約数なら K / dK の約数であるという事実を利用しましょう。d を √K までループさせれば十分です。また xy がちゃんと整数になるかどうかも忘れずチェックしましょう。コードは次のようになります(C++)。

#include <iostream>
#include <vector>
using namespace std;

long long Sqrt(long long n) {
    long long s, t;
    if (n <= 0) return 0;
    s = 1;  t = n;
    while (s < t) { s <<= 1;  t >>= 1; }
    do {
        t = s;
        s = (n / s + s) >> 1;
    } while (s < t);
    return t;
}

// nの約数の一覧を返す関数
vector<long long> Divisors(long long n) {
    long long s = Sqrt(n);
    vector<long long> v;
    for (long long i = 1; i <= s; i++) {
        if (n % i != 0) continue;
        v.push_back(i);
        v.push_back(n / i);
    }
    return v;
}

int main()
{
    long long a, b;
    cin >> a >> b;
    long long answer = 0;
    long long z = b * b - a * a;
    vector<long long> v = Divisors(z);
    for (int i = 0; i < v.size(); i++) {
        long long p = v[i], q = z / v[i];
        if ((p + q) % 2 != 0 || p <= q) continue;
        long long x = (p + q) / 2, y = (p - q) / 2;
        answer += x + y;
    }
    cout << answer << endl;
}

本問の条件では上記のコードでクリアできます。

さらなる高速化

本問ではさらなる高速化が可能です。ここでも因数分解が威力を発揮します。右辺を次のように変形しましょう。

    (xy) (xy) = (ba) (ba)

一つ前のやり方では、b2 - a2 をさまざまな約数の候補 d で除算していましたが、ここではその必要はありません。右辺は x + yx - y に分解できるわけですから、まずは x + y の約数、x - y の約数をそれぞれ求めましょう。そして、それぞれの約数どうしの組み合わせについて、積を求めれば、それが b2 - a2 の約数になると分かります。

一つ前のやり方では、約数の探索のために √(b2 - a2) 回ほどの除算を要していました。今回の方法では、√(b + a) + √(b - a) 回ほどの除算で約数が求められます。計算量を一気に抑えることが可能になります。

以上のコードを以下に記します。ab が 1010 位の大きさであっても、十分短い実行時間で答えが得られます(あまり大きいと最終的な答えが 64 ビット整数の範囲を超えるため、多倍長整数をサポートする python で書いてます)。また同じ約数がダブってしまわないよう、set を用いています。

def sqrt(n):
    if n <= 0: return 0
    s, t = 1, n
    while s < t:
        s, t = s << 1, t >> 1
    while True:
        s, t, = (n / s + s) >> 1, s
        if s >= t: break
    return t

def divisors(n):
    s = sqrt(n)
    v = set()
    for i in range(1, s+1):
        if n % i != 0: continue
        v.add(i)
        v.add(n / i)
    return v

a, b = map(int, raw_input().split())
v = divisors(b + a)
u = divisors(b - a)
w = set()
answer = 0
for p in v:
    for q in u:
        w.add(p * q)
for p in w:
    q = (b * b - a * a) / p
    if (p + q) % 2 != 0 or p <= q: continue
    x, y = (p + q) / 2, (p - q) / 2
    answer += x + y

print answer

みんなのコード

本問の挑戦者で、コードを公開して頂いた方のコードを Togetter でまとめました。(随時追加していきます。)こちらもぜひチェックして下さい。

CodeIQ「スクエア・カルテット」問題 みんなのコード

CodeIQ運営事務局より

Kawazoeさん、ありがとうございました!
Kawazoeさんの次の問題にご期待ください!

  • 10
  • このエントリーをはてなブックマークに追加

■この記事を書いた人

avatar

@riverplus Kawazoe

数学・パズル好きのソフトウェアエンジニアです。某社で研究開発をしています。Project Euler の問題開発チームメンバー(最近はあまり活動せず)。twitter: @riverplus Web Site: riverplus.net

■関連記事

【謎解きプログラム】この処理は?【コードを読もう】解答と解説... 【謎解きプログラム】この処理は?【コードを読もう】 本問題は、表題のテーマで、プログラムにちなんだ謎を解くというものでした。 それでは以下、各問題とその解答を見ていきましょう。 問題のオープニング ある日、出社すると、あなたのPCのログイン画面に、謎の挑戦状が表示されていた。 「24時間以...
数学の問題をプログラミングで解こう!「ループ・トラッキング」問題解説... 問題のおさらい 自然数 n に対し、関数 Fn(x) を次のように定義します(floor():床関数)。 例えば n=10, x=1 のとき、F10(1) = floor(4×1×9÷10) = 3 です。 さて、整数 k(0 ≦ k ≦ n)に対して、関数 Fn による変換を繰り返し行い...
【謎解きプログラム】どんな結果になる?【アロー関数】解答と解説... 【謎解きプログラム】どんな結果になる?【アロー関数】 本問題は、表題のテーマで、プログラムにちなんだ謎を解くというものでした。 それでは以下、各問題とその解答を見ていきましょう。 問題のオープニング ある日、出社すると、あなたのPCのログイン画面に、謎の挑戦状が表示されていた。 「24時間...
数学の問題をプログラミングで解こう!「カウント・スリー」問題解説... 問題のおさらい 自然数を 1 から順に書き並べていきます。 このとき、3 の数字が現れる回数を数えます。 自然数 n に対し、ちょうど n 個目の 3 の数字が現れたときに書いている数を F(n) と定義します。 例えば F(10)=35 です。 下の通り、10 個目の 3 は、35 を書いて...
【息抜き】カードを上手く並べよう【言語不問】解答と解説... 【息抜き】カードを上手く並べよう【言語不問】 本問題は、表題のテーマで、簡単なプログラムを書くものです。 それでは以下、問題とその解答を見ていきましょう。 問題 あなたは、11から99までの、89枚のカードを持っています。問題では、横幅と高さの整数が与えられます。この横幅と高さで作られるマス...
【コードミステリ】数字に隠されたメッセージ【言語不問】解答と解説... 【コードミステリ】数字に隠されたメッセージ【言語不問】 本問題は、表題のテーマで、簡単なプログラムを書くものです。 喜屋武ちあきさんによるCodeIQ MAGAZINEでのブックレビューに合わせて、『顔貌売人』(文藝春秋)とのコラボ問題として出題されたものです。 それでは以下、問題とその解...

今週のPickUPレポート

新着記事

週間ランキング

CodeIQとは

CodeIQ(コードアイキュー)とは、自分の実力を知りたいITエンジニア向けの、実務スキル評価サービスです。

CodeIQご利用にあたって
関連サイト
codeiq

リクルートグループサイトへ