計算機プログラムの構造と解釈 第二版（第一章）

木構造再帰

#フィボナッチ数列
fib = ->(n) {
  case n
  when 0 then 0
  when 1 then 1
  else
    fib[n - 1] + fib[n - 2]
  end
}

p fib[20]    #=>6765
p fib[100]   #フリーズ！


fib = ->(n) {
  fib_iter = ->(a, b, count) {
    count.zero? ? b : fib_iter[a + b, a, count - 1]
  }
  fib_iter[1, 0, n]
}

p fib[100]    #=>354224848179261915075

#amountセントの金額を、50, 25, 10, 5, 1セントの硬貨で両替するしかたは何通りあるか
count_change = ->(amount) {
  cc = ->(amount, kinds_of_coins) {
    first_denomination = ->{
      case kinds_of_coins
      when 1 then 1
      when 2 then 5
      when 3 then 10
      when 4 then 25
      when 5 then 50
      end
    }
    
    if amount.zero?
      1
    elsif amount < 0 or kinds_of_coins.zero?
      0
    else
      cc[amount, kinds_of_coins - 1] + cc[amount - first_denomination[], kinds_of_coins]
    end
  }
  cc[amount, 5]
}

p count_change[100]    #=>292

問題1.11。

f = ->(n) {
  (n < 3) ? n : f[n - 1] + 2 * f[n - 2] + 3 * f[n - 3]
}

p f[10]    #=>1892
p f[100]   #フリーズ！


f = ->(n) {
  f_iter = ->(a, b, c, n) { n.zero? ? a : f_iter[a + 2 * b + 3 * c, a, b, n - 1] }
  (n < 3) ? n : f_iter[2, 1, 0, n - 2]
}

p f[100]    #=>11937765839880230562825561449279733086

問題1.12。パスカルの三角形。

pascal = ->(n) {
  each_cons = ->(ar) {
    if ar.size < 2
      []
    else
      [ar[0] + ar[1]] + each_cons[ar.drop(1)]
    end
  }
  if n.zero?
    [1]
  else
    [1] + each_cons[pascal[n - 1]] + [1]
  end
}

6.times {|i| p pascal[i]}
#[1]
#[1, 1]
#[1, 2, 1]
#[1, 3, 3, 1]
#[1, 4, 6, 4, 1]
#[1, 5, 10, 10, 5, 1]

べき乗

expt = ->(b, n) {
  n.zero? ? 1 : b * expt[b, n - 1]
}

p expt[2, 5]     #=>32


fast_expt = ->(b, n) {
  is_even = ->(n) { (n % 2).zero? }
  square = ->(x) { x * x }
  
  if n.zero?
    1
  elsif is_even[n]
    square[fast_expt[b, n / 2]]
  else
    b * fast_expt[b, n - 1]
  end
}

p fast_expt[2, 100]    #=>1267650600228229401496703205376

問題1.16。これは思いつかない。ここを見た。

fast_expt_iter = ->(b, n) {
  expt_iter = ->(b, n, a) {
    if n.zero?
      a
    elsif n.even?
      expt_iter[b * b, n / 2, a]
    else
      expt_iter[b, n - 1, a * b]
    end
  }
  expt_iter[b, n, 1]
}

p fast_expt_iter[2, 10]    #=>1267650600228229401496703205376

問題1.17。

fast_mul = ->(a, b) {
  double = ->(n) {n * 2}
  halve = ->(n) {n / 2}
  
  if b.zero?
    0
  elsif b.even?
    double[fast_mul[a, halve[b]]]
  else
    a + fast_mul[a, b - 1]
  end
}

p fast_mul[423568, 12457984]    #=>5276803366912

問題1.18。

fast_mul_iter = ->(a, b) {
  double = ->(n) {n * 2}
  halve = ->(n) {n / 2}
  
  mul_iter = ->(a, b, m) {
    if b.zero?
      m
    elsif b.even?
      mul_iter[double[a], halve[b], m]
    else
      mul_iter[a, b - 1, a + m]
    end
  }
  mul_iter[a, b, 0]
}

p fast_mul_iter[423568, 12457984]    #=>5276803366912

問題1.19。対数的ステップ数のフィボナッチ数列。

fib = ->(n) {
  fib_iter = ->(a, b, p, q, count) {
    if count.zero?
      b
    elsif count.even?
      fib_iter[a, b, p * p + q * q, 2 * p * q + q * q, count / 2]
    else
      fib_iter[b * q + a * q + a * p, b * p + a * q, p, q, count - 1]
    end
  }
  fib_iter[1, 0, 0, 1, n]
}

p fib[1000]
#=>434665576869374564356885276750406258025646605173717804024817290895
#  365554179490518904038798400792551692959225930803226347752096896232
#  398733224711616429964409065331879382989696499285160037044761377951
#  66849228875

これはこういうこと。
　　
とすると
　　
ただし
　　
とする。
　

素数性のテスト

smallest_divisor = ->(n) {
  find_divisor = ->(test_divisor) {
    divides = ->(a, b) {(b % a).zero?}
    
    if test_divisor ** 2 > n
      n
    elsif divides[test_divisor, n]
      test_divisor
    else
      find_divisor[test_divisor + 1]
    end
  }
  find_divisor[2]
}

is_prime = ->(n) { n == smallest_divisor[n] }

2.upto(20) {|i| p [i, is_prime[i]]}

問題1.21。

p smallest_divisor[199]    #=>199
p smallest_divisor[1999]   #=>1999
p smallest_divisor[19999]  #=>7

問題1.22。

timed_prime_test = ->(n) {
  start_prime_test = ->(start_time) {
    report_prime = ->(t) { puts "\t *** %3.6f sec" % t }
    
    report_prime[Time.now - start_time] if is_prime[n]
  }
  print n.to_s
  start_prime_test[Time.now]
}

search_for_primes = ->(a, b) {
  timed_prime_test[a]
  search_for_primes[a + 1, b] unless a == b
}

結果。

1000    ~
     1009 *** 0.000036 sec
     1013 *** 0.000031 sec
     1019 *** 0.000029 sec
10000   ~
    10007 *** 0.000111 sec
    10009 *** 0.000084 sec
    10037 *** 0.000106 sec
100000  ~
   100003 *** 0.000568 sec
   100019 *** 0.000288 sec
   100043 *** 0.000227 sec
1000000 ~
  1000003 *** 0.002030 sec
  1000033 *** 0.000813 sec
  1000037 *** 0.000784 sec

それぞれ平均値で
  3.135416666666667
  3.5980066445182723
  3.3490304709141276
倍で、√10 = 3.1622776601683795 倍に近い。

問題1.23。

smallest_divisor = ->(n) {
  nxt = ->(i) {(i == 2) ? 3 : i + 2}
  divides = ->(a, b) {(b % a).zero?}
  
  find_divisor = ->(test_divisor) {
    if test_divisor ** 2 > n
      n
    elsif divides[test_divisor, n]
      test_divisor
    else
      find_divisor[nxt[test_divisor]]
    end
  }
  find_divisor[2]
}

check_primes = ->(ar) {
  return if ar.empty?
  timed_prime_test[ar.first]
  check_primes[ar.drop(1)]
}

ar = [1009, 1013, 1019, 10007, 10009, 10037,
      100003, 100019, 100043, 1000003, 1000033, 1000037]
check_primes[ar]
puts

結果。

1009	 *** 0.000021 sec
1013	 *** 0.000011 sec
1019	 *** 0.000017 sec
10007	 *** 0.000057 sec
10009	 *** 0.000027 sec
10037	 *** 0.000024 sec
100003	 *** 0.000174 sec
100019	 *** 0.000088 sec
100043	 *** 0.000087 sec
1000003	 *** 0.000536 sec
1000033	 *** 0.000282 sec
1000037	 *** 0.000278 sec

Fermatテスト。

expmod = ->(base, exp, m) {
  if exp.zero?
    1
  elsif exp.even?
    expmod[base, exp / 2, m] ** 2 % m
  else
    base * expmod[base, exp - 1, m] % m
  end
}

fermat_test = ->(n) {
  try_it = ->(a) {expmod[a, n, n] == a}
  try_it[rand(1..n - 1)]
}

is_fast_prime = ->(n, times) {
  if times.zero?
    true
  elsif fermat_test[n]
    is_fast_prime[n, times - 1]
  else
    false
  end
}

ar = [1009, 1013, 1019, 10007, 10009, 10037,
      100003, 100019, 100043, 1000003, 1000033, 1000037]
ar.each {|x| p [x, is_fast_prime[x, 1]]}

問題1.24。

timed_prime_test = ->(n) {
  start_prime_test = ->(start_time) {
    report_prime = ->(t) { puts "\t *** %3.6f sec" % t }
    
    report_prime[Time.now - start_time] if is_fast_prime[n, 3]
  }
  print n.to_s
  start_prime_test[Time.now]
}

ar.each {|x| timed_prime_test[x]}

結果。

1009	 *** 0.000040 sec
1013	 *** 0.000024 sec
1019	 *** 0.000020 sec
10007	 *** 0.000025 sec
10009	 *** 0.000031 sec
10037	 *** 0.000021 sec
100003	 *** 0.000033 sec
100019	 *** 0.000026 sec
100043	 *** 0.000024 sec
1000003	 *** 0.000031 sec
1000033	 *** 0.000027 sec
1000037	 *** 0.000031 sec

問題1.27。

carmichael = [561, 1105, 1729, 2465, 2821, 6601]
carmichael.each {|x| p [x, is_fast_prime[x, 3], is_prime[x]]}
#=>
#[561 , true, false]
#[1105, true, false]
#[1729, true, false]
#[2465, true, false]
#[2821, true, false]
#[6601, true, false]

引数としての手続き

#数列の和
sum = ->(term, a, nxt, b) {
  if a > b
    0
  else
    term[a] + sum[term, nxt[a], nxt, b]
  end
}

inc = ->(n) {n + 1}
cube = ->(n) {n ** 3}

sum_cubes = ->(a, b) {
  sum[cube, a, inc, b]
}

p sum_cubes[1, 10]    #=>3025

identity = ->(x) {x}

sum_integers = ->(a, b) {
  sum[identity, a, inc, b]
}

p sum_integers[1, 10]   #=>55

pi_sum = ->(a, b) {
  pi_term = ->(x) {1.0 / (x + 2) / x}
  pi_next = ->(x) {x + 4}
  
  sum[pi_term, a, pi_next, b]
}

p 8 * pi_sum[1, 1000]    #=>3.139592655589783

#積分
integral = ->(f, a, b, dx) {
  add_dx = ->(x) {x + dx}
  
  dx * sum[f, dx / 2.0 + a, add_dx, b]
}

p integral[cube, 0, 1, 0.01 ]    #=>0.24998750000000042
p integral[cube, 0, 1, 0.001]    #=>0.249999875000001
#正確な値は 0.25

問題1.29。Simpson の公式。

simpson = ->(f, a, b, n) {
  coef = ->(i) {
    if i == 0 or i == n
      1
    else
      i.even? ? 2 : 4
    end
  }
  h = (b - a) / n.to_f
  y = ->(k) {coef[k] * f[a + k * h]}
  sum[y, 0, inc, n] * h / 3
}

p simpson[cube, 0, 1, 100 ]    #=>0.24999999999999992
p simpson[cube, 0, 1, 1000]    #=>0.2500000000000002

問題1.30。反復的に。

sum = ->(term, a, nxt, b) {
  iter = ->(i, result) {
    if i > b
      result
    else
      iter[nxt[i], result + term[i]]
    end
  }
  iter[a, 0]
}

p sum_integers[1, 10]    #=>55

問題1.31。

#再帰的
product = ->(term, a, nxt, b) {
  if a > b
    1
  else
    term[a] * product[term, nxt[a], nxt, b]
  end
}

factorial = ->(n) {
  product[identity, 1, inc, n]
}

p factorial[10]    #=>3628800


square = ->(n) {n ** 2}

pi = ->(n) {
  term = ->(x) {square[x].to_f / square[x + 1]}
  pi_next = ->(x) {x + 2}
  4 * product[term, 2, pi_next, n - 2] * n / 2
}

p pi[5000]    #=>3.141906828561547

#反復的
product = ->(term, a, nxt, b) {
  iter = ->(i, result) {
    if i > b
      result
    else
      iter[nxt[i], result * term[i]]
    end
  }
  iter[a, 1]
}

p pi[5000]    #=>3.14190682856154

問題1.32。

#再帰的
accumulate = ->(combiner, null_value, term, a, nxt, b) {
  accum = ->(i) {
    if i > b
      null_value
    else
      combiner[term[i], accum[nxt[i]]]
    end
  }
  accum[a]
}

add = ->(x, y) {x + y}

sum = ->(term, a, nxt, b) {
  accumulate[add, 0, term, a, nxt, b]
}

sum_integers = ->(a, b) {
  sum[identity, a, inc, b]
}

p sum_integers[1, 10]    #=>55

#反復的
accumulate = ->(combiner, null_value, term, a, nxt, b) {
  iter = ->(i, result) {
    if i > b
      result
    else
      iter[nxt[i], combiner[result, term[i]]]
    end
  }
  iter[a, null_value]
}

mul = ->(x, y) {x * y}

product = ->(term, a, nxt, b) {
  accumulate[mul, 1, term, a, nxt, b]
}

p factorial[10]    #=>3628800

問題1.33。

filtered_accumulate = ->(combiner, null_value, term, a, nxt, b, filter) {
  filtered = ->(x) {filter[x] ? x : filtered[nxt[x]]}
  iter = ->(i, result) {
    i = filtered[i]
    if i > b
      result
    else
      iter[nxt[i], combiner[result, term[i]]]
    end
  }
  iter[a, null_value]
}

#区間a, bの素数の二乗の和
problem133a = ->(a, b) {
  filtered_accumulate[add, 0, square, a, inc, b, is_prime]
}

p problem133a[2, 10]    #=>87


gcd = ->(a, b) { b.zero? ? a : gcd[b, a % b] }

#nと互いに素で, nより小さい正の整数の積
problem133b = ->(n) {
  filter = ->(x) {gcd[x, n] == 1}
  filtered_accumulate[mul, 1, identity, 1, inc, n - 1, filter]
}

p problem133b[10]    #=>189

一般的方法としての手続き

区間二分法。

average = ->(a, b) {(a + b) / 2.0}

is_colse_enough = ->(x, y) {(x - y).abs < 0.001}

search = ->(f, neg, pos) {
  midpoint = average[neg, pos]
  if is_colse_enough[neg, pos]
    midpoint
  else
    test_value = f[midpoint]
    if test_value > 0
      search[f, neg, midpoint]
    elsif test_value < 0
      search[f, midpoint, pos]
    else
      midpoint
    end
  end
}

half_interval_method = ->(f, a, b) {
  a_value = f[a]
  b_value = f[b]
  if a_value < 0 and b_value > 0
    search[f, a, b]
  elsif b_value < 0 and a_value > 0
    search[f, b, a]
  else
    raise "Values are not of opposite sign #{a} #{b}"
  end
}

sin = Math.method(:sin)
p half_interval_method[sin, 2.0, 4.0]    #=>3.14111328125

p half_interval_method[->(x) {x ** 3 - 2 * x - 3}, 1.0, 2.0]
#=>1.89306640625

関数の不動点の探索。

tolerance = 0.00001

fixed_point = ->(f, first_guess) {
  is_close_enough = ->(v1, v2) {(v1 - v2).abs < tolerance}
  try = ->(guess) {
    nxt = f[guess]
    if is_close_enough[guess, nxt]
      nxt
    else
      try[nxt]
    end
  }
  try[first_guess]
}

cos = Math.method(:cos)
p fixed_point[cos, 1.0]    #=>0.7390822985224024

p fixed_point[->(y) {sin[y] + cos[y]}, 1.0]    #=>1.2587315962971173

#平均緩和法
sqrt = ->(x) {
  fixed_point[->(y) {average[y, x / y]}, 1.0]
}

p sqrt[2]    #=>1.4142135623746899

問題1.35。黄金比。

p fixed_point[->(x) {1 + 1 / x}, 2.0]    #=>1.6180327868852458

問題1.36。

fixed_point = ->(f, first_guess) {
  is_close_enough = ->(v1, v2) {(v1 - v2).abs < tolerance}
  try = ->(guess) {
    p nxt = f[guess]
    if is_close_enough[guess, nxt]
      nxt
    else
      try[nxt]
    end
  }
  try[first_guess]
}

log = Math.method(:log)
p fixed_point[->(x) {log[1000] / log[x]}, 10.0]
#=>
#2.9999999999999996
#6.2877098228681545
#3.7570797902002955
#5.218748919675316
#4.1807977460633134
#4.828902657081293
#4.386936895811029
#4.671722808746095
#4.481109436117821
#4.605567315585735
#4.522955348093164
#4.577201597629606
#4.541325786357399
#4.564940905198754
#4.549347961475409
#...
#4.555532257016376
#ステップ数は34

#平均緩和法
p fixed_point[->(x) {average[x, log[1000] / log[x]]}, 10.0]
#=>
#6.5
#5.095215099176933
#4.668760681281611
#4.57585730576714
#4.559030116711325
#4.55613168520593
#4.555637206157649
#4.55555298754564
#4.555538647701617
#4.555536206185039
#4.555536206185039
#ステップ数は11

問題1.37。連分数。

#再帰的
cont_frac = ->(n, d, k) {
  fcf = ->(i) {
    if i == k
      n[k] / d[k]
    else
      n[i] / (d[i] + fcf[i + 1])
    end
  }
  fcf[1]
}

#反復的
cont_frac = ->(n, d, k) {
  iter = ->(i, result) {
    if i.zero?
      result
    else
      iter[i - 1, n[i] / (d[i] + result)]
    end
  }
  iter[k - 1, n[k] / d[k]]
}

#黄金比
check = ->(k) {
  a = 1 / cont_frac[->(i) {1.0}, ->(i) {1.0}, k]
  if (a - 1.6180327868852458).abs < 0.00001
    puts k
    a
  else
    check[k + 1]
  end
}

p check[1]
#=>
#13
#1.6180257510729614

問題1.38。自然対数の底。

e_2 = ->(j) {
  n = ->(i) {1.0}
  d = ->(i) {
    if (i - 1) % 3 == 1
      2 * ((i - 1) / 3 + 1)
    else
      1
    end
  }
  cont_frac[n, d, j]
}

puts
p e_2[10]  + 2    #=>2.7182817182817183
p e_2[100] + 2    #=>2.7182818284590455

問題1.39。タンジェント（正接）。

tan_cf = ->(x, k) {
  n = ->(i) {(i == 1) ? x : -x ** 2}
  d = ->(i) {2 * i - 1}
  cont_frac[n, d, k]
}

p tan_cf[Math::PI / 3, 10]    #=>1.732050807568877