2013年02月15日 13:00 [Edit]

Math - おまいら素因数分解どうしてますか?

みなさんは素因数分解の必要にせまられたとき、どうしてますか?

たとえば、こんなとき。

挑戦者求む!【アルゴリズム】チョコの量を減らせ! by The Essence of Programming 結城 浩│CodeIQ
与えられた個数の立方体を組み上げて、できるだけ表面積の小さな直方体を作りましょう

まあ小学生で習うぐらいですし、都度書き下ろしても大した事なさそうにも思えます。

/* works reliably only for 32-bit integer */
var primes = (function(sqrtmax) {
  var result = [2];
  loop: for (var n = 3; n <= sqrtmax; n += 2) {
    for (var i = 0; i < result.length; i++) {
      var p = result[i];
      if (n % p === 0) continue loop;
      if (p * p > n) break;
    }
    result.push(n);
  }
  return result;
})(0xffff);
var factor = function(n) {
  if (n < 2) return undefined;
  var result = [];
  for (var i = 0, l = primes.length; i < l; i++) {
    var p = primes[i];
    while (n % p === 0) {
        result.push(p);
        n /= p;
    }
    if (n === 1) return result;
  }
  if (n !== 1) result.push(n);
  return result;
};
var is_prime = function(n) { return n < 2 ? undefined : factor(n).length === 1 };

var fact = function(n){
    var result = 1;
    while(n){ result *= n--; }
    return result;
};
p(JSON.stringify(factor(fact(13))));
p(JSON.stringify(factor(Math.pow(2,32)+1))); /* Fermat's F5 */

言語によっては標準装備されてたりもします。こんな風に。

require 'prime'
def fact(n) (1..n).inject(:*); end
fact42 = fact(42)
p Prime.prime_division(fact42)

しかし、因数分解すべき数字が78桁だったりしたらどうでしょう?

試してみればわかりますが、たいていはたかだか 64bit 整数でも根を上げます。

% irb1.9 
irb(main):001:0> require 'prime'
=> true
irb(main):002:0> 4294967291.prime?
=> true
irb(main):003:0> 18446744073709551557.prime?
^CIRB::Abort: abort then interrupt!
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:354:in `call'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:354:in `block in succ'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:353:in `loop'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:353:in `succ'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:287:in `block in each'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:286:in `loop'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:286:in `each'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:162:in `prime?'
	from /opt/local/lib/ruby1.9/1.9.1/prime.rb:34:in `prime?'
	from (irb):3
	from /opt/local/bin/irb1.9:12:in `
'

こういう場合みなさんどうしてますか?

msieve

普通のパソコンで動く専用ソフトとしては、これが目下最強そうです。

Msieve | Free Security & Utilities software downloads at SourceForge.net
Msieve is a C library implementing a suite of algorithms to factor large integers. It contains an implementation of the SIQS and GNFS algorithms; the latter has helped complete some of the largest public factorizations known

以下はFreeBSDのportsで入れたものですが、瞬殺ですね。

% msieve -e -p -q -v 280671392065546467397265294532969672241810318954163887187279320454220348884327


Msieve v. 1.50 (SVN unknown)
Fri Feb 15 12:38:27 2013
random seeds: d21ba56a c31d19a5
factoring 280671392065546467397265294532969672241810318954163887187279320454220348884327 (78 digits)
p9 factor: 162425297
p9 factor: 215940091
p9 factor: 358456949
p9 factor: 369941863
p9 factor: 369941863
p9 factor: 479871607
p9 factor: 706170617
prp18 factor: 481362815814826159
elapsed time 00:00:01

PARI/GP

もう少し一般的なソフトとしては、PARI/GPなんかが有名どころです。

PARI/GP Development Headquarters
PARI/GP is a widely used computer algebra system designed for fast computations in number theory (factorizations, algebraic number theory, elliptic curves...), but also contains a large number of other useful functions to compute with mathematical entities such as matrices, polynomials, power series, algebraic numbers etc., and a lot of transcendental functions. PARI is also available as a C library to allow for faster computations.
% gp
                  GP/PARI CALCULATOR Version 2.5.3 (released)
          i386 running darwin (x86-64/GMP-5.0.5 kernel) 64-bit version
     compiled: Feb 12 2013, gcc-4.2 (clang-425.0.24) (based on LLVM 3.2svn)
                 (readline v6.2 enabled, extended help enabled)

                     Copyright (C) 2000-2011 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500509
? factor(280671392065546467397265294532969672241810318954163887187279320454220348884327)
%1 = 
[162425297 1]

[215940091 1]

[358456949 1]

[369941863 2]

[479871607 1]

[706170617 1]

[481362815814826159 1]

? 

こちらも瞬殺でした。

perl - Math::Prime::Util

とはいえ、出来れば普段使っている言語から使えた方がうれしい。PARI/GPはMath::PariというPerl Bindingがすでにあるのですが、いまいち使いづらい。

この手の問題を解くのにちょうどよいものはないかと探して見つけたのが、 Math::Prime::Util でした。llevalでやってもタイムアウトしない程度には優秀ですし、Math::Pariのようにデータのやりとりに特別なオブジェクトを用いることもないので、結果を再利用するのに最適です。

use v5.14;
use bignum;
use Math::Prime::Util qw/:all/;
say for factor(280671392065546467397265294532969672241810318954163887187279320454220348884327);

他の言語だとどうなんだろう。Pythonとか得意そうなんだけど、ちょっとぐぐったぐらいではわからなかった。一番上のJavaScriptの例程度のやつならいくらでも見つかるのだけど。「ワシの言語ではこうやってる」という事例を緩募。

追記:

Wolfram|Alpha

しかし圧巻なのはWolfram|Alphaでしょう。なんと数字を入れるだけ。それが整数なら頼みもしないのに素因数分解してくれちゃいます。

明示的に素因数分解してほしかったら、数字の前に'factor'と入れるだけ。

あとは角砂糖を積み上げてチョコレートを塗るだけなのですが、実行する前にこれをよく見てください。

280671392065546467397265294532969672241810318954163887187279320454220348884327 - Wolfram|Alpha
≈ 0.0028 × the number of atoms in the visible universe (≈ 10^80)

(つд⊂)ゴシゴシ…え?えーっ!

全宇宙の元素の0.0028倍!?

ショ糖の分子量は342。仮に角砂糖が分子一個でできていたとしても、宇宙がまるまる一個必要になる計算です。

愛は地球どころか、宇宙より重かった…

Dan the Prime Blogger


この記事へのトラックバックURL

この記事へのコメント
>coreutil なのにgmpに依存するというのはいかがなもんかなあ…
おそらく2**64未満ならばgmpとリンクされていないビルドでも素因数分解できると思います。

確かに、factorがcoreutilsに入った経緯とか調べると面白いかもしれませんね--;
Posted by https://me.yahoo.com/a/zddE0JA.xd5_j1Kf_iW3aSCTxohwwbzE_ug-#b7178 at 2013年02月20日 03:30
> #define HAVE_GMP 1でビルドされたcoreutils factorなら通ります。
coreutil なのにgmpに依存するというのはいかがなもんかなあ…
そもそもcoreにふさわしいのかなfactor(1)って…

Dan the Man with Too Many Dependencies to Resolve
Posted by dankogai at 2013年02月19日 09:33
twitterにリプ致しましたが、見落とされたのでしょうか…。

/opt/local/bin/gfactor 280671392065546467397265294532969672241810318954163887187279320454220348884327
280671392065546467397265294532969672241810318954163887187279320454220348884327: 162425297 215940091 358456949 369941863 369941863 479871607 706170617 481362815814826159

#define HAVE_GMP 1でビルドされたcoreutils factorなら通ります。Linuxだとソースからリビルドが必要かもしれません(Debainでは必要でした)。

cielavenir.
Posted by https://me.yahoo.com/a/zddE0JA.xd5_j1Kf_iW3aSCTxohwwbzE_ug-#b7178 at 2013年02月19日 01:38