「40 - 32 / 2 = 4!」

Twitterでフォローしてくださった@rbtnnさんのブログ「http://b-rabbit-alice.blogspot.com/」を拝読していたら、面白いものを見つけてしまいました。

• http://b-rabbit-alice.blogspot.com/2011/07/40-32-2-41000.html

こういうの大好き（笑）。

自分なりにどこまでいけるか、挑戦してみました。
そのせいでいささか長ーいエントリになってます。

ちなみに。ここのコードはGitHubにも登録してあります。

1.基準を確認

まずは、オリジナルのコードの動作の確認。スタイルだけ自分の好みに合わせてありますが、内容は一緒です。これをリファレンスにしてなにをどうしてどうなったかを見ていきます。

#include <iostream>
#include <string>
#include <boost/format.hpp>

void print(double a, double b, double c)
{
    const double m = (a - b) / c;
    const double n = a - b / c;
    std::cout << (boost::format("(%1% - %2%) / %3% = %4%") % a % b % c % m) << std::endl;
    std::cout << (boost::format(" %1% - %2%  / %3% = %4%! = %5%") % a % b % c % m % n) << std::endl;
    std::cout << std::endl;
}

double factorial(double n)
{
    return (n > 0) ? (n * factorial(n - 1)) : 1;
}

int main(int, char* [])
{
    for(double a = 1; a < 1000; ++a)
    {
        for(double b = 1; b < 1000; ++b)
        {
            for(double c = 2; c < 100; ++c)
            {
                if((a > b) && ((a - (b / c)) == factorial((a - b) / c)))
                {
                    print(a, b, c);
                }
            }
        }
    }

    return 0;
}

実行結果。無事おなじ結果が得られました。

(25 - 5) / 5 = 4
 25 - 5  / 5 = 4! = 24

(30 - 18) / 3 = 4
 30 - 18  / 3 = 4! = 24

(40 - 32) / 2 = 4
 40 - 32  / 2 = 4! = 24

(138 - 108) / 6 = 5
 138 - 108  / 6 = 5! = 120

(230 - 220) / 2 = 5
 230 - 220  / 2 = 5! = 120

(728 - 416) / 52 = 6
 728 - 416  / 52 = 6! = 720

(731 - 473) / 43 = 6
 731 - 473  / 43 = 6! = 720

(735 - 525) / 35 = 6
 735 - 525  / 35 = 6! = 720

(748 - 616) / 22 = 6
 748 - 616  / 22 = 6! = 720

(756 - 648) / 18 = 6
 756 - 648  / 18 = 6! = 720

(765 - 675) / 15 = 6
 765 - 675  / 15 = 6! = 720

(816 - 768) / 8 = 6
 816 - 768  / 8 = 6! = 720

(833 - 791) / 7 = 6
 833 - 791  / 7 = 6! = 720

(952 - 928) / 4 = 6
 952 - 928  / 4 = 6! = 720

実行時間を計測してみました。ハードウェアはいつも通りPowerPCのMac、iBook G4 1.33GHz、OSはMac OS X 10.5.8、コンパイルはGCC 4.0.1、オプションは -ansi -Wall で最適化は指定していません。単位は秒でtimeコマンドで計測しています。コンソール出力は/dev/nullに捨てています。

	1回目	2回目	3回目	平均	高速化
1	40.255	40.254	40.242	40.250	1

「高速化」はココが基準なので「1」とします。

2.メモ化(memoization)

階乗の効率化といえば、定番のメモ化(memoization)。

• Memoization - Wikipedia

factorial関数をクラスにして事前に計算した結果は保存して再利用するようにします。
ついでに、浮動小数点演算をなくし、割り切れる場合にだけ式が成り立つかどうかを検証するようにしてみました。
もうひとつついでに、パラメータをコマンドラインから指定できるようにしています。

#include <iostream>
#include <string>
#include <vector>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

class Factorial
{
public:
    Factorial() : factorial_() {}

    int operator [] (int n) const
    {
        if(factorial_.size() <= n)
        {
            factorial_.push_back((n > 0) ? (n * (*this)[n - 1]) : 1);
        }
        return factorial_[n];
    }

private:
    mutable std::vector<int> factorial_;
};

bool isDivisible (int dividend, int divisor)
{
    return (dividend % divisor) == 0;
}

void print(int a, int b, int c)
{
    const int m = (a - b) / c;
    const int n = a - b / c;
    std::cout << (boost::format("(%1% - %2%) / %3% = %4%") % a % b % c % m) << std::endl;
    std::cout << (boost::format(" %1% - %2%  / %3% = %4%! = %5%") % a % b % c % m % n) << std::endl;
    std::cout << std::endl;
}

int main(int argc, char* argv[])
{
    if(argc != 2)
    {
        std::cerr << "usage: " << argv[0] << " <count>\n";
        return 0;
    }

    const int max = boost::lexical_cast<int>(argv[1]);

    Factorial factorial;

    for(int a = 1; a < max; ++a)
    {
        for(int b = 1; b < max; ++b)
        {
            for(int c = 2; c < max; ++c)
            {
                if((a > b) && isDivisible (b, c) && isDivisible (a - b, c) && ((a - b / c) == factorial[(a - b) / c]))
                {
                    print(a, b, c);
                }
            }
        }
    }

    return 0;
}

実行結果は省略…と思ったんですが、ひとつ新しい解を見つけてしまいました。オリジナルではcの探索範囲が1〜100だったのですが、その範囲を広げてみたらcが100より大きい場合の解があるのを発見してしまいました。

(721 - 103) / 103 = 6
 721 - 103  / 103 = 6! = 720

	1回目	2回目	3回目	平均	高速化
2	34.817	34.801	34.803	34.807	1.15

上に書いた通りcの探索範囲を広げた影響か、あんまり早くなってないです。

3.探索する範囲を絞ろう

やっぱり探索する範囲が広過ぎるのがネックになっているようなので、これを絞り込むことを考えます。
まず、a、b、cはいずれも正の整数で、(a - b) / cの結果も正の整数とするとb < aでないとなりません。すでにふるい分けの条件文の中にこの式がありますが、ループ自体から取り除くことを考えます。また、(a - b) / cが正の整数ということは、分子より分母が小さくないといけないのでc < a - bのはずです。

これらを元に書き直したのが次のコード。

#include <iostream>
#include <string>
#include <vector>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

class Factorial
{
public:
    Factorial() : factorial_() {}

    int operator [] (int n) const
    {
        if(factorial_.size() <= n)
        {
            factorial_.push_back((n > 0) ? (n * (*this)[n - 1]) : 1);
        }
        return factorial_[n];
    }

private:
    mutable std::vector<int> factorial_;
};

bool isDivisible (int dividend, int divisor)
{
    return (dividend % divisor) == 0;
}

void print(int a, int b, int c)
{
    const int m = (a - b) / c;
    const int n = a - b / c;
    std::cout << (boost::format("(%1% - %2%) / %3% = %4%") % a % b % c % m) << std::endl;
    std::cout << (boost::format(" %1% - %2%  / %3% = %4%! = %5%") % a % b % c % m % n) << std::endl;
    std::cout << std::endl;
}

int main(int argc, char* argv[])
{
    if(argc != 2)
    {
        std::cerr << "usage: " << argv[0] << " <count>\n";
        return 0;
    }

    const int max = boost::lexical_cast<int>(argv[1]);

    Factorial factorial;

    for(int a = 1; a < max; ++a)
    {
        for(int b = 1; b < a; ++b)
        {
            const int c_max = a - b;
            for(int c = 2; c < c_max; ++c)
            {
                if(isDivisible (b, c) && isDivisible (a - b, c) && ((a - b / c) == factorial[(a - b) / c]))
                {
                    print(a, b, c);
                }
            }
        }
    }

    return 0;
}

	1回目	2回目	3回目	平均	高速化
3	8.570	8.661	8.667	8.633	4.66

そこそこ早くなってきました。

4.もっと絞れないか？

cの範囲ですが、a - (b / c)も正の整数になることするとc < bでもあるはずです。ということはa - bとbの小さい方よりもcは小さくなるはず。

というわけで書き換えたのが以下のコード。変更が少ないので部分だけ。

これを、

            const int c_max = a - b;
            for(int c = 2; c < c_max; ++c)

こう変更。条件式の演算子が<=になっていることに注意。

            const int c_max = std::min(a - b, b);
            for(int c = 2; c <= c_max; ++c)

また、min関数テンプレートを使うためヘッダファイルを追加で読み込み。

#include <algorithm>

	1回目	2回目	3回目	平均	高速化
4	4.502	4.518	4.545	4.188	9.61

まぁまぁ、早くなってきました。

5.探索方法を見直してみる

とはいえ。ここまでのやり方では三重ループなのはかわらないので、せいぜい定数倍でしか早くなりません。扱う数が大きくなるにつれ手に負えなくなるのはかわりません。
根本的に探索方法を変える必要がありそうです。

a - b / c = n!という式を見たとき、目につくのはnの値の小ささ。n!は急激に値が大きくなるので、n自体はそれほど大きくなりません。1000まで探索してもせいぜい6どまり。ならばということで、nをループの基準においてみます。

一番外側、aは任意の正の整数とします。

a - b / c = n!ですから、a = n! + b / cで、b / cも正の整数ですから、a > n!が成り立ちます。a > n!となるnのうち最大のものをnmaxとすれば、1≦n≦nmaxです。

またa - b / c = n!から、b / c = a - n!になります。
cが1のときのbの値をb1とすると、b / c = b1 / 1 = b1 = a - n!となります。実際にはcは正の整数なのでbはb1の倍数、つまりa - n!の倍数になります。

さらに(a - b) / c = nで、nは正の整数ですから、(a - b) / cも分母より分子が大きくなければなりません。ゆえにc < a - bです。

結果として、

aのループ

nのループ

bの（b1の倍数の）ループ

となります。内側のふたつのループがかなり小さくなるはずなので、効率化が期待できます。

実装。

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

class Factorial
{
public:
    Factorial() : factorial_() {}

    int operator [] (int n) const
    {
        if(factorial_.size() <= n)
        {
            factorial_.push_back((n > 0) ? (n * (*this)[n - 1]) : 1);
        }
        return factorial_[n];
    }

    int getMaxLessOrEqual(int a) const
    {
        int result = 0;
        while((*this)[result] < a)
        {
            ++result;
        }
        return --result;
    }

private:
    mutable std::vector<int> factorial_;
};

bool isDivisible (int dividend, int divisor)
{
    return (dividend % divisor) == 0;
}

void print(int a, int b, int c)
{
    const int m = (a - b) / c;
    const int n = a - b / c;
    std::cout << (boost::format("(%1% - %2%) / %3% = %4%") % a % b % c % m) << std::endl;
    std::cout << (boost::format(" %1% - %2%  / %3% = %4%! = %5%") % a % b % c % m % n) << std::endl;
    std::cout << std::endl;
}

int main(int argc, char* argv[])
{
    if(argc != 2)
    {
        std::cerr << "usage: " << argv[0] << " <count>\n";
        return 0;
    }

    const int max = boost::lexical_cast<int>(argv[1]);

    Factorial factorial;

    for(int a = 20; a < max; ++a)
    {
        int n_max = factorial.getMaxLessOrEqual(a);

        for(int n = 1; n <= n_max; ++n)
        {
            int b1 = a - factorial[n];
            for(int c = 2; (b1 * c) < a; ++c)
            {
                const int b = b1 * c;
                if(isDivisible(a - b, c) && (((a - b) / c) == n))
                {
                    print(a, b, c);
                }
            }
        }
    }

    return 0;
}

	1回目	2回目	3回目	平均	高速化
5	0.022	0.022	0.022	0.022	1830

みごとに高速化できました。

6.集計

	1回目	2回目	3回目	平均	高速化
1	40.255	40.254	40.242	40.250	1
2	34.817	34.801	34.803	34.807	1.15
3	8.570	8.661	8.667	8.633	4.66
4	4.502	4.518	4.545	4.188	9.61
5	0.022	0.022	0.022	0.022	1830

不満があるとすれば、数学的に解いてみたかったというところ。
演算量を見積もってそれにそってコードを改善していくというのは、数値計算の考えからいえば正攻法だと思うんですが、できれば規則性とか見つけてみたいなー、と。
いまは、そこまで頭が働いてませんです、はい。

7.おまけ：C++でできるならRubyでもできるよね…？

ここまで高速にできるならインタプリタでも充分実用になるよね、ということでRubyに翻訳。直訳なのでRubyっぽくないかもしれません。

class Factorial
  def initialize
    @factorial = []
  end

  def get(n)
    @factorial << ((n > 0) ? (n * get(n - 1)) : 1) if @factorial.size() <= n
    return @factorial[n]
  end

  def lower_bound(a)
    result = 0
    while get(result) < a
      result += 1
    end
    result - 1
  end
end

def is_divisible(dividend, divisor)
  ((dividend / divisor) * divisor) == dividend
end

def print_abc(a, b, c)
  m = (a - b) / c
  n = a - b / c
  puts "(#{a} - #{b}) / #{c} = #{m}"
  puts " #{a} - #{b}  / #{c} = #{m}! = #{n}"
  puts
end

Max = ARGV.shift.to_i

factorial = Factorial.new

(1...Max).each do |a|
  m = factorial.lower_bound a
  (1..m).each do |n|
    b1 = a - factorial.get(n)
    c  = 2
    while (b1 * c) < a
      b = b1 * c
      if is_divisible(a - b, c) && (((a - b) / c) == n)
        print_abc(a, b, c)
      end
      c += 1
    end
  end
end

実行結果は省略。

8.おまけ、その２：ぢゃ、Haskellでは？

行き詰まりました（爆）。すっかり腕が鈍っていて、書いているうちに、わけわかんなくなりました。Haskell篇は後日ということで…。