素因数分解(Ruby)

素因数分解は結局素数で順に割っていくしかないのだけれど、素数をわざわざ求めてというのは却って大変ですね。

うまいやり方としては、2 および 3 以上の奇数で割るという方法があります。これは多少のムダは出るけど、全部の数で割るより効率的ですね。さらにうまいやり方があって、2, 3 以上は 5, 7, 11, 13, 17, 19, 23, 25, 29 .... と、ステップ幅を 2 と 4 の交代にして割っていくという方法があります。かしこいですね。これは例えば Ruby でこんなふうに書けます。
prime_factorize.rb

def prime_factorize(n)
  if n <= 1 or !n.class.ancestors.include?(Integer)
    raise ArgumentError, "#{n} incorrect."
  end
  result = {}
  divide = ->(i) {
    while (n % i).zero?
      result[i] = result.has_key?(i) ? result[i] + 1 : 1
      n /= i
    end
  }
  divide[2]
  divide[3]
  k, step = 5, 2
  guard = Math.sqrt(n).to_i
  while n != 1
    divide[k]
    k += step
    return result if k > guard    #高速化
    step = 6 - step
  end
  result
end

if __FILE__ == $0
  p prime_factorize(15141523320)    #=>{2=>3, 3=>3, 5=>1, 7=>2, 11=>1, 19=>1, 37=>2}
end

結果はハッシュで返ります。

なお、素因数分解Ruby では標準添付ライブラリにあるので、実際はそれを使えばよいですね。ちなみに例えば 69316536495422 の素因数分解をやると、上のプログラムではフリーズしますが(笑)、ライブラリのメソッドだと苦もなく瞬時に素因数分解します。
※追記:2行付け加えただけで劇的に高速化しました。これだと 69316536495422 の素因数分解でも一瞬です。(4/29))

 

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

L-system で描いた再帰曲線の例

obelisk.hatenablog.com
前回のエントリで Ruby で実装してみた L-system を使って、実際に再帰曲線を描いてみます。Gem 'kaki-lsystem' が必要です。


C曲線。
20180323002343
コード。

require 'kaki/lsystem'

l = Lsystem.new(600, 500, "C curve")
l.move(-130, 100)
l.prologue do
  clear(color(0x2b5, 0x60ff, 0x5d89))
  color(0xfc9e, 0xe0aa, 0x7ac7)
end
l.set("-") {left(45)}
l.set("+") {right(45)}
l.set("F") {forward(8)}
l.init("F")
l.rule("F", "+F--F+")
l.draw(10)

以下、require 'kaki/lsystem' は書きません。


コッホ・アイランド。
20180322043252
コード。

l = Lsystem.new(600, 600, "Koch island")
l.move(-150, -150)
l.prologue do
  clear(color(0xffff, 0xf600, 0xf540))
  color(0x28a5, 0xa3c8, 0xca7)
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("F") {forward(8)}
l.init("F-F-F-F")
l.rule("F", "F+FF-FF-F-F+F+FF-F-F+F+FF+FF-F")
l.draw(2)

 

ペンローズ・タイル
20180322042852
コード。

l = Lsystem.new(600, 600, "Penrose Tiling")
l.prologue do
  clear(color(0xe013, 0xfee3, 0xfee3))
  color(0xffff, 0x5feb, 0xe51f)
end
l.set("-") {left(36)}
l.set("+") {right(36)}
walk = proc {forward(10)}
l.set("M", &walk)
l.set("N", &walk)
l.set("O", &walk)
l.set("P", &walk)
l.set("A", &walk)
l.init("[N]++[N]++[N]++[N]++[N]")
l.rule("M", "OA++PA----NA[-OA----MA]++")
l.rule("N", "+OA--PA[---MA--NA]+")
l.rule("O", "-MA++NA[+++OA++PA]-")
l.rule("P", "--OA++++MA[+PA++++NA]--NA")
l.rule("A", "")
l.draw(5)

 

Sierpinski arrowhead。続けると「シェルピンスキーのギャスケット」になります。
20180324150434
コード。

n = 6
length = n.zero? ? 500 : 500.0 / 2 ** n

l = Lsystem.new(520, 500, "Sierpinski arrowhead")
l.move(-250, -240)
l.prologue do
  clear(color(0xffff, 0xef5b, 0xf677))
  color(0x12de, 0x5320, 0xb4a)
end
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("L") {}
l.set("R") {}
l.set("F") {forward(length)}
l.init("LF")
l.rule("L", "-RF+LF+RF-")
l.rule("R", "+LF-RF-LF+")
l.rule("F", "")
l.draw(n)

 
 
20180324150916
コード。

l = Lsystem.new(500, 500)
l.move(50, -230)
l.prologue do
  clear(color(0xd467, 0xfbab, 0xd0c7))
  color(0xe713, 0x5742, 0xff23)
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("F") {forward(10)}
l.init("F-F-F-F")
l.rule("F", "FF-F-F-F-F-F+F")
l.draw(3)

 
 
草っぽいの。
20180323002349
コード。

l = Lsystem.new(500, 500)
l.move(-150, -240)
l.prologue do
  clear(color(0xf7d3, 0xffff, 0xc99b))
  color(0x30e3, 0x8636, 0x10d2)
end
l.dir = Vector[0, 1]
l.set("-") {left(10)}
l.set("+") {right(15)}
l.set("F") {forward(13)}
l.set("X") {}
l.set("Z") {}
l.init("X")
l.rule("F", "FX[FX[+XF]]")
l.rule("X", "FF[+XZ++X-F[+ZX]][-X++F-X]")
l.rule("Z", "[+F-X-F][++ZX]")
l.draw(4)

 

ペアノ・ゴスペル曲線。
20180325031914
コード。

l = Lsystem.new(600, 550, "Peano-Gosper curve")
l.move(80, 250)
l.prologue do
  clear(color(0x5968, 0x5968, 0x5968))
  color(0x8024, 0xffff, 0x653)
end
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("F") {forward(9)}
l.set("X") {}
l.set("Y") {}
l.init("X")
l.rule("X", "X+YF++YF-FX--FXFX-YF+")
l.rule("Y", "-FX+YFYF++YF+FX--FX-Y")
l.draw(4)

再帰曲線を描く言語「L-system」を Ruby で実装した

L-system って何でしょうか。Wikipedia にはこうあります。

L-system(エルシステム、Lindenmayer system)は形式文法の一種で、植物の成長プロセスを初めとした様々な自然物の構造を記述・表現できるアルゴリズムである。自然物の他にも、反復関数系(Iterated Function System; IFS)のようないわゆる自己相似図形やフラクタル図形を生成する場合にも用いられる。

https://ja.wikipedia.org/wiki/L-system

 

つまりは、再帰曲線を描く言語ですね。例としてまず「コッホ曲線」を描いてみましょう。
まず、前に一定の長さ進むコマンドを「F」とします。で、開始手続きも「F」としましょう。これだけで
20180322132620
と直線が一本引かれます。
つぎに、再帰的な「置き換え」の規則を定めます。それを「F→F-F++F-F」とします。これは開始手続きの「F」を「F-F++F-F」で置き換えるということになります。「-」は左60度、「+」は右60度の回転とします。よって、新たな手続きも「F-F++F-F」となります。これを図示すると
20180322132618
となります。
置き換えを続けましょう。今度は前の手続き「F++F-F」の F をすべて F++F-F に置き換えます。すると「F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F」となります(ややこしいですね)。これは
20180322132615
となります。
もう一度置き換えてみます。今度は「F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F-F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F-F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F」ですよ。これは
20180322132613
となりますね。まさしく再帰的な「コッホ曲線」が描けています。


さて、これはこんなコードで描けます。Gem化しているので、Gem 'kaki-lsystem' をインストールして下さい。

require 'kaki/lsystem'

n = 3
length = 500.0 / 3 ** n

l = Lsystem.new(550, 200)
l.move(-250, -70)
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("F") {forward(length)}
l.init("F")
l.rule("F", "F-F++F-F")
l.draw(n)

init("F") は開始手続きの指定、rule("F", "F-F++F-F") は置き換えの規則の指定です。draw(3) で変換を 3回施したものを描画します。

set("-") {left(60)} はコマンド「-」が何をするのか、ブロック内に記述します。つまり left(60) は左60度回転ということですね。forward(20) ならまっすぐ 20 進むということです(ちなみに forward(20, false) なら、描画せずに 20 進みます)。このブロック内は、じつは Gem 'oekaki' のタートルグラフィックスの機能を使っているので、そのメソッドを書きます。よく使うのは left(deg), right(deg), forward(length, draw=true) あたりでしょうか。これについては

などを参考にして下さい。なお、set() は使用したすべてのコマンドについて設定して下さい。何もしない場合でも空のブロックが必要です。また、rule() で変換規則が指定されていないコマンドは、そのまま持ち越されます。


置き換えるコマンドを二種類以上にすることもできます。例えば

require 'kaki/lsystem'

l = Lsystem.new(400, 350, "Dragon curve")
l.move(-50, -50)
l.prologue do
  clear(color(0xf7d3, 0xffff, 0xc99b))    #背景色
  color(0x143f, 0x832d, 0xe783)           #ペンの色
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("A") {forward(20)}
l.set("B") {forward(20)}
l.init("A")
l.rule("A", "A+B+")
l.rule("B", "-A-B")
l.draw(7)

では、開始手続きが A で、置き換えの規則は「A→A+B+」「B→-A-B」と二種類になっています。これを実行すると
20180322141013
となります。いわゆる「ドラゴン曲線」ですね。他のメソッドも簡単に説明しておきましょう。move(x, y) はペンの開始位置を指定します。デフォルトは (0, 0) で、これはキャンバスの中心です。なお、数学のグラフと同じで y方向は上がプラスの向きです。prologue {...} はブロック内で初期処理をします。背景色や色の指定をすることが多いです。なお、ここでは残念ながら left(), right() が使えない(というか2度呼ばれてしまうので正確に指定できません)ので、ペンの最初の方向はアクセサーメソッド dir= で指定することになります(あとで説明します)。


コマンド「[」「]」でスタックのプッシュ、ポップを表します。つまり [ で状態を保存し、 ] で元の状態に復帰します。例えば

require 'kaki/lsystem'

l = Lsystem.new(400, 640)
l.move(50, -310)
l.prologue {color(0x28a5, 0xa3c8, 0xca7)}
l.dir = Vector[0, 1]
l.set("-") {left(15)}
l.set("+") {right(15)}
l.set("F") {forward(10)}
l.init("F")
l.rule("F", "F[+F-F-F]F[--F+F+F]")
l.draw(4)

こんな風です。dir= で開始時の方向をベクトルで指定しています。ただし、ベクトルは単位ベクトル(長さが 1 )になるようにして下さい。この結果は
20180322143459 20180322143456 20180322142514
となります。左からそれぞれ 1, 2, 4 回の繰り返しになります。


ブロックはクロージャなので、変数を使ったサンプル。「ヒルベルト曲線」の色を描画しながら替えていきます。

require 'kaki/lsystem'

n, f = 0, 0
col = [[0xffff, 0, 0], [0, 0xffff, 0]]

l = Lsystem.new(500, 500, "Hilbert curve")
l.move(-230, 230)
l.set("+") {left(90)}
l.set("-") {right(90)}
l.set("A") {}
l.set("B") {}
l.set("F") do
  if n.zero?
    f = (f + 1) % 2
    color(*col[f])
  end
  forward(15)
  n = (n + 1) % 16
end
l.init("A")
l.rule("A", "-BF+AFA+FB-")
l.rule("B", "+AF-BFB-FA+")
l.draw(5)

結果。
20180322144224
 

Gem 'kaki-lsystem' のソースコードは Gist にあります。この Gem は内部で Gem 'oekaki' を使っています。
https://gist.github.com/obelisk68/1c3d278cb5fe2778a11ba719d44af743


※参考にしたページ(ありがとうございます!)
L-system
L-system入門 - 人工知能に関する断創録
HiiragiCompany -L-System Tips-
L-systemを利用した植物の描画 | Makoto Hirahara
 
 
※L-system の描画例
obelisk.hatenablog.com

ジュリア集合を描いてみる(Ruby)

ジュリア集合(Julia set)を Ruby で描いてみました。
20180318023736
 
描画には自作の Gem 'oekaki' を使っています。
oekaki | RubyGems.org | your community gem host
GTK+でお絵かきしてみた(Ruby) - Camera Obscura
 
コード。
julia_set.rb

require 'oekaki'

Width, Height = 600, 400

Oekaki.app width: Width, height: Height, title: "Julia set" do
  draw do
    clear
    red   = color(0xffff, 0, 0)
    green = color(0, 0xffff, 0)
    blue  = color(0, 0, 0xffff)
    
    a, jmid = 4.0 / Width, Height / 2
    Width.times do |i|
      jmid.times do |j|
        x, y = a * i - 2.0, a * j
        catch(:jump) do
          16.times do
            x2, y2 = x * x, y * y
            d = x2 + y2
            break if d < 1e-10
            d = d * d
            x, y = (1.0 / 3) * (2 * x + (x2 - y2) / d), (2.0 / 3) * y * (1 - x / d)
            ya = sqrt(3) / 2
            
            drawing = proc do |x1, y1, c1, c2|
              if x1 * x1 + y1 * y1 < 0.0025
                point(i, jmid + j, c1)
                point(i, jmid - j, c2)
                throw(:jump)
              end
            end
            
            drawing[x - 1, y, red, red]
            drawing[x + 0.5, y + ya, green, blue]
            drawing[x + 0.5, y - ya, blue, green]
          end
        end
      end
    end
  end
end

 

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

これおもしろい本です。おすすめ!

あみだくじの横線(Ruby)

アルゴリズム・パズルです。

問題:
1 ~ 7 の数字を並べてあみだくじを作ることを考えます。引く横線の数を 10本とするとき、下に得られる数字の並びは何とおりになるでしょうか。ただし、下に得られる数字の並びが 10本より少ない横線のあみだくじで既に得られる場合、それは除くことにします。

 
Ruby で解いてみました。
q57.rb

require 'set'

N = 7
List = "1234567"

def exchange(i, ls, num)
  ls[i], ls[i + 1] = ls[i + 1], ls[i]
  if num == 10
    @result << ls
  else
    @store << ls
    (N - 1).times do |j|
      next if i == j
      exchange(j, ls.dup, num + 1)
    end
  end
end

@store = Set[List]
@result = Set.new
(N - 1).times {|i| exchange(i, List.dup, 1)}
puts (@result - @store).size    #=>573

答えはいちおう 573通りと求められたのですが、時間が約 16秒もかかってしまいました。exchange() の呼び出しは 14648436回もあります。力任せのやり方で、これではダメですね。

模範解答の一。

v, h = 7, 10
total = 0

[*0..v - 1].permutation do |final|
  cnt = 0
  v.times do |i|
    cnt += final.take_while {|j| j != i}.count {|k| k > i}
  end
  total += 1 if cnt == h
end
puts total

超シンプルなコードなのだが、これで 0.1秒。さらに高速なアルゴリズムも紹介されている…。すごいな。
 

公平に分けられたケーキ(Ruby)

アルゴリズム・パズルです。

問題:
16cm × 12cm の長方形のケーキがあります。二人で交互にケーキを切るのですが、ケーキを切った人が小さい方のケーキを取り、交代して残りをさらに切るというようにします。(ただし切り方は、辺に平行に直線的に切り、しかも切ってできたケーキの辺の長さは cm 単位になるようにします。)これを最小が 1cm × 1cm になるまで続けます。
さて、切り終えたとき二人の取ったケーキのそれぞれの総量が等しくなるとします。このような切り方は何とおりもありますが、その中で切った「切り口」の長さの総和を出すとき、その最小値を求めなさい。

 
Ruby で解いてみました。コードは以下です。
q56.rb

require 'set'

#(x, y)のケーキを切ってできる全ての (面積, 切った長さの合計) のペアを返す
#ただし面積はここでカットした人の総計
def cut(x, y)
  x, y = y, x if x < y
  return @memo[[x, y]] if @memo.has_key?([x, y])
  return Set[[1, 1]] if x == 2 and y == 1
  
  result = Set.new
  
  #前の人がカットしたところから自分の分を求める
  calc = proc do |s, t|
    cut(s, t).each do |sqar, ln|
      square = x * y - sqar
      next if square > X * Y / 2    #高速化
      length = ln + s
      result << [square, length]
    end
  end
  
  #可能なあらゆる場合でケーキを切る
  (x / 2).times {|i| calc.call(y, x - i - 1)}
  (y / 2).times {|i| calc.call(x, y - i - 1)}
  
  @memo[[x, y]] = result
end

X, Y = 16, 12
@memo = {}

result = cut(X, Y)
p result.select {|sqar, _| sqar == X * Y / 2}.sort_by {|_, ln| ln}[0]    #=>[96, 47]

答えは 47cm です。これはかなりむずかしかったです。あまりいい解き方とも思えません。実行時間は自分の環境で 1.5秒ほどでした。もっとも深いループの回数は 1386985回です。

なお、標準添付ライブラリの set を使っていますが、これを配列でやると飛んでもないことになってしまいます。


※追記(3/7)
ハッシュを使ってほんの少し書き替えたところ、劇的に高速化しました。実行時間は 0.05秒ほどになりました。じつに 30倍の高速化です。アルゴリズムはまったく同じですが、探索の際にうまく枝刈りをする方法を思いつきました。もっとも深いループの回数は 33192回なので、ほぼ 1 / 40 になっています。コードはこちら。30cm × 30cm の場合でも 0.5秒足らずで終了します。
 
※再追記(4/14)
あとから見なおしたらよくこんなこと思いついたなと思いました。関数 calc で前の面積から次の面積を求めるところが我ながら巧妙かと感じます。
 

GTK+ でカラーパレットを作る(Ruby)

20180305002252
いつも色がうまく選べないので、自分用に GTK+ を使って色を試すことのできるコマンドを作ってみました。Gem 'thor', 'gtk2' が必要です。

コード。
color_p

require 'thor'
require 'gtk2'

class MyWindow < Gtk::Window
  def initialize
    super("color palette")
    set_default_size(250, 250)
    set_resizable(false)
    
    colorsel = Gtk::ColorSelection.new
    entry  = Gtk::Entry.new
    
    vbox = Gtk::VBox.new
    add(vbox)
    vbox.pack_start(colorsel, false, true, 0)
    vbox.pack_start(entry, true, true, 0)
    colorsel.current_color = Gdk::Color.new(0, 0, 0)
    
    # カラーパレット変更時
    colorsel.signal_connect("color-changed") do |w|
      c = w.current_color
      t = "(#{c.red.to_s(16)}, #{c.green.to_s(16)}, #{c.blue.to_s(16)})"
      entry.set_text(t)
    end
    
    signal_connect("destroy") {Gtk.main_quit}
    show_all
  end
end

class ColorPalette < Thor
  default_command :select
  desc 'select', 'Run color palette.'
  def select
    w = MyWindow.new
    Gtk.main
  end
end

ColorPalette.start

 
下のサイトを参考にさせて頂きました。ありがとうございます。
Ruby Entry markup: いち雑記