『情報科学入門 Rubyを使って学ぶ』の練習問題の解答例の他、授業スライドの練習問題の解答例や、プログラミングのポイント等を掲載します。他クラスの方の利用歓迎。改良・訂正箇所ありましたらご指摘ください。

2010/12/22

[教]練習問題(第六章)

科書「第六章・数値計算」の練習問題の解答例です。


数学が終わっている私は、泣きそうになりながら解いてます。


練習6.1 (台形公式の正確さ)





練習6.2 (Simpson公式による積分)
a) def  simpson(xs,xe,n)
       deltax=(xe-xs)/(2.0*n)
       sum=f(xs)+4*f(xs+deltax)+f(xe)
       for  i  in 1..(n-1)
         sum=sum+2*(xs+2*i*deltax)+4*f(xs+(2*i+1)*deltax)
       end
       deltax*sum/3.0
     end


b)


練習6.3 (Monte Carlo法の正確さ)

練習6.4 (乱数列となる条件) ※ここの解答は今までで一番適当です
a) 地域によって、夏や冬などある程度規則があるので、言えない、と思う><

b) 言える。言えなきゃ宝くじにならないぜ!

c) 各数字が4枚までと決まっているので、乱数とは言えない。


練習6.5 (浮動小数点数の表現)
a) (2)は、二進数であることを示す
 符号 s = 1
 指数 d'1..d'8 = 11111110
 仮数 d1..d23 = 11000000000000000000000
また、単精度のバイアス値 b = 127 より
 (-1)^s × 1.d1..d23(2) × 2^(d'1..d'8(2) - b)
 = (-1)^1 × 1.11000000000000000000000(2) × 2^(11111110(2) - 127)
 = (-1) × { 2^0 + 2^(-1) + 2^(-2) } × 2^(254 - 127)
 = -2.977471 × 10^38

b) (a) s = 0 , d'1..d'8 = 01111111, d1..d23 = 00000000000000000000000

     (b) 3145728 = 2^20 × 3 = (-1)^0 × 1.5 × 2^21 = (-1)^0 × 1.5 × 2^(148 - 127) より
           s = 0 , d'1..d'8 = 10010100 , d1..d23 = 10000000000000000000000

     (c) 5/65536 = 5 / 2^16 = (-1)^0 × 1.25 × 2^(-14) = (-1)^0 × 1.25 × 2^(113 - 127) より
           s = 0 , d'1..d'8 = 01110001, d1..d23 = 01000000000000000000000


練習6.6 (不動小数点の大小)
えっ
符号のとこ見るといっぱつで反証になっちゃうんだけど

x < 0 < y のとき
sx = 1 かつ sy = 0 であるから
(xの浮動小数点数) > (yの不動小数点数)

・・・あれー?
一応符号を無視してやってみます。
でも、私数学苦手ですから。もっとすっきり証明もとむ。


十進数で表されたx,y(0<x<y)はそれぞれ
 x = p × 2^u
 y = q × 2^v
とできる。ただし、
 1 ≦ p,q < 2 かつ -127 ≦ u,v ≦ 128 ー※
である。
I) u=vと仮定すると(つまり指数部が一致すると仮定)
 x<yより、
  p<q ⇔ p-1<q-1 ⇔ dx1..dx23(2)<dy1..dy23(2)
 であるから、
 (xの浮動小数点数)<(yの浮動小数点数)
II) u>vと仮定すると
  y = q × 2^v < 2^(v+1) ≦ 2^u < q × 2^v = x (∵※)
 よって矛盾。
III) u<vと仮定すると
  dx'1..dx'8(2) < dy'1..dy'8(2)
 であるから、
 (xの浮動小数点数)<(yの浮動小数点数)
以上より、証明できた。




練習6.7 (丸め誤差)
0.1*5-0.5 = 0.0
となる。これは、  ※(36)は36桁分の省略があることを示す
 0.1(10) = 1.100110011..(36)..0011010(2) × 2^(-4)(10)
  5(10) = 1.010000000..(36)..0000000(2) × 2^2(10)
 0.5(10) = 1.000000000..(36)..0000000(2) × 2^(-1)(10)
であることを考えると、
 0.1(10) × 5(10) = 10.00000000000..(36)..0000010(2) × 2^(-2)(10)
                         桁落ち
          = 1.000000000..(36)..0000000(2) × 2^(-1)(10)
となり、0.5(10)と一致することから成り立つ。




練習6.8 (倍精度表現の精度)
  仮数部のビット数:52
10進数有効桁数 (約):log10{2^(52+1)} ≒ 16
  指数部のビット数:11
     最小の指数:2^(-1022) ≒ 6.6 × 10^309
     最大の指数:2^(1023) ≒ 8.9 × 10^307
       最大値:2 × 8.9 × 10^307 ≒ 1.8 × 10^308




練習6.9 (Pivotingの完成)
load("./abs.rb")
def  maxrow(a,k)
  s=k
  t=abs(a[k][k])
  if  k!=a.length()-1
    for  y  in k..(a.length()-2)
      if  t < abs(a[y+1][k])
        s=y+1
        t=abs(a[y+1][k])
      end
    end
  end
  s
end

gjp(a)は教科書P137にあるので省略しますよ。


練習6.10 (部分の通過距離)
load("./abs.rb")
include(Math)
def  penetration(theta,d,x,y)
  q = abs(x*cos(theta) + y*sin(theta) - d)
  if  q<0.5
    sqrt(1-4*q**2)
  else
    0.0
  end
end


練習6.11 (係数行列の作成)
def  coefficients(l,n)
  m = make2d(l*n,l**2)
  for  k  in 0..(n-1)
    for  d  in  0..(l-1)
      for  y  in  0..(l-1)
        for  x  in  0..(l-1)
          m[l*k+d][l*y+x] = penetration(3.1415926535897932384*k/n,d,x,y)
        end
      end
    end
  end
  m
end


練習6.12 (係数行列と定数ベクトルの合成)
def  make_equations(l,n,m,s)
  a=make2d(l*n,l**2+1)
  for j in 0..(l*n-1)
    for i in 0..(l**2-1)
      a[j][i]=m[j][i]
    end
    a[j][l**2]=s[j/l][j%l]
  end
  a
end


練習6.13 (画像の復元)
def  equations_to_image(l,solved)
  a=make2d(l,l)
  for j in 0..(l-1)
    for i in 0..(l-1)
      a[j][i]=solved[l*j+i][l**2]
    end
  end
  a
end





練習6.14 (解析的な積分が難しい関数の数値積分)
まずはg(x)を定義。ちなみにsinやlogが出てくるので、数学関数を含むことを宣言しましょう。
Include(Math)
def g(x)
  sin(x)*1.0/log(x)
end

a) def trapezoid_sinlog(xs,xe,n)
      deltax = (xe-xs)*1.0/n
      sum = (g(xs)+g(xe))/2.0
      for i in 1..(n-1)
        sum = sum + g(xs+i*deltax)
      end
      deltax * sum
    end

b) def simpson_sinlog(xs,xe,n)
      deltax = (xe-xs)/(2.0*n)
      sum = g(xs)+4*g(xs+deltax)+g(xe)
      for i in 0..(n-1)
        sum = sum + 2*g(xs+2*i*deltax)+4*g(xs+(2*i+1)*deltax)
      end
      deltax*sum/3.0
    end


練習6.15 (Monte Carlo法による3次元の積分)
include(Math)   ←乱数rand()は数学関数に含まれる
def  montecarlo3d(n)

  m=0
    for  i  in 1..n
      x=rand()
      y=rand()
      z=rand()
      if  x*x+y*y+z*z<1.0
        m=m+1
      end
    end
  m*8.0/n
end
誤差は自分で比べてみましょう。


練習6.16 (Monte Carlo法による積分)
include(Math) 
def f(xs,x)
  (x+xs)*1.0/((x+xs+1)*(x+xs+2))
end

def montecarlo_f(xs,xe,n)
  m=0
    for i in 1..n
      x=rand()*(xe-xs)
      y=rand()*(3-2*sqrt(2))
      if y<f(xs,x)
        m=m+1
      end
    end
  m*(xe-xs)*(3-2*sqrt(2))/n
end

2010/12/15

[ス]練習問題(第六章)

授業スライド「06数値計算」の練習問題の解答例です。 

えーーーーーーー
皆さん今日は。S23-7 情報科学シケタイです。
突然ですが私数学がとっても苦手でして、この章に突入してからあわあわな感じです。
それを考慮の上ご覧ください><




  • 以下の積分(省略w)を台形公式で近似して円周率を求めよ。n を引数として円周率を返す関数を定義せよ。
trapezoid.rb の f を定義しなおし、
def pi(n)
4.0 * trapezoid(....)
end

def  f(x)
  1.0/(1+x**2)
end


def  pi(n)
  4.0*trapezoid(0,1,n)
end




  •  以下の積分をシンプソン公式で近似して円
    周率を求めよ。n を引数として円周率を返
    す関数を定義せよ
def  f(x)
  4.0/(1+x**2)
end


def  simpson(xs,xe,n)
  deltax=(xe-xs)/(2.0*n)
  sum=f(xs)+4*f(xs+deltax)+f(xe)
  for  i  in 1..(n-1)
    sum=sum+2*(xs+2*i*deltax)+4*f(xs+(2*i+1)*deltax)
  end
  deltax*sum/3.0
end


  • 4*montecarlo(n) を色々な引数で実行し、円周率が近似できることを確かめよ。
  • show(mcplot(n)) を実行して、点がばらまかれる様子を観察せよ。– 配列の引数に浮動小数点数を与えると切り捨てられる。

2010/12/03

[管理情報] 教科書第四章の解答について

教科書第四章の練習問題の解答で重大なミスが発見されました。

練習4.6、練習4.7で、「繰り返しを用いて」と指定されているにも関わらず、
再帰を用いた解答を掲載していました。
今現在訂正を急いでおります。


12/04 0:12 練習4.6を訂正いたしました。
12/04 0:47 練習4.7を訂正いたしました。

2010/12/01

[基]ターミナルへの命令

Rubyを起動する前、ターミナルへの命令をいくつかあげます。

pwd    今プログラムが働いている階層(フォルダ)を確認する
cd ruby  今いる階層内のフォルダ"ruby"に移動する
cd ~    最初にプログラムが働いていた階層に戻る
is     今いる階層内のファイル名を表示する
cat bmi.rb 今いる階層内のファイル"bmi.rb"の内容を表示する
覚えましょう。

2010/11/29

[ス]練習問題(第四章)

授業スライド「04関数から計算へ」の練習問題の解答例です。


  • xがyで割り切れることを判定する関数 divisible(x,y) を定義せよ。

def  divisible(x,y)
  x%y==0
end



  •  素数とは、1 と自分自身しか約数がないよう な数である。上で定義した関数 sod を使って n が素数のときにのみ true, そうでないときに false となるような関数 prime(n) を定義せよ。

def  prime(n)
  sod(n,n-1) ==1
end



  •  n 個から k 個を選ぶ組み合わせ数 nCk を求め る combination(n,k) を定義せよ。

def  combination(n,k)
  if  k>n
    0
  else
    if  k==0
      1
    else
      combination(n-1,k-1) + combination(n-1,k)
    end
  end
end



  •  関数 tnpo(n) は n が偶数なら 1/2, 奇数なら 3 倍して 1 加えた数を求めるものだった。数学 者 Collatz はどんな整数 n が与えられたときでも、この関数を使って数を変化させてゆくと、 いずれ1 になると予想した。例えば 3 から始めた場合は3⇒10⇒5⇒16⇒8⇒4⇒2⇒1 と いった具合に予想通りになっていることが確められる。そこで n から上の手順で数を変化 させて 1 になるまでの回数を collatz(n) とする。 例えば collatz(5) = 5, collatz(16) = 4 である。
A)    collatz(n) と collatz(tnpo(n)) の関係を書け。
collatz(n)=collatz(tnpo(n))+1ですね。
B)    collatz(n) を求める関数 collatz(n) を定義せよ。 

def  collatz(n)


  if   n==1


    0


  else
    collatz(tenpo(n))+1
  end
end




  • Sierpinski のカーペット   
    –n 次のカーペットは、縦横が3**n で、n-1次ののカーペットを8 枚敷き詰めて作られる。真ん中は空いている。0次のカーペットは、縦横1 の黒い正方形とする。(実際のプログラムでは、白黒が反転する。)

def  sub(a,n,x,y)
  if  n==0
    a[y][x]=1
  else
    for  i  in 0..2
      for  j  in 0..2
        if  !(i==1&&j==1)
          sub(a,n-1,x+j*3**(n-1),y+i*3**(n-1))
        end
      end
    end
  end
end


def  cantor_dust(n)
  a=make2d(3**n,3**n)
  sub(a,n,0,0)
  a
end


  • Sierpinskiの三角形  –大きさn*n で、i 行j 列目がiCjを2 で割った余りになっているような配列を作る関数sierpinski(n) を定義せよ。(見易さのために白黒を逆にしている。)

def  sierpinski(n)
  a=make2d(n,n)
  a[0][0]=1
  for  i  in  1..(n-1)
    a[i][0]=1
    for  j  in  1..(n-1)
      a[i][j]=(a[i-1][j-1]+a[i-1][j])%2
    end
  end
  a
end


  • match の定義ではs 中にp が必ず現われることを仮定していた。p が現われない場合に-1 と答えるmatch_safe(s,p) を定義せよ。

def  match_safe(s,p)
  i=0
  w=p.length()
  while submatch(s,i,p,w)<w && i<=s.length()-w
    i=i+1
  end
  if  i>s.length()-w
    -1
  else
    i
  end
end

2010/11/17

[基]数列を作る(第二章)

1次元数列を作る
試しに,a0=1,a1=2,a2=3,a3=4,a4=5という数列を作ってみましょう。
a = [1,2,3,4,5]
はい終わり。


以下簡単に
数列を参照する→ a
数列のある項を参照する→ a[2]↓ (これでa2を参照できる。3と表示されるはず)
数列のある項を書き換える→ a[2] = 10↓ (これでa2が10になる)


2次元数列(画像)を作る
Rubyでお絵描きをしてみましょう。まずはirbを終了して、isrbを起動します。
isrbの起動は、"isrb"と入力してReturnを押すだけ。
さて、isrbでは各項の要素を0〜1の数値にすることで、数列を黒〜白のモノトーン画像として出力することができます。
試しに作ってみましょうか。

b = [ [0,0,0],
         [1,1,1],
         [0,0,0] ]

数列の中に数列を突っ込んでるんですね。これで、二次元数列bが出来ました。
これを画像として出力してみましょう。isrb独自の関数showを使います。

show(b)

こうなりましたか?
これで0が黒になることがわかりましたね。
0.1や0.6なんかを要素にして画像出力してみましょう。グレーになりますよ。





二次元数列に対する操作は、ほとんど一次元数列と同じなのですが、二次元なので項数の表現が違います。
b[1][2]だったら、b数列の1行目の2列目の項、ということになります。(通常b2,1など横の座標→縦の座標と表現するものですが、Ruby上では逆です。)
そして注意してほしい点は、この行と列の数え方は、どちらも0から始まっているということ。つまり一番左上の項はb[0][0]になります。
そしてさきほどのb[1][2]は、真ん中の列の一番右の数字のことになります。

以下簡単に

数列を参照する→ b
数列のある項を参照する→ b[1][2]↓ (これでb2,1を参照できる。1と表示されるはず)
数列のある項を書き換える→ b[1][2] = 0.5↓ (これでb2,1が0.5になる)


3次元数列(カラー画像)を作る
さて、今度はカラー画像を作ってみましょう。
二次元数列で、0~1の数字を入れていたところに、さらにRGB数列を入れて、色を表現します。
RGB数列は以下の様になっています。
 [(Red値),(Green値),(Blue値)]
それぞれのカラー値は、0~1で表します。
主な色のRGB数列をのせます。
[1,0,0]…赤
[0,1,0]…緑
[0,0,1]…青
[1,1,0]…黄色
[0,1,1]…水色
[1,0,1]…紫
[0,0,0]…黒
[1,1,1]…白
少数を使って、自分の出したい色をうまく表現してみましょう!

では肝心のこのRGB数列の使い方を。
c = [ [ [1,0,0], [0.5,0.5,0], [0,1,0] ],
        [ [0.5,0,0.5], [0.5,0.5,0.5], [0,0.5,0] ],
        [ [0,0,1], [0,0,0.5], [0,0,0] ] ]
まあこんな感じで、数列の中に数列を突っ込んで、さらにRGB数列を突っ込むんですね。
これで3×3のカラー画像ができたので、
show(c)
で表示してみましょう。

以下簡単に
数列を参照する→ c
数列のある項を参照する→ c[1][2]↓ (これでc2,1を参照できる。[0,0.5,0]と表示されるはず)
数列のある項を書き換える→ c[1][2] = [0,1,1]↓ (これでc2,1が[0,1,1]になる)

以上でっす。

[基]関数を定義する(第一章)

rubyで関数を定義する方法を簡単に説明します。

次のような関数を定義したいとします。

2点(x,y),(u,v)の距離を求める関数
distance(x,y,u,v)={(x-u)^2+(y-v)^2}^(1/2)

distance(x,y,u,v)は,f(x)みたいなものですね。

ではさっそく入力に入ります。この文字が入力部分です。
(以下,returnキーを押すことを↓と示します)

まずは,数学関数(ここではルート)を使うことを宣言しましょう。


include(Math)

次の行に"Object"と表示されたら完了。
関数定義に入りましょう。

def distance(x,y,u,v)↓   ←"def"は"definition"の事
  sqrt((x-u)**2+(y-v)**2)↓  ←関数の内容(sqrt(〜)がルートだったのは覚えてる?)
end            ←定義終わり

次の行に"nil"って表示されたら完了です☆ミ

次に,定義した関数を使ってみましょう。

原点(0,0)と点(3,4)の距離を調べましょう。
入力は
distance(0,0,3,4)
だけ^^

さあ次の行に5って表示されたら成功。
違っていたらどこかでミスしてます。


☆まとめ☆

関数の定義の仕方
def [関数名]
 [関数式]
end

2010/11/14

[教]練習問題(第四章)

教科書「第四章・関数から計算へ」の練習問題の解答例です。




練習4.1 (素数の判定)
def  prime(n)
  sod(n,n-1) ==1   ←sod(n,k)...1〜kのうちnの約数である数の和を求める関数

end
この場合、n=1のときtrueとなる


別解
def  prime(n)
  sod(n,n) == 1+n
end
この場合、n=1のときfalseとなる


練習4.2 (組み合わせ数)
def  combination(n,k)
  if  k>n
    0
  else
    if  k==0
      1
    else
      combination(n-1,k-1) + combination(n-1,k)
    end
  end
end




練習4.3 (素数の判定の改良)
def  prime2(n)
  n==1 ||  gd(n,n-1)==1
end




練習4.4 (再起と繰り返しの比較)




練習4.5 (RNA塩基配列の探索)
実行してみるとわかりますが、このアルゴリズムにはすごい無駄があります。
でもこれ以上がみつからないので、一応掲載しておきます。
def  rna(s)
  w=s.length()
  a=Array.new(w/3)
  i=0
  j=0
  while i<=w-3
    if s[i..(i+2)]=="AUG"
      a[j]=i
      j=j+1
    end
    i=i+1
  end
  a
end




練習4.6 (繰り返しによって値を集める)
a) def  factorial_loop(n)
     f=1
     for  i  in 1..n
       f=f*i
     end
     f
   end


b) def  power2_loop(n)
     p=1
     for  i  in 1..n
       p=p*2
     end
     p
   end

c) def  power_loop(x,n)
     p=1
     for  i  in 1..n
       p=p*x
     end
     p
   end


d) def  taylor_e_loop(x,n)
     t=1
     for  i  in 1..n
       t=t+power_loop(x,i)/factorial_loop(i)
     end
     t
   end



練習4.7 (繰り返しによって条件を満たす値を求める)
a) def  nod_loop(k,n)
     d=0
     for  i  in 1..n
       if  k%i==0
         d=d+1
       end
     end
     d
   end


b) def  nop_loop(n)
     p=0
     for  i  in 1..n
       if  prime2(n)
         p=p+1
       end
     end
     p
   end

c) def  msod_loop(n)
     m=1
     for  i  in 2..n
       if  sod(i,i)>=m
         m=sod(i,i)
       end
     end
     m
   end



練習4.8 (繰り返しによって条件を満たす値を探す)
a) load("./prime.rb")    ←素数を判定する関数(練習4.1で定義)
   def np_loop(n)
     while  !prime2(n)
       n=n+1
     end
     n
   end


b) load("./prime.rb")
   def nth_prime_loop(p,n)
     j=0           ←jはpから何番目に大きい素数かを表す局所関数
     while  j<n
       p=np_loop(p+1)   ←pを「元のpの次に大きい素数」に書き換える
        j=j+1
     end
     p
   end


c) def  collatz_loop(n)
     j=0
     while n!=1
       n=tnpo(n)
       j=j+1
     end
     j
   end

d) def  perfect_loop(n)
     sod(n,n-1)==n
   end
   def  next_perfect_loop(n)
     while  !perfect_loop(n)
       n=n+1
     end
     n
   end


練習4.9 (再帰的に値を集める)
a) def  factorial(n)
     if  n==1
       1
     else
       factrial(n-1)*n
     end
   end

b) def  power2(n)
     if  n==0
       1
     else
       power2(n-1)*2
     end
   end

c) def  power(x,n)
     if  n==0
       1
     else
       power(x,n-1)*x
     end
   end

d) def  taylor_e(x,n)
     if  n==0
       1
     else
       taylor_e(x,n-1)+power(x,n)/factorial(n)
     end
   end


練習4.10 (再帰的に条件を満たす値を集める)
a) def  nod(k,n)
     if  n==1
       1
     else
       if  k%n==0
         nod(k,n-1)+1
       else
         nod(k,n-1)
       end
     end
   end


b) def  nop(n)
     if  n==1
       1
     else
       if  prime2(n)
         nop(n-1)+1
       else
         nop(n-1)
       end
     end
   end


c) def  msod(n)
     if  n==1
       1
     else
       if  sod(n,n)<=msod(,n-1)
         msod(n-1)
       else
         sod(n,n)
       end
     end
   end




練習4.11 (再帰的に条件を満たす値を探す)
a) def  np(n)
     if  prime2(n)
       n
     else
       np(n+1)
     end
   end


b) def nth_prime(p,n)
     if  n==1
       np(p+1)
     else
       nth_prime(np(p+1),n-1)
     end
   end


c) def  collatz(n)
     if   n==1
       0
     else
       collatz(tenpo(n))+1
     end
   end


d) def  perfect_loop(n)
     sod(n,n-1)==n
   end
   def  next_perfect_loop(n)
     if  perfect_loop(n)
       n
     else
       next_perfect_loop(n+1)
     end
   end


練習4.12 (配列を作る)
a) def  make3d(a,b,c)
     v=Array.new(a)
  for  i  in 0..(a-1)
    v[i]=make2d(b,c)
  end
end

b) def  makend(n,m)
     a=Array.new(m)
     if  n==1
       for  i  in 0..(m-1)
         a[i]=0
       end
     else
       for  i  in 0..(m-1)
         a[i]=makend(n-1,m)
       end
     end
     a
   end

c) def sub(d,i)
     if  i==0
       a=make1d(d[0])
     else
       a=Array.new(d[i])
       for j in 0..(d[i]-1)
         a[j]=sub(d,i-1)
       end
      end
      a
    end

   def  makearray(d)
     sub(d,(d.length()-1))
   end

 



練習4.13 (配列の並べ替え)
a) def  sub_ascending(a)   ←補助関数。昇順になってないとこを数える
       j=0
       for  i  in 0..(a.length()-2)
         if  a[i]>a[i+1]
           j=j+1
         end
       end
       j
     end


     def  is_ascending(a)
       sub_is_ascending(a)==0
     end

b)  def swap_ascending(a,i)
   if a[i]>a[i+1]
   swap(a,i,i+1)
   end
  end

c) def  swap_ascending_all(a)
      for  i  in 0..(a.length()-2)
        swap_ascending(a,i)
      end
      a
    end

d) def  swapsort(a)
      while  !is_ascending(a)
         a=swap_ascending_all(a)
      end
      a
    end




練習4.14 (Eratosthenesのなんとか)
load("./make1d.rb")
def  eratosthenes(n)
  a=make1d(n-1)
  for  i  in 0..(n-2)
    if  a[i]==0
      for  j  in (i+1)..(n-2)
        if  (j+2)%(i+2)==0
          a[j]=1
        end
      end
    end
  end
  a
end


練習4.15 (Sierpinskiの三角形)
def  sierpinski(n)
  a=make2d(n,n)
  a[0][0]=1
  for  i  in  1..(n-1)
    a[i][0]=1
    for  j  in  1..(n-1)
      a[i][j]=(a[i-1][j-1]+a[i-1][j])%2
    end
  end
  a
end