スキップしてメイン コンテンツに移動

Project Euler - Problem 3

問題

  • 原文

    What is the largest prime factor of the number 600851475143 ?

  • 日本語訳

    600851475143 の素因数のうち最大のものを求めよ。

解答

いきなり難易度が上がりました。素因数分解の問題です。 結局のところ片っ端から割ってみるしかないのですが、探索範囲を狭めることはできます。

正の整数の範囲で考えると、l = m * nかつm >= nであれば、floor(sqrt(l)) >= nであることが直感的に分かります。ここでfloorは実数を0の方向に丸めて整数にする関数、sqrtは平方根を返す実数関数です。 nが分かれば、mは除算で求められます。つまりlが与えられた時、それを2つの因数に分解するのに探索する範囲は高々2 .. floor(sqrt(l))で十分ということです。 2つに分けられたならこっちのもの、得られたmとnを再帰的に分解して、その結果を合わせれば素因数分解の結果が得られます。 アルゴリズムを書き下すと次のようになります:

  1. n <- floor(sqrt(l))とする。
  2. n = 1なら、これ以上は分解できない(lは素数である)のでlをそのまま返す。
  3. nがlの因数なら、m <- l / nとし、mとnに対してこのアルゴリズムを再帰的に適用し、その結果を連結して返す。
  4. 因数でないなら、n <- n - 1とし、2.から繰り返す。

これで計算は可能です。しかしもう少し最適化を検討してみましょう。

偶数は一般に2mの形で表せます。同様に奇数は2m + 1と書けます。 偶数同士の積は2m * 2n = 4mnとなり、l = 2mnとおけば2lなのでこれも偶数です。 奇数同士の積は(2m + 1)(2n + 1) = 4mn + 2(m + n) + 1であり、l = 2mn + m + nとすると2l + 1なのでこれも奇数です。 偶数と奇数の積は2m(2n + 1) = 4mn + 2mなので、l = 2mn + mとおいて2lなので偶数です。 つまり、奇数の因数は必ず奇数のみから成るということが分かります。先ほどのアルゴリズムでは、4.のステップでn <- n - 1としていました。これではlとnが両方とも奇数の時、1回不要な繰り返しが生じることになります。つまりこの場合を特別扱いして、n <- n - 2とすれば繰り返しを減らすことができます。 幸い、問題で与えられた600851475143という数値は奇数なので、この最適化の恩恵を十分に受けることができます。

この最適化を施したところ、40%ほど高速化しました。上の考察におけるlはn、nはdivという名前になっています。

(define (divisor-of? n m) (zero? (modulo m n)))
(define (factorize n)
  (define div (inexact->exact (floor (sqrt n))))
  (define (factorize-1 n div)
    (cond
     ((= div 1) (list n))
     ((and (odd? n) (even? div)) (factorize-1 n (- div 1)))
     ((divisor-of? div n) (append (factorize (/ n div)) (factorize div)))
     ((odd? n) (factorize-1 n (- div 2)))
     (else (factorize-1 n (- div 1)))))
  (factorize-1 n div))
(define (solve)
  (apply max (factorize 600851475143)))
(define (main argv)
  (display (solve))
  (newline))

追記

コメントの匿名氏のコードを基にしたところ500倍高速になりました。よって上記の議論は忘れて手続き的に書いた方が良さげ。 多値を扱うのにSRFI-8のreceiveマクロを使っています。

(define (divisor-of? m n) (zero? (modulo n m)))
(define (factorize n)
  (define (factorize-1 n i result)
    (cond
     ((< n (* i i)) (reverse (if (= n 1) result (cons n result))))
     ((divisor-of? i n) (factorize-1 (/ n i) i (cons i result)))
     (else (factorize-1 n (+ i 2) result))))
  (receive (n result)
           (let prepare-loop ((n n) (result '()))
             (if (divisor-of? 2 n) (prepare-loop (/ n 2) (cons 2 result))
                 (values n result)))
           (factorize-1 n 3 result)))
(define (solve)
  (apply max (factorize 600851475143)))
(define (main argv)
  (display (solve))
  (newline))

コメント

  1. (define factorize
       (lambda (n)
     
          (define result '()
     
          (define check
             (lambda (n i)
                (let loop ([n n] [c 0])
                   (if [zero? (remainder n i)]
                         (loop (/ n i) (+ c 1))
                         (begin
                            (if [> c 0] (set! result (cons (cons i c) result)))
                            n)))))
     
          (let loop ([n (check n 2)] [i 3])
             (cond
                ([= n 1] (reverse result))
                ([> (* i i) n] (reverse (cons (cons n 1) result)))
                (else (loop (check n i) (+ i 2)))))))
     
    (caar (reverse (factorize 600851475143)))

    返信削除
  2. 速い!当社比400倍ですね。
    checkが何をしているのか少し悩みました。乗数をカウントして連想リストにしているのか。
    副作用を使ったコードは避けたいのですが、ここまで違うと参るなあ。参考にさせてもらいます。

    返信削除

コメントを投稿

このブログの人気の投稿

部分継続チュートリアル

この文書についてこれはCommunity Scheme Wikiで公開されているcomposable-continuations-tutorial(2010年09月30日版)の日本語訳です。誤字脱字・誤訳などがありましたらコメントあるいはメールで御指摘いただけると幸いです。本訳は原文のライセンスに基づきCreative Commons Attribution-ShareAlike 2.0 Genericの下で公開されます。Original text: Copyright© 2006-2010 Community Scheme WikiJapanese translation: Copyright© 2011 SATOH Koichi本文部分継続(Composable continuation)は継続区間を具象化することで制御を逆転させるものです。 ウンザリするほど複雑な概念を表す長ったらしいジャーゴンのように聞こえますが、実際はそうではありません。今からそれを説明します。resetとshiftという2つのスペシャルフォームを導入するところから始めましょう[1]。 (reset expression)は特別な継続を作るなりスタックに目印を付けるなりしてからexpressionを評価します。簡単に言えば、expressionが評価されるとき、あとから参照できる評価中の情報が存在するということです。 実際にはshiftがこの情報を参照します。(shift variable expression)は目印のついた場所、つまりresetを使った場所にジャンプし、その場所からshiftを呼び出した場所までのプログラムの断片を保存します; これはプログラムの区間を「部分継続」として知られる組み合わせ可能な手続きに具象化し、この手続きにvariableを束縛してからexpressionを評価します。組み合わせ可能(Composable)という語はその手続きが呼び出し元に戻ってくるため、他の手続きと組み合わせられることから来ています。 Composable continuationの別名として例えば限定継続(Delimited continuation)や部分継続(Partial continuation)もありますが、ここでは一貫して「組み合わせ可能」という用語を使います(訳注: …

Perl 5 to 6 - 列挙型

これはMoritz Lenz氏のWebサイトPerlgeek.deで公開されているブログ記事"Perl 5 to 6" Lesson 16 - Enumsの日本語訳です。原文はCreative Commons Attribution 3.0 Germanyに基づいて公開されています。本エントリにはCreative Commons Attribution 3.0 Unportedを適用します。Original text: Copyright© 2008-2010 Moritz LenzJapanese translation: Copyright© 2011 SATOH KoichiNAME"Perl 5 to 6" Lesson 16 - 列挙型SYNOPSISenum bit Bool <False True>; my $value = $arbitrary_value but True; if $value { say "Yes, it's true"; # 表示される } enum Day ('Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'); if custom_get_date().Day == Day::Sat | Day::Sun { say "Weekend"; } DESCRIPTION列挙型は用途の広い獣です。定数の列挙からなる低レベルのクラスであり、定数は典型的には整数や文字列です(が任意のものが使えます)。これらの定数は派生型やメソッド、あるいは通常の値のようにふるまいます。 but演算子でオブジェクトに結びつけることができ、これによって列挙型を値に「ミックスイン」できます:my $x = $today but Day::Tue; 列挙型の型名を関数のように使うこともでき、引数として値を指定できます:$x = $today but Day($weekday); こうするとオブジェクトは列挙型の名前(ここではDay)をメソッド名として持ちます:say $x.D…

去る6月に Perl 5.32.0 がリリースされたので差分を把握するために perldelta を読んだ件

要旨Perl 5 メジャーバージョンアップの季節がやって来たのでまともな Perl プログラマの嗜みとして perldelta を読んだ。今回は有り体に言えばルーティン的なリリースで、言語コアの拡張は他言語にも見られる構文が実験的に入ったくらいで大きな変化はない。新機能は RegExp の拡充が主である。比較的重要と思われる変更点を抜粋する。新機能isa 演算子実験的機能。Python とか Java における isinstance とか instanceof。これまでも UNIVERSAL::isa があったが、これはメソッドなのでレシーバにオブジェクトでもクラスでもない値 (i.e., 未定義値 / bless されていないリファレンス) を置くと実行時エラーが起きるのが問題だった:package Foo {use Moo;}package Bar {use Moo; extends 'Foo';}package Baz {use Moo;}use feature qw/say/;sub do_something_with_foo_or_return_undef {my ($foo) = @_;# Returns safely if the argument isn't an expected instance, in mind.returnunless$foo->isa('Foo'); ...;}# OK.do_something_with_foo(Bar->new);# |undef| is expected in mind, but actually error will be thrown.do_something_with_foo(undef); これを避けるために今までは Scalar::Util::blessed を併用したりしていたわけだが、isa 演算子は左辺が何であっても意味のある値を返すのでよりシンプルになる:# True+(bless +{} => 'Foo') isa Foo;# Falseundef isa Foo;# False+{} isa Foo;# False+(bless +{} => 'Baz') isa Foo; 比較演算子…