Shoeisha Technology Media

CodeZine(コードジン)

記事種別から探す

<numeric>がらみの数値演算に触発されて、「巨大整数のかけ算」を行う方法

フーリエ変換と畳み込み定理

  • LINEで送る
  • このエントリーをはてなブックマークに追加
2015/08/24 14:00

 前回<numeric>をひとめぐりして以来数値演算が面白く感じられ、Webの海に漁に出かけてました。その中でちょっと気になるpaperを見つけましてね:「Big Integer Multiplication with CUDA FFT(cuFFT) Library」。CUDA FFT(cuFFT)を使って巨大整数の掛け算を行うおはなし。10pageにも満たないpaperながら計算の手順が分かりやすく述べられていて、ものは試しとcuFFTのかわりにFFTWで実装してみたらあっさりできちゃいました。数万桁の巨大な整数の掛け算を高速に行う方法をご紹介します。

目次

紙と鉛筆での掛け算

 僕の愛用するVisual C++で扱える最大の整数はunsigned long long:符号なし64bit整数です。1bitつまり2進での1桁は10進では約0.3桁に相当するので、unsigned long longで表現できる最大の整数は10進で(64*0.3=)19桁、9桁どうしの掛け算まではオーバーフローせずに求めることができます。……さて、それ以上に大きな数千/数万桁の整数の掛け算を正しく行うにはどうしましょうね、と。

 整数の掛け算を紙と鉛筆で行う"筆算"は小学校で習いました。計算機向けにちょっとアレンジしておさらいします。

              4   5   6
 ×               7   8
-------------------------
             32  40  48   (400 + 50 + 6)× 8
 +      28  35  42       (400 + 50 + 6)×70
-------------------------
         28  67  82  48  
      3   5   5   6   8  繰り上がり

 456×78を筆算で計算するには、まずそれぞれの数を桁ごとにバラし、総当たりで掛け合わせたのち各桁の和を求めます。出てきたコタエは 28 67 82 48、各桁が10未満となるよう、1の位から順に繰り上げを行うと 3 5 5 6 8、これが掛け算の結果です。

 小学校の先生は僕たちにすばらしい知恵を授けてくれました。紙と鉛筆(そして時間と根気)さえあればどんなに大きな整数の掛け算だろうが、このやり方で答えが求まります。

 筆算による掛け算を実装してみましょう、入力となる整数と返り値(結果)は10進表記での文字列:std::stringとします。

list-1 onpaper_multiplies.cpp
#include <vector>
#include <string>
#include <algorithm>

using namespace std;

string onpaper_multiplies(const string& sa, const string& sb) {

  int digits = (int)(sa.size() + sb.size());

  vector<int> a(sa.size(),0);
  vector<int> b(sb.size(),0);
  vector<int> c(digits,0);

  // 下位の桁から順にセット
  auto char2int = [](char ch) { return ch - '0';}; 
  transform(sa.rbegin(), sa.rend(), a.begin(), char2int); 
  transform(sb.rbegin(), sb.rend(), b.begin(), char2int); 

  // a,b各桁の全組み合わせで掛け算を行う
  for ( size_t i = 0; i < a.size(); ++i ) {
    for ( size_t j = 0; j < b.size(); ++j ) {
      c[i+j] += a[i] * b[j];
    }
  }

  // 桁あふれの繰り上げ
  int carry = 0;
  transform(c.begin(), c.end(), c.begin(), 
                 [&](int n) { int t = n + carry; carry = t /10; return t % 10; });

  // 文字列化(ゼロ・サプレスを行う)
  string result;
  result.reserve(digits);
  bool zero = false;
  for_each(c.rbegin(), c.rend(), 
                [&](int ch) { if ( ch != 0 && !zero ) { zero = true; }
                              if ( zero ) { result.push_back(ch+'0'); } });
  return result;
}

#ifdef TRIAL
#include <iostream>
using namespace std;

int main() {
  string sa = "456";
  string sb =  "78";
  string sc = onpaper_multiplies(sa, sb);
  cout << sa << " * " << sb << " = " << sc << endl;   
}
#endif

 正しい結果が得られたみたいだけど、念のためそこそこ大きな数を食わせて確認しておきます。.NET Framework: System.Numerics.BigInteger の助けを借りて2つの整数とその積を生成し、onpaper_multiplies での結果とを照合します。

list-2 make_random_multiplies.cpp(C++/CLI)
#include <string>
#include <random>
#include <algorithm>
#include <msclr/marshal_cppstd.h>

using namespace System;
using namespace System::Numerics;
using namespace std;
using namespace msclr::interop;

int main(array<System::String ^> ^args) {

  if ( args->Length < 2 ) {
    Console::WriteLine(L"make_random_multiplies <digits> <rows>");
    return 1;
  }
  int digits = Int32::Parse(args[0]); // digits桁の2整数とその積を
  int rows   = Int32::Parse(args[1]); // rows組出力する
  mt19937 gen;
  uniform_int_distribution<> dist(0,9);
  auto random_digit = [&]() -> wchar_t { return dist(gen) + L'0'; };

  wstring sa(digits, L'0');
  wstring sb(digits, L'0');

  for ( int i = 0; i < rows; ++i ) {
    generate_n(begin(sa), digits, random_digit);
    generate_n(begin(sb), digits, random_digit);
    BigInteger a = BigInteger::Parse(marshal_as<String^>(sa));
    BigInteger b = BigInteger::Parse(marshal_as<String^>(sb));
    BigInteger c = BigInteger::Multiply(a, b);
    Console::WriteLine(L"{0} {1} {2}", a.ToString(), b.ToString(), c.ToString());
  } 

  return 0;
}
list-3 test_onpaper_multiplies.cpp
#include <iostream>
#include <string>

using namespace std;

// 所要時間の計測を行うよう、少し修正した
string onpaper_multiplies(const string& sa, const string& sb, float& duration);

int main() {
  string sa, sb, expected, actual;
  float elapsed;
  float duration = 0.0f;
  int times = 0;
  while ( cin >> sa >> sb >> expected ) {
    actual = onpaper_multiplies(sa, sb, elapsed);
    duration += elapsed;
    if ( expected != actual ) { 
      cerr << L"oops!" << endl; 
      return 1; 
    }
    ++times;
  }
  cerr << "ok. takes " << duration << "[ms] for " << times << " multiplies." << endl;
  return 0;
}
fig-1
fig-1

 1000桁の掛け算を一万回やって同じ結果が得られてますから、間違っちゃいないみたいです。

 筆算アルゴリズムで書かれたonpaper_multipliesには1つ大きな問題があります、それは処理時間。

 この計算法のキモとなるのは各桁総当たりの掛け算部分。見てのとおり二重のfor-loopですから、N桁どうしの掛け算に必要な時間はNの二乗に比例します。時間計算量:O(N^2) てことです。

 Nが小さいうちはどってことないのですが、大きくなるにつれてその遅さに我慢できなくなります。マルチスレッド化したところで高々数倍、多少マシにはなりますが O(N^2)には歯が立ちません。


  • LINEで送る
  • このエントリーをはてなブックマークに追加

著者プロフィール

  • επιστημη(エピステーメー)

    C++に首まで浸かったプログラマ。 Microsoft MVP, Visual C++ (2004.01~) だったり わんくま同盟でたまにセッションスピーカやったり 中国茶淹れてにわか茶人を気取ってたり、 あと Facebook とか。 著書: - STL標準講座 (監修) -...

おすすめ記事

All contents copyright © 2006-2017 Shoeisha Co., Ltd. All rights reserved. ver.1.5