コッホ曲線を描く（Python, Ruby）

自己相似図形であるコッホ曲線を Python と Ruby で描いてみました。

Python では手軽にタートル・グラフィックスが使えるので、これを利用するのが簡単です。

　
3次のコッホ曲線を描きます。

from turtle import *

def draw(length, depth):
    if depth == 0:
        forward(length)
    else:
        draw(length / 3, depth - 1)
        left(60)
        draw(length / 3, depth - 1)
        right(120)
        draw(length / 3, depth - 1)
        left(60)
        draw(length / 3, depth - 1)

color('firebrick')
up()
setx(-250)
down()

draw(500, 3)
input()

　
　
Ruby では簡単なタートル・グラフィックスを実装して描くことにします（Turtle クラス）。描画は自作の Gem 'oekaki' で行っています。

　
4次のコッホ曲線を描きます。

require 'oekaki'

Width, Height = 600, 300

class Turtle
  def initialize(x, y, ob)
    @pen_po = Vector[x, y]
    @pen = ob
    @dir = Vector[1, 0]
  end
  
  def left(deg)
    θ = PI * deg / 180
    @dir = Matrix[[cos(θ), -sin(θ)], [sin(θ), cos(θ)]] * @dir
  end
  
  def right(deg)
    left(-deg)
  end
  
  def forward(length)
    next_po = @pen_po + @dir * length
    @pen.line(Width / 2 + next_po[0], Height / 2 - next_po[1],
       Width / 2 + @pen_po[0], Height / 2 - @pen_po[1])
    @pen_po = next_po
  end
end


Oekaki.app width: Width, height: Height, title: "Koch curve" do
  draw do
    color(0, 0, 0)
    rectangle(true, 0, 0, Width, Height)
    
    t = Turtle.new(-250, -50, self)
    
    drawing = proc do |length, depth|
      if depth.zero?
        t.forward(length)
      else
        drawing[length / 3, depth - 1]
        t.left(60)
        drawing[length / 3, depth - 1]
        t.right(120)
        drawing[length / 3, depth - 1]
        t.left(60)
        drawing[length / 3, depth - 1]
      end
    end
    
    color(0, 65535, 0)
    drawing[500.0, 4]
  end
end

Python でも Ruby でもやっていることはほぼ同じです。

Gem 'oekaki' については以下。
oekaki | RubyGems.org | your community gem host
GTK+でお絵かきしてみた（Ruby） - Camera Obscura
　
誰かえらい人、Ruby でタートル・グラフィックスを実装しないですかね。初心者が遊ぶのに手頃だと思うのだけれど。自分でやれるといいのだけれど、Python みたいにインタラクティブに描画させるようにすると自分のスキルではむずかしそう。
→自分で実装してみた。
　
　
※参考
GTK+でヒルベルト曲線（Ruby） - Camera Obscura
Python の Turtle でヒルベルト曲線 - Camera Obscura

Ruby でローレンツアトラクタを描画する

Ruby + gnuplot でローレンツアトラクタを描いてみました。
　

　
全体的にここなどを参考にしました。微分方程式はオイラー法（参考）で数値計算しています。

gnuplot での描画は numo/gnuplot という Gem を使っています。

require 'numo/gnuplot'

fx = lambda {|x, y, z, r, p, b| -p * x + p * y}
fy = lambda {|x, y, z, r, p, b| -x * z + r * x - y}
fz = lambda {|x, y, z, r, p, b|  x * y - b * z}

dt = 1e-3
x, y, z = 1, 1, 1
p, r, b = 10, 28, 8 / 3.0

ax, ay, az = [], [], []

100000.times do
  x += dt * fx[x, y, z, r, p, b]
  y += dt * fy[x, y, z, r, p, b]
  z += dt * fz[x, y, z, r, p, b]
  ax << x
  ay << y
  az << z
end

Numo.gnuplot do
  unset :key
  splot ax, ay, az, w: :dots
end
gets    #終了待ち

1Ωの抵抗10個で黄金比の値に近づける（Ruby）

問題：

1Ω の抵抗 10個を使い、合成抵抗が黄金比 1.6180339887..Ωにもっとも近づく場合の値を、少数第10位まで求めよ。

　
aΩ と bΩ の抵抗をつなげる場合、直列つなぎにすれば合成抵抗はたんに a + b Ω になりますが、並列つなぎの場合はそれらの逆数の和の逆数、つまり
　　
になるのが重要なところです。抵抗10個できわめて複雑な組み合わせをつくることができます。

Ruby で求めてみました。

require 'set'

def product(a, b)
  @ar[a].each do |i|
    @ar[b].each do |j|
      @ar[a + b] << i + j
      @ar[a + b] << (i * j) / (i + j)
    end
  end
end


@ar = Array.new(11) {Set.new}
@ar[1] << Rational(1, 1)
for i in 2..10
  for j in 1..(i / 2)
    product(j, i - j)
  end
end

ans = 10
@ar[10].each do |x|
  ans = x if (x - 1.6180339887).abs < (ans - 1.6180339887).abs
end
puts "%.10f" % ans.to_f
puts ans

答えは 1.6181818182 となります。分数だと 89/55 ですね。
ちなみに、模範解答よりもずっとシンプルなコードだと思います。

さて、これがどのような回路なのかですが、それも求めるとなるとさらに手を入れないといけないですね。だいぶ複雑になりそうです。
　
　

プログラマ脳を鍛える数学パズル シンプルで高速なコードが書けるようになる70問

ナップサック問題（Ruby）

次のような問題があるとします。

学校でクラブ活動をするのに、選んだクラブの人数の合計が 150人以下になるようにするとします。クラブの想定人数とそれが使う敷地面積は次のように与えられています。

クラブを幾つか適当に選ぶとき、必要な敷地面積の総和の最大値を求めなさい。

　
Ruby で解いてみました。総当り法でも解けますが、動的計画法を使いました。

number = [40, 30, 24, 20, 14, 16, 15, 40, 10, 12]
square = [11000, 8000, 400, 800, 900, 1800, 1000, 7000, 100, 300]

area = Array.new(151, 0)
for i in 0..9
  tmp = area.dup
  area.each_with_index do |sq, j|
    next if sq.zero?
    index = j + number[i]
    next if index > 150
    a = sq + square[i]
    tmp[index] = a if a > tmp[index]
  end
  tmp[number[i]] = square[i] if square[i] > tmp[number[i]]
  area = tmp
end
puts area.max

答えは 28800平方メートルになります。
　

プログラマ脳を鍛える数学パズル シンプルで高速なコードが書けるようになる70問

　
　

追記

配列ではなくハッシュを使うべきかも知れません。この程度ならそれほど実行時間に差はでませんが。

number = [40, 30, 24, 20, 14, 16, 15, 40, 10, 12]
square = [11000, 8000, 400, 800, 900, 1800, 1000, 7000, 100, 300]

area = {}
for i in 0..9
  tmp = area.dup
  area.each do |num, sq|
    index = num + number[i]
    next if index > 150
    a = sq + square[i]
    tmp[index] = a if !tmp[index] or a > tmp[index]
  end
  tmp[number[i]] = square[i] if !tmp[number[i]] or square[i] > tmp[number[i]]
  area = tmp
end
puts area.values.max

別解（2020/11/14）

典型的な動的計画法だとこうなります。

club = [[11000, 40], [8000, 30], [400, 24], [800, 20], [900, 14],
        [1800, 16], [1000, 15], [7000, 40], [100, 10], [300, 12]]
N = 150
C = club.size

area = Array.new(C + 1) {Array.new(N + 1, 0)}

club.each_with_index do |(square, number), i|
  (0..N).each do |num|
    if num >= number
      area[i + 1][num] = [area[i][num - number] + square, area[i][num]].max
    else
      area[i + 1][num] = area[i][num]
    end
  end
end

puts area[C][N]

右折禁止（アルゴリズム・パズル）

アルゴリズム・パズルを Ruby で解いてみました。

　

class TurnLeft
  class Position
    def initialize(x, y)
      @x, @y = x, y
    end
    attr_accessor :x, :y
    
    def +(dir)
      Position.new(@x + dir[0], @y + dir[1])
    end
  end
  
  class Field
    def initialize
      @yoko = Array.new(Height + 1) {Array.new(Width , 0)}
      @tate = Array.new(Width  + 1) {Array.new(Height, 0)}
    end
    
    def set(po, dir)
      if dir[1].zero?
        x = (dir[0] > 0) ? po.x : po.x - 1
        @yoko[po.y][x] = 1
      else
        y = (dir[1] > 0) ? po.y : po.y - 1
        @tate[po.x][y] = 1
      end
    end
    private :set
    
    def get(po, dir)
      if dir[1].zero?
        x = (dir[0] > 0) ? po.x : po.x - 1
        return 1 if x >= Width or x < 0
        a = @yoko[po.y][x]
      else
        y = (dir[1] > 0) ? po.y : po.y - 1
        return 1 if y >= Height or y < 0
        a = @tate[po.x][y]
      end
      set(po, dir) if a.zero?
      a
    end
    
    def dup
      Marshal.load(Marshal.dump(self))
    end
  end
  
  def go
    puts main(Position.new(0, 0), [1, 0], Field.new)
  end
  
  def main(po, dir, field)
    co = 0
    nxt = po + dir
    return 1 if nxt.x == Width and nxt.y == Height   #ゴールか？
    return 0 if nxt.x > Width or nxt.x < 0 or nxt.y > Height or nxt.y < 0
    return 0 unless field.get(po, dir).zero?         #すでに通っているか？
    
    co += main(nxt, dir, field.dup)                  #直進
    co += main(nxt, [-dir[1], dir[0]], field.dup)    #左折
  end
  private :main
end

Width, Height = 6, 4

TurnLeft.new.go

答えは 2760通りです。かかった時間は自分の環境で約 3.7秒でした。

※追記（2018/3/13）
別の方法で解いてみました。

def yoko?(x, y, x1)
  return false if x1 < 0 or x1 > Width
  x2 = [x, x1].min
  return false if @yoko[y][x2] == 1
  x2
end

def tate?(x, y, y1)
  return false if y1 < 0 or y1 > Height
  y2 = [y, y1].min
  return false if @tate[x][y2] == 1
  y2
end

def solve(x, y, dir)
  walk = lambda do |dir|
    dirs = [[1, 0], [0, -1], [-1, 0], [0, 1]]
    x1 = x + dirs[dir][0]
    y1 = y + dirs[dir][1]
    if (dir == 0 or dir == 2) and (x2 = yoko?(x, y, x1))
      @yoko[y][x2] = 1
      solve(x1, y1, dir)
      @yoko[y][x2] = 0
    elsif (dir == 1 or dir == 3) and (y2 = tate?(x, y, y1))
      @tate[x][y2] = 1
      solve(x1, y1, dir)
      @tate[x][y2] = 0
    end
  end
  
  if x == Width and y == 0
    @co += 1
  else
    walk.call(dir)            #直進
    walk.call((dir + 1) % 4)  #左折
  end
end

Width, Height = 6, 4

@yoko = Array.new(Height + 1) {Array.new(Width , 0)}
@tate = Array.new(Width  + 1) {Array.new(Height, 0)}
@co = 0

solve(0, Height, 0)
puts @co    #=>2760

解くのにかかった時間は 0.3秒ほどです。
　

プログラマ脳を鍛える数学パズル シンプルで高速なコードが書けるようになる70問

GTK+ で落書き 8（Ruby）

Gem 'oekaki' で落書きです。
oekaki | RubyGems.org | your community gem host
GTK+でお絵かきしてみた（Ruby） - Camera Obscura

　
再帰的な tree を描いてみました。画像では動きはないですが、じつはアニメーションです。

　

require 'oekaki'

Width, Height = 400, 300
Turn = 20
Length = 40
Ratio = 0.9
Depth = 8

θ = PI * Turn / 180
left  = Matrix[[cos( θ), -sin( θ)], [sin( θ), cos( θ)]]
right = Matrix[[cos(-θ), -sin(-θ)], [sin(-θ), cos(-θ)]]


Oekaki.app width: Width, height: Height, title: "tree" do
  draw do
    color(0, 0, 0)
    rectangle(true, 0, 0, Width, Height)
  end
  
  ar = [[Vector[Width / 2, 1], Vector[0, Length]]]
  branch_num = 0
  
  id = timer(80) do
    ob = ar.shift
    v1 = left  * ob[1]
    v2 = right * ob[1]
    p1 = ob[0] + v1
    p2 = ob[0] + v2
    
    color(0, 65535, 0)
    line(ob[0][0], Height - ob[0][1], p1[0], Height - p1[1])
    line(ob[0][0], Height - ob[0][1], p2[0], Height - p2[1])
    
    ar << [p1, v1 * Ratio]
    ar << [p2, v2 * Ratio]
    
    branch_num += 2
    Gtk.timeout_remove(id) if branch_num >= 2 ** (Depth + 1) - 2
    true
  end
end

いわゆる「スライドパズル（15パズル）」もどきを Ruby で解いてみる

スライドパズルっていうおもちゃがありますよね。4×4 のマス目に空白がひとつあって（残りはコマ）、コマは空白にスライドさせて動かすしかなく、それを特定のパターンに揃えるというものです。ここではルールを少し改変して、特定のマスを別の特定のマスに移動させるのに最短で何手かかるかを考えます。

例えば 3×2 のパズルでいちばん右下のマスが空白であり、いちばん左上のマスにあるコマをいちばん右下へもってくるには、最短で 9手かかります。

●●●
● ○

●●●
●○ 

●● 
●○●

● ●
●○●

●○●
● ●

●○●
 ●●

 ○●
●●●

○ ●
●●●

○● 
●●●

○●●
●● 

9手で終了

ただしこれは、下の方から見ていって下さい（プログラムの関係でそうなっています）。

僕が見たアルゴリズムパズルの本では、10×10 のパズルでやるのと同等な問題が与えてありました（駐車場で車を脱出させるという体裁になっていました）。これは自分の書いたプログラムだと、69手ということになります（調べたパターンの数は 8927個）。計算時間はこれだと（自分の環境で）およそ 41秒と、かなり微妙な数字になりました。瞬殺は無理でしたね。これ以上広いパズルだと、根本から考えなおさないとダメかも知れません。

コードは以下です。幅優先探索で求めています。また、既に出たパターンに逢着した場合はスキップしています。

class Parking
  def initialize
    start = Pattern.new
    start.car.x, start.car.y = 0, 0
    start.space.x, start.space.y = Width - 1, Height - 1
    @@patterns = [start]
  end
  
  def solve
    buffer = @@patterns.dup
    while (po = buffer.shift)
      [po.up, po.down, po.left, po.right].each do |i|
        next unless i
        if i[1]    #true ならば終了
          @@co = 0
          i[0].show
          puts "#{@@co - 1}手で終了"
          puts "調べたパターン数は#{@@patterns.size}"
          exit
        end
        buffer << i[0]
      end
    end
  end
end

class Pattern < Parking
  class Car   < Struct.new(:x, :y); end
  class Space < Struct.new(:x, :y); end
  
  def initialize
    @car   = Car.new
    @space = Space.new
    @parent = nil
  end
  attr_accessor :car, :space, :parent
  
  def up
    x, y = space.x, space.y - 1
    return false if y < 0
    check(x, y)
  end
  
  def down
    x, y = space.x, space.y + 1
    return false if y >= Height
    check(x, y)
  end
  
  def left
    x, y = space.x - 1, space.y
    return false if x < 0
    check(x, y)
  end
  
  def right
    x, y = space.x + 1, space.y
    return false if x >= Width
    check(x, y)
  end
  
  def check(x, y)
    nxt = Pattern.new
    nxt.space.x, nxt.space.y = x, y
    nxt.car = (@car.x == x and @car.y == y) ? @space : @car
    
    @@patterns.each do |pa|
      return false if nxt.car == pa.car and nxt.space == pa.space 
    end
    nxt.parent = self
    
    @@patterns << nxt
    [nxt, (nxt.car.x == Width - 1 and nxt.car.y == Height - 1)]
  end
  private :check
  
  def putout
    Height.times do |i|
      st = "●" * Width
      st[@car.x]   = "○" if @car.y   == i
      st[@space.x] = " " if @space.y == i
      puts st
    end
    puts
  end
  
  def show
    @@co += 1
    putout
    @parent.show if @parent
  end
  
  def inspect
    "<Pattern:car(#{@car.x},#{@car.y}), space(#{@space.x},#{@space.y})>"
  end
end

Width = 4; Height = 4

Parking.new.solve

なお、このプログラムは Pattern#putout で String リテラルの破壊的変更を行っています。
　
ちなみに「15パズル（4×4）」の場合、21手で終了します。時間は約 0.1秒です。

プログラマ脳を鍛える数学パズル シンプルで高速なコードが書けるようになる70問