読者です 読者をやめる 読者になる 読者になる

初代Masteries

きっとモヒカンにもなれないお前たちに告げる!!!

オイラー王になれませんでした!

perl euler

弊社(GaiaX)にはプログラミング部と呼ばれる部活動があります.
始業前に「朝練」と呼んでいる自主勉強会を開催したり, Project Eulerの問題を解く会を開いたりして, お互いの技能を高め合ったりしています*1.

この部活動ですが, 今年度に入ってから, 「Peroject Eulerの問題を解く会」を更に発展(?)させて, 「Euler王決定戦」を開催しています.
こんな感じで紹介してもらったりしているようですが, 要するにProject Eulerの問題を題材にした, 社内での競技プログラミング大会... みたいな取り組みです.

先月開催された第1回のEuler王選手権では, @kazuphさんが優勝して, Euler王の証のトロフィーを授与されていました.
自分は現在大阪住まいではありますが, プログラミング部部長の@tyokomineさんご好意で, 大学から遠隔参加(?)させて頂いて, @kazuphさんからトロフィーを奪還すべく頑張った... のですが, さすがに2時間では無理でした.



...が, 8時間悩みに悩んだ末, なんとか答えを導くことができたので貼っておきます.

Problem 108

use strict;
use warnings;
use Math::Prime::XS qw/primes/;

my @primes = primes(17);
my @answer = (@primes, @primes);

my $min_divisor = 2001;
my $max = mul(map {$_ ** 2} @primes);

search_factors(2, 2);
print sqrt(mul(@answer)) . "\n";

sub search_factors {
    my @factors = @_;

    my $divisor = divisor(@factors);
    if (mul(@factors) >= $max) {
        return;
    }
    if ($divisor >= $min_divisor && $divisor <= divisor(@answer)) {
        @answer = @factors;
        $max = mul(@factors);
        return;
    }

    my $next_prime = next_prime($factors[-1]);
    my $before_prime = before_prime($factors[-1]);

    search_factors(@factors, $next_prime, $next_prime) if $next_prime != 0;
    if ($factors[-1] == 2
        || $before_prime != 0
        && count(\@factors, $before_prime) > count(\@factors, $factors[-1])
    ) {
        my $bottom = $factors[-1];
        search_factors(@factors, $bottom, $bottom);
    }
}

sub mul {
    my @factors = @_;

    my $mul = 1;
    $mul *= $_ for @factors;
    return $mul;
}

sub divisor {
    my @factors = @_;

    my $divisors = {};

    for my $factor (@factors) {
        $divisors->{$factor}++;
    }

    my $divisor = 1;
    $divisor *= $_ for map { $_ + 1 } values %{$divisors};

    return $divisor;
}


sub count {
    my ($factors, $num) = @_;

    my $count;
    for my $factor (@{$factors}) {
        $count++ if $factor == $num;
    }
    return $count;
}

sub before_prime {
    my $prime = shift;

    return 0 if $prime == 2;

    my $index = 0;
    for my $p (@primes) {
        return $primes[$index-1] if defined $primes[$index-1] && $prime == $p;
        $index++;
    }
    return 0;
}

sub next_prime {
    my $prime = shift;

    my $index = 0;
    for my $p (@primes) {
        return $primes[$index+1] if defined $primes[$index+1] && $prime == $p;
        $index++;
    }
    return 0;
}

1 / x + 1 / y = 1 / n という式を変形すると, (y + x) / xy = 1 / n となります.
両辺に n/xy を掛けて, n * (y + x) = xy. 変形して, xy - nx - ny = 0.
ここで両辺にn^2を足すと, xy - nx - ny + n^2 = n^2 となり, (x - n)(y - n) = n ^ 2 となります.

こうすることで, この問題は, n ^ 2 の約数の個数を求める問題ととらえることができます.
但し, (x, y) = (y, x)になる場合と, x = yになる場合があるので, 約数の個数は2001個以上でなければなりません.
逆に言えば, (n ^ 2 の約数の個数 / 2) + 1 が, nにおける約数の数となります.

例えば, 問題の解説には, n = 4 には3つの異なる解がある, とありました.
n ^ 2 = 4 ^ 2 = 16 で, 16には約数が5個あります. (5 / 2) + 1 = 3, というわけで正しいですね.

次に, n ^ 2 の約数ですが, 約数は素因数分解の結果から求められることを利用して解いていきます.
例えば, 先ほどの n = 4, n ^ 2 = 16 の場合, 16を素因数分解すると 2 * 2 * 2 * 2 = 2 ^ 4です.
約数は, Σ(素因数分解の結果の指数 + 1) で求めることができるので, 2 ^ 4 の場合, 約数は5個になります.
また, 360の場合, 素因数分解をすると 2 * 2 * 2 * 3 * 3 * 5 = 2 ^ 3 * 3 ^ 2 * 5 ^ 1 なので, (3 + 1) * (2 + 1) * (1 + 1) = 24個の約数を持ちます.

また, n ^ 2 を素因数分解した結果では, 素数の指数はかならず偶数でなければなりません.
そうでなければ(指数が奇数の素数があれば), nが整数にならないからです.
スクリプト中で, n ^ 2 を求める為に素因数分解の候補を大量に生成していますが, この特徴を使うことで, 素数を格納する配列に対して, 一度に2個の要素を挿入することができます.

更に, n ^ 2 を素因数分解した結果は, 2, 3, 5, 7... と, 2から始まる連続した素数で構成され, その指数は数が大きくなるにつれ, 減っていきます.
例えば, 2 ^ 2 * 3 ^ 2 と, 2 ^ 2 * 5 ^ 2 は, どちらも約数は6個ですが, 2 ^ 2 * 5 ^ 2 の方が大きな数となります.
今回求めたいのは, 1 / x + 1 / y = 1 / n について, 異なる解の数が1000を越える, 最小の n を求めなければならないので, 2 ^ 2 * 3 ^ 2 でなければなりません.
2 ^ 2 * 3 ^ 1 と 2 ^ 1 * 3 ^ 2 についても同様で, 約数はどちらも同じですが, 2 ^ 1 * 3 ^ 2 の方が大きな数になるので, 2 ^ 2 * 3 ^ 1 のようにならなければなりません.

答えとなり得る n ^ 2 の最大値とその際の約数ですが, 先ほどの規則を適用すると, 3 ^ 6 = 729, 3 ^ 7 = 2187 より約数が2187個以下, n ^ 2 <= 2 ^ 2 * 3 ^ 2 * 5 ^ 2 * 7 ^ 2 * 11 ^ 2 * 14 ^ 2 = 260620460100 となります.
また, これより必要となる素数は17以下でよいこともわかります.

これらの特徴を用いて, 約数が2001個以上の n ^ 2 について, n ^ 2 の約数の数が最も小さくなる n が答えとなります.
なお, n ^ 2 の約数の数が同数になる場合がありますが, その場合は n ^ 2 そのものの値が一番小さいものを選べばよいです.

まとめ

クッソ汚いコードですが, 出てきた答えをProject Eulerのページで登録すると「正解!」という感じだったので, もう... ゴールしてもいいよね...

とりあえず, 疲れたので, 寝ます. また起きてから加筆修正するかもしれません...

*1:自分も何度かお邪魔させて頂いたことがあります.