浮動小数点演算

スライド「実数データ表現」を体験するために，符号部を1ビット，指数部を5 ビット，仮数部を10ビットで表現した16ビット浮動小数点表現を考える(最近の PCのGPU(Graphic Processing Unit)ではこの16ビット浮動小数点表現が使われる こともある)．この3つ 組を大きさ3の配列にして表現したのが以下のプログラム( float16.rb)である． 
include Math

# 2を底にした対数を取り切り捨てる
def log2(x)
  (log(x)/log(2)).floor
end

# 16ビットの浮動小数点形式への変換
# 大きさ3の配列で
# [符号(0: 正, 1: 負), 指数部(1-62), 仮数部(0-1023)]
# 指数部が63の時はNaN(Not a Number 表現できる範囲外の数)
# 指数部が0の時は0
def toFloat16(f)
  if f<0
    sign=1
    f=-f;
  else
    sign=0
  end
  if f==0
    [0,0,0]
  else
    exp=log2(f)
    if(exp+31<=0)
      [sign,0,0]
    elsif(exp+31>=63)
      [sign,63,0]
    else
      f=f/2**exp
      [sign,exp+31,(f* 2**10).floor-1024]
    end
  end
end
# NaNに対応する定数
NAN=0.0/0.0
# 16ビットの浮動小数点形式からの変換
def fromFloat16(x)
  if(x[1]==0)
    (1-2*x[0])*0.0
  elsif(x[1]==63)
    NAN
  else
    (1-2*x[0])* 2.0**(x[1]-31) * (x[2]+1024) * 2.0**-10
  end
end
# Float16同士の加算
def add(x,y)
  toFloat16(fromFloat16(x)+fromFloat16(y))
end
# Float16同士の減算
def sub(x,y)
  toFloat16(fromFloat16(x)-fromFloat16(y))
end
# Float16同士の積算
def mul(x,y)
  toFloat16(fromFloat16(x)*fromFloat16(y))
end
# Float16同士の除算
def div(x,y)
  toFloat16(fromFloat16(x)/fromFloat16(y))
end
これを，ファイルにダウンロードして保存し，irbの中から 
irb
>> load "./float16.rb"
load "./float16.rb"
=> true
>> toFloat16(0.1) 
toFloat16(0.1)
=> [0, 27, 614]
>> fromFloat16([0,28,614])
fromFloat16([0,28,614])
=> 0.199951171875
>> fromFloat16(add(toFloat16(0.1),toFloat16(0.2)))
fromFloat16(add(toFloat16(0.1),toFloat16(0.2)))
=> 0.2998046875
のようにして使うことができる．

演習

丸め誤差

16ビット浮動小数点表現では丸め誤差を容易に体験できる． たとえば，0.4は一度16ビット浮動小数点表現で表して，Rubyの浮動小数点表現に戻してみると， 
fromFloat16(toFloat16(0.4))
=> 0.39990234375
のようになる．そのため，誤差なく表現できていないことが分かる．一方，0.25は 
fromFloat16(toFloat16(0.25))
=> 0.25
のように元に戻るので，「誤差はあるがたまたま64ビット表現でも区別できな いほど」というレアなケースを考えなければ，誤差無く表現できていると判断 できる．

演習

以下の数の中で，16ビット浮動小数点表現で表した時に誤差無く表現できない数は何か?

桁落ち誤差

スライド「桁落ち誤差の例」を見て，以下の「自然対数のx=1における微係数を求める」プログラムを動かしてみて途中から 誤差が増えていくことを確認しなさい． 
def f(x)
  return Math.log(x);   # f(x)=ln(x)
end

x=1.0;
for i in 1..15
  h =0.1**i     #h=0.1, 0.01, ..., 1e-15
  df=(f(x+h)-f(x))/h
  err=(df-1.0).abs;
  print "h=",h,"\t f'(",x,")=",df," err=",err,"\n"
end
以下の「sin(x)のx=0における微係数を求める」プログラムを動かしてみて，さ きほどの例との違いを観察しなさい．この差の原因を考えてみなさい(正確に答えるには当初考えたよりも深い考察が必要な問題だと分かったので，答えられる範囲で解答するように)． 
def f(x)
  return Math.sin(x);   # f(x)=ln(x)
end
x=0.0;
for i in 1..15
  h =0.1**i     #h=0.1, 0.01, , 1e-15
  df=(f(x+h)-f(x))/h
  err=(df-1.0).abs;
  print "h=",h,"\t f'(",x,")=",df," err=",err,"\n"
end

Euler法とRunge-Kutta法

dy/dx=1/2y
という常微分方程式を考える．y(1)=1という初期条件のもとで，y(2)の 値を求めることを試みる(こたえは√2になるはずである)．オイラー法で求め るプログラムは 
def euler(n)
  y = 1.0
  h = 1.0/n
  for i in 0..n-1
     y = y + (1.0/(2*y))*h
  end
  y
end
のようになる．一方，Runge-Kutta法では 
def runge_kutta(n)
  y = 1.0
  h = 1.0/n
  for i in 0..n-1
     k1= (1.0/(2*y))*h
     k2= (1.0/(2*(y+k1)))*h
     y=y+(k1+k2)/2.0
  end
  y
end
のようになる．

演習

nの値を1から10倍ずつ変化させていって，結果が√2と小数点以下何桁一致するかを 調べてみましょう． 
for i in (0..6)
 print "      euler(10**",i,")=",euler(10**i),"\n"
 print "runge_kutta(10**",i,")=",runge_kutta(10**i),"\n"
end