Ruby - 多桁乗算（Toom-Cook 法 (3-way)）！

前回は、C++ による「多桁乗算（Toom-Cook 法 (3-way)）」のアルゴリズムを紹介しました。

• C++ - 多桁乗算（Toom-Cook 法 (3-way)）！

今日は、同じアルゴリズムを Ruby で実現してみました。
Ruby では桁数（整数型の範囲）をあまり気にしなくても、メモリの許される限り計算できますが、それでも多桁同士の乗算では不都合になることもあるでしょう。 アルゴリズムについては、上記リンクの過去記事等を参照してください。

実際、ほとんど同じです。

以下、Ruby によるサンプルスクリプトです。

0. 前提条件

• Linux Mint 14 Nadia (64bit) での作業を想定。
• Ruby 2.0.0-p0 を使用。

1. Ruby スクリプト作成

スクリプト中の随所にあるテスト用コメントを解除すると、１ループあたりの乗算回数をカウントしたり、100 ループした場合の処理時間を計測するようになっている。
また、計算可能な最大桁数の設定を 3 のべき乗桁にしている。

File: multiply_toom_cook_3.rb

#! /usr/local/bin/ruby
#*********************************************
# 多倍長乗算 ( by Toom-Cook 法 (3-way) )
#   - 多倍長 * 多倍長
#   - 最下位の桁を配列の先頭とする考え方
#*********************************************
#
class MultiplyToomCook3
  D_MAX = 729  # 計算可能な最大桁数 ( 3 のべき乗 )
  D     = 729  # 実際に計算する桁数 ( D_MAX 以下 )

  def initialize
    # ====================================== #
    # テストなので、被乗数・乗数は乱数を使用 #
    # ====================================== #
    @a, @b = Array.new, Array.new
    rnd_a, rnd_b = Random.new(0), Random.new(1)
    # 被乗数・乗数設定
    0.upto(D - 1) do |i|
      @a << rnd_a.rand(10)
      @b << rnd_b.rand(10)
    end
  end

  # 計算 ( Toom-Cook 法 )
  def calc_toom_cook
    # 配列初期設定 ( コンストラクタで生成した配列を使用 )
    a, b = @a, @b  # 被乗数配列、乗数配列
    z = Array.new  # 計算結果用配列

    begin
      # 最大桁に満たない部分は 0 を設定
      (D_MAX - a.size).times { |i| a << 0 }
      (D_MAX - b.size).times { |i| b << 0 }
      # ====[ TEST ]===>
      # t1 = Time.now      # 計算開始時刻
      # 100.times do |l|  # 100 回 LOOP
      # @cnt_mul = 0       # 乗算回数リセット
      # <===[ TEST ]====
      # 乗算 ( Toom-Cook 法 (3-way) )
      z = multiply_toom_cook_3(a, b)
      # ====[ TEST ]===>
      # end
      # t2 = Time.now      # 計算終了時刻
      # @tt = t2 - t1      # 計算時間
      # <===[ TEST ]====
      # 繰り上がり処理
      z = do_carry(z)
      # 結果出力
      display(a, b, z)
    rescue => e
      raise
    end
  end

  private

  # 乗算 ( 標準(筆算)法 )
  def multiply_normal(a, b)
    # 配列サイズ取得
    a_len, b_len = a.size, b.size
    # 計算結果初期化
    z = Array.new(a_len + b_len, 0)
    # 各配列を各桁とみなして乗算
    b_len.times do |j|
      a_len.times do |i|
        z[j + i] += a[i] * b[j]
        # ====[ TEST ]===>
        # @cnt_mul += 1  # 乗算カウント
        # <===[ TEST ]====
      end
    end
    return z
  rescue => e
    raise
  end

  # 乗算 ( Toom-Cook 法 (3-way) )
  #   結果用配列は以下のように配置し、
  #     +----+----+----+----+----+----+----+----+----+----+
  #     |   c0    |   c2    |   c4    |   c1    |   c3    |
  #     +----+----+----+----+----+----+----+----+----+----+
  #   最後に、c1, c3 を所定の位置に加算する。
  #     +----+----+----+----+----+----+
  #     |   c0    |   c2    |   c4    |
  #     +----+----+----+----+----+----+
  #          +----+----+----+----+
  #          |   c1    |   c3    |
  #          +----+----+----+----+
  def multiply_toom_cook_3(a, b)
    # 各種宣言
    a_m1, a_m2      = Array.new, Array.new
    a_0, a_1, a_inf = Array.new, Array.new, Array.new
    b_m1, b_m2      = Array.new, Array.new
    b_0, b_1, b_inf = Array.new, Array.new, Array.new
    c_m1, c_m2      = Array.new, Array.new
    c_0, c_1, c_inf = Array.new, Array.new, Array.new
    c0, c1, c2      = Array.new, Array.new, Array.new
    c3, c4          = Array.new, Array.new

    begin
      # 配列サイズ取得
      t_len = a.size
      # ９桁（配列９個）になった場合は標準乗算
      return multiply_normal(a, b) if t_len <= 9
      # 配列分割
      a2 = a[(t_len * 2 / 3)..-1]
      a1 = a[(t_len / 3)..(t_len * 2 / 3 - 1)]
      a0 = a[0..(t_len / 3 - 1)]
      b2 = b[(t_len * 2 / 3)..-1]
      b1 = b[(t_len / 3)..(t_len * 2 / 3 - 1)]
      b0 = b[0..(t_len / 3 - 1)]
      # a(-2) = 4 * a2 - 2 * a1 + a0, b(1) = 4 * b2 - 2 * b1 + b0 (by シフト演算)
      (t_len / 3).times do |i|
        a_m2 << (a2[i] << 2) - (a1[i] << 1) + a0[i]
        b_m2 << (b2[i] << 2) - (b1[i] << 1) + b0[i]
      end
      # a(-1) = a2 - a1 + a0, b(1) = b2 - b1 + b0
      (t_len / 3).times do |i|
        a_m1 << a2[i] - a1[i] + a0[i]
        b_m1 << b2[i] - b1[i] + b0[i]
      end
      # a(0) = a0, b(0) = b0
      a_0, b_0 = a0, b0
      # a(1) = a2 + a1 + a0, b(1) = b2 + b1 + b0
      (t_len / 3).times do |i|
        a_1[i] = a2[i] + a1[i] + a0[i]
        b_1[i] = b2[i] + b1[i] + b0[i]
      end
      # a(inf) = a2, b(inf) = b2
      a_inf, b_inf= a2, b2
      # c(-2) = a(-2) * b(-2)
      c_m2 = multiply_toom_cook_3(a_m2, b_m2)
      # c(-1) = a(-1) * b(-1)
      c_m1 = multiply_toom_cook_3(a_m1, b_m1)
      # c(0) = a(0) * b(0)
      c_0 = multiply_toom_cook_3(a_0, b_0)
      # c(1) = a(1) * b(1)
      c_1 = multiply_toom_cook_3(a_1, b_1)
      # c(inf) = a(inf) * b(inf)
      c_inf = multiply_toom_cook_3(a_inf, b_inf)
      # c4 = 6 * c(inf) / 6
      c4 = c_inf
      # c3 = -c(-2) + 3 * c(-1) - 3 * c(0) + c(1) + 12 * c(inf) / 6
      ((t_len / 3) * 2).times do |i|
        c3[i]  = -c_m2[i]
        c3[i] += (c_m1[i] << 1) + c_m1[i]
        c3[i] -= (c_0[i] << 1) + c_0[i]
        c3[i] += c_1[i]
        c3[i] += (c_inf[i] << 3) + (c_inf[i] << 2)
        c3[i] /= 6
      end
      # c2 = 3 * c(-1) - 6 * c(0) + 3 * c(1) - 6 * c(inf) / 6
      ((t_len / 3) * 2).times do |i|
        c2[i]  = (c_m1[i] << 1) + c_m1[i]
        c2[i] -= (c_0[i] << 2) + (c_0[i] << 1)
        c2[i] += (c_1[i] << 1) + c_1[i]
        c2[i] -= (c_inf[i] << 2) + (c_inf[i] << 1)
        c2[i] /= 6
      end
      # c1 = c(-2) - 6 * c(-1) + 3 * c(0) + 2 * c(1) - 12 * c(inf) / 6
      ((t_len / 3) * 2).times do |i|
        c1[i]  = c_m2[i]
        c1[i] -= (c_m1[i] << 2) + (c_m1[i] << 1)
        c1[i] += (c_0[i] << 1) + c_0[i]
        c1[i] += (c_1[i] << 1)
        c1[i] -= (c_inf[i] << 3) + (c_inf[i] << 2)
        c1[i] /= 6
      end
      # c0 = 6 * c(0) / 6
      c0 = c_0
      #  z = c4 * x^4 + c3 * x^3 + c2 * x^2 + c1 * x + c0
      z = c0 + c2 + c4
      ((t_len / 3) * 2).times { |i| z[i + t_len / 3] += c1[i] }
      ((t_len / 3) * 2).times { |i| z[i + t_len    ] += c3[i] }
      return z
    rescue => e
      raise
    end
  end

  # 繰り上がり処理
  def do_carry(a)
    cr = 0  # 繰り上がり

    begin
      a.size.times do |i|
        a[i] += cr
        cr = a[i] / 10
        a[i] -= cr * 10
      end
      # オーバーフロー時
      puts "[ OVERFLOW!! ] #{cr}" unless cr == 0
      return a
    rescue => e
      raise
    end
  end

  # 結果出力
  def display(a, b, z)
    # 上位桁の不要な 0 を削除するために、配列サイズを取得
    len_a, len_b, len_z = D_MAX, D_MAX, D_MAX * 2
    while a[len_a - 1] == 0 do len_a -= 1 if a[len_a - 1] == 0 end
    while b[len_b - 1] == 0 do len_b -= 1 if b[len_b - 1] == 0 end
    while z[len_z - 1] == 0 do len_z -= 1 if z[len_z - 1] == 0 end
    # a 値
    puts "a ="
    (len_a - 1).downto(0) do |i|
      print a[i]
      print " " if (len_a - i) % 10 == 0
      puts if (len_a - i) % 50 == 0
    end
    puts
    # b 値
    puts "b ="
    (len_b - 1).downto(0) do |i|
      print b[i]
      print " " if (len_b - i) % 10 == 0
      puts if (len_b - i) % 50 == 0
    end
    puts
    # z 値
    puts "z ="
    (len_z - 1).downto(0) do |i|
      print z[i]
      print " " if (len_z - i) % 10 == 0
      puts if (len_z - i) % 50 == 0
    end
    puts; puts
    # ====[ TEST ]====
    # puts "Counts of multiply / 1 loop = #{@cnt_mul}"     # 乗算回数
    # puts "Total time of all loops     = #{@tt} seconds"  # 処理時間
    # ====[ TEST ]====
  rescue => e
    raise
  end
end

if __FILE__ == $0
  begin
    # 計算クラスインスタンス化
    obj = MultiplyToomCook3.new
    # 乗算 ( Toom-Cook 法 )
    obj.calc_toom_cook
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
  end
end

• Gist - Ruby script to multiply big-digit values with Toom-Cook 3-way method.

2. 実行

まず、実行権限を付与。

$ chmod +x multiply_toom_cook_3.rb

そして、実行。

$ ./multiply_toom_cook_3.rb
a =
4345367000 0015558046 5992851545 8392322526 6661657161 
2940646222 3693630633 2878928428 2931075189 4459942962 
0472047143 1713204503 7429609429 7320661140 5711218714 
8538514523 9923578665 1579065598 3315649556 2636468353 
2850188038 2090122465 2291725695 6708063266 7932721494 
7942319492 9163030723 9638085573 9503746353 7468715593 
8041892365 3487011431 2123396480 0444130541 4830691976 
6375317951 4967144872 8029471638 9868078876 6991300561 
7275325293 3147458058 7336424609 4074817783 0788334840 
6091902995 9132761677 5255524121 1237390448 8301180552 
0485182215 7793237999 7335435793 2856480756 3632388619 
6682988756 9408135690 9580896340 3482838802 3288833606 
0032864192 1432332927 2095702302 4210503951 0557348443 
4464495302 7639971189 4148655400 2723740991 0733318320 
5303498951 8776188674 253973305
b =
7501324065 2738206911 5039490164 4021680779 1593935265 
4517107481 1563095300 4247865767 6579175671 6550946828 
8849310928 2032373166 1809936890 4395890008 8505100323 
3262151146 5641179972 5027654271 6778904491 2603452658 
0104477210 0681082763 6231568633 0763932142 2664695603 
6038377177 6784204672 5050316333 9319575588 5651066222 
3741358300 5240023426 8409350815 7894695159 1170270888 
7057292645 2461606103 7406689505 8225079836 1284499328 
0452193943 0049246899 9592674603 6763426074 0219495812 
5062777866 5259459760 8835743093 0760692006 9696613714 
0412820715 7653833044 2353763341 8309248038 2706034215 
3024307811 8585455729 8864419390 7598704512 1560641560 
7726629062 7213289104 8971702094 7744875265 4680311778 
0379772086 3954773968 9772940293 0418439156 3789388101 
9679960719 7742425429 671005985
z =
3259600604 9558376990 1785049682 6721702532 5044456662 
7821751560 5036432472 4374806508 0089563122 9955094382 
1333781499 2055893283 6740957683 1430717530 3027924720 
6720795606 4453432868 2206490351 5092735143 8650604451 
9077876082 7869878408 8411792866 7296210361 5290067315 
3875194716 5716569026 5501861714 0787645952 9360021510 
7591546225 7321872079 8588558834 8221039668 1905946789 
5755898057 2475290889 3925276026 1115682662 2595034942 
1989465979 9459641436 1822383391 4701051731 6723842532 
9085112016 0805208659 7380301971 0575139883 5276329813 
9051257188 5386696747 9476206813 2405000008 5014000385 
2228602411 7170667435 4071623256 2250293920 8872537198 
9011614093 1973290891 6383126505 4534641138 9553720898 
2096240322 2778414222 8348650739 7672465365 8649120844 
4683221619 8595677967 5204218518 2006601879 7132768680 
9939113203 3472472733 7952464367 0925045889 8276313966 
4012354898 8510115703 2912673298 1529033667 5058577540 
7142956305 7507040532 3467386175 8890359078 9911328777 
6785349324 6971158815 4968664190 8788296990 4226330624 
8768393813 4049426571 5426114747 1205478779 0601988669 
8632035327 0009838500 9172310931 7743393092 6846552160 
2161967107 5702127657 0750291819 1733614937 3394110123 
6021198173 0409567089 7387534596 6112485570 8480721294 
2320621217 6853778480 9999266165 8428469221 5392120811 
5298716572 9275964985 1149125056 3067787170 4139504582 
2131626469 1343571643 5438670288 6985060164 1273634363 
7550116727 3486556597 9439666434 5684036730 1379710374 
4199031303 4300325496 4186410270 3714149988 1526135971 
2877396833 4115425352 6249841425 3502915946 8328043426 
85230425

「標準（筆算）法」や「Karatsuba 法」で計算した結果と一致した。

3. 検証

上記のコードを利用して、乗算回数が何回になるか、計算(100回ループ)にどれくらい時間がかかるかを検証してみた。
参考のため「Karatsuba 法」での乗算テスト結果も掲載している。計算桁数が異なるため単純には比較できないが、1万桁を超えると乗算回数は「Karatsuba 法」より微減し、計算に要する時間も少しずつ短くなっているのが分かる。
だが、やはり C++ 版による計算より莫大に時間がかかっている。
ちなみに、C++ 版では 1000 回ループの時間を計測したが、Ruby 版では時間がかかりすぎるため 100 回ループの時間を計測した。

---

今のところ「Karatsuba 法」や「Toom-Cook 法」を使用する局面はありませんが、理解しておいて損にはならないでしょう。

以上。