From: Tadashi Saito <shiba@...2.accsnet.ne.jp>
Date: 2003-07-08T17:47:10+09:00
Subject: [ruby-dev:20552] Re: [BigDecimal] renaming proposal

This is a multi-part message in MIME format.

------=_NextPart_000_00E9_01C34578.F9AA5E70
Content-Type: text/plain;
	charset="iso-2022-jp"
Content-Transfer-Encoding: 7bit

斎藤です。なにか誤解していたみたいです。

----- Original Message ----- 
From: "Shigeo Kobayashi" <shigeo@tinyforest.gr.jp>
Sent: Tuesday, July 08, 2003 4:59 PM
Subject: [ruby-dev:20550] Re: [BigDecimal] renaming proposal


> > 「現状のbigdecimal.cの数学関数の部分は削除して、Rubyで
> > 　書き直し、それを標準添付の一部としたい」
> >
> > ということだと思ったのですが、そうではない、ということでしょうか。
> いや、そういうことを意味したかったんですが．．．
>
> > 数学関数を標準添付したくないということですか?
> 標準添付できたらいいなと思っています。

了解です。自分が誤解していたらしく、すみませんでした。

とりあえず今まで作った分を、添付してみました。
[ruby-dev:20547]で添付したテストも、少なくともN=700くらいまでは
OKになることを確認済みです。

ただ自分は小林さんと違い、あまり数学に造詣が深くありません。
むしろ苦手だった方です。
今回のスクリプトも、奥村先生のアルゴリズム辞典や
数学の本を片手にほぼそのまま、コピーしただけです。
ですので、いろいろと修正して下さればありがたいです。
(自分も安直にテーラー展開(正確にはマクローリン展開?)した
だけなので…)

> 要するにBigDecimal::E(n)等が、正確に n 桁を保証できてい
> ないということだと思います。（計算内部の）中間変数の精度
> をもう少し上げればましになるとは思いますが、x**n/(n!) のよう
> な計算を繰り返し実行する単純なテーラー展開では所詮こん
> なものかもしれませんね。　やはり、削除したくなってきた！

いやいや、落ち着いてください ^^;;

自分も、正確な n 桁を保証するという方法のが分かりませんでした。
そこでとりあえず、ループ中の足し算の中で
「もし以前の結果と足した後の結果の n 桁が等しければ
　ループを抜ける」
という方法を取りました。が、数学的にはとてもあやしいと思います。

> > 手元で作ったスクリプトはもちろんそうしてありますが、
> > BigDecimal自体について言っているつもりでした。
> もちろん、 bigdecimal.c 内でより良いものがサポートできれば
> 申し分ないです。ただ、現状の関数は我ながらお粗末なので、
> できれば bigdecimal.c から E や sincos 等を削除したいと思っ
> ているわけです。

いやいや、そこまでおっしゃるほど悪いものではないと思いますよ。

添付したものは結構冗長なものになってしまいました。
EとPIはまだ実装していません。
とりあえずのたたき台バージョンということで。
おまけに、BigDecimal#to_fを再定義しています。
この話題は、また別のスレッドで。

---
斎藤ただし

------=_NextPart_000_00E9_01C34578.F9AA5E70
Content-Type: application/octet-stream;
	name="bigdecimal-math.rb"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
	filename="bigdecimal-math.rb"

require 'bigdecimal'

class BigDecimal
  def to_f
    f, x, y, z = self.split
    s = '0.' + x
    f * s.to_f * y ** z
  end
end

module Math
  alias sqrt! sqrt
  alias log! log
  alias sin! sin
  alias cos! cos
  alias exp! exp

  def sqrt(x, n = nil)
    maxnr = 100

    unless x.class == BigDecimal
      raise ArgumentError, 'wrong number of arguments(2 for 1)' unless n.nil?
      return sqrt!(x)
    end

    raise ArgumentError, "argument out of domain" if x < 0
    return x if x.zero? || x.infinite? || x.nan?
    
    if n.nil?
      n = [x.prec[0], x.exponent].max
    elsif n <= 0
      raise ArgumentError, "zero or negative argument"
    end
    
    two = BigDecimal.new('2')
    unless x > Float::MAX
      f = sqrt!(x.to_f)
      s = BigDecimal.new(f.to_s)
    else
      one = BigDecimal.new('1')
      s = (x > one) ? x : one
    end
    exp = s.exponent

    t = nil
    maxnr.times do
      last = s
      s = (x/s + s) / two
      t = s.truncate(n-exp)
      break if s == last || (s.prec[0] >= n && t == last.truncate(n-exp))
    end
    t
  end

  def log(x, n = nil)
    maxnr = 100

    unless x.class == BigDecimal
      raise ArgumentError, 'wrong number of arguments(2 for 1)' unless n.nil?
      return log!(x)
    end

    raise ArgumentError, "argument out of domain" if x <= 0
    return x if x.infinite? || x.nan?
    
    if n.nil?
      n = [x.prec[0], x.exponent].max
    elsif n <= 0
      raise ArgumentError, "zero or negative argument"
    end
    
    u = (x - 1) / (x + 1)
    u2 = u * u
    s = 2 * u
    exp = s.exponent

    t = nil
    i = 1
    maxnr.times do
      last = s
      u *= u2
      i += 2
      s += 2 * u / i
      t = s.truncate(n-exp)
      break if s == last || (s.prec[0] >= n && t == last.truncate(n-exp))
    end
    t
  end

  def sin(x, n = nil)
    maxnr = 100

    unless x.class == BigDecimal
      raise ArgumentError, 'wrong number of arguments(2 for 1)' unless n.nil?
      return sin!(x)
    end
    
    zero = BigDecimal.new('0')
    nan = zero / zero
    return nan if x.infinite? || x.nan?
    
    if n.nil?
      n = [x.prec[0], x.exponent].max
    elsif n <= 0
      raise ArgumentError, "zero or negative argument"
    end
    
    one = BigDecimal.new('1')
    two = BigDecimal.new('2')
    x1 = x
    x2 = x * x
    sign = one
    y = x
    m = one
    d = one

    z = nil
    maxnr.times do
      last = y
      sign = -sign
      x1 *= x2
      m += two
      d *= (m-one) * m
      y += sign * x1 / d
      z = y.truncate(n)
      break if y == last || (y.prec[0] >= n && z == last.truncate(n))
    end
    z
  end

  def cos(x, n = nil)
    maxnr = 100

    unless x.class == BigDecimal
      raise ArgumentError, 'wrong number of arguments(2 for 1)' unless n.nil?
      return cos!(x)
    end
    
    zero = BigDecimal.new('0')
    nan = zero / zero
    return nan if x.infinite? || x.nan?
    
    if n.nil?
      n = [x.prec[0], x.exponent].max
    elsif n <= 0
      raise ArgumentError, "zero or negative argument"
    end
    
    one = BigDecimal.new('1')
    two = BigDecimal.new('2')
    x1 = one
    x2 = x * x
    sign = one
    y = one
    m = zero
    d = one

    z = nil
    maxnr.times do
      last = y
      sign = -sign
      x1 *= x2
      m += two
      d *= (m-one) * m
      y += sign * x1 / d
      z = y.truncate(n)
      break if y == last || (y.prec[0] >= n && z == last.truncate(n))
    end
    z
  end

  def exp(x, n = nil)
    maxnr = 100

    unless x.class == BigDecimal
      raise ArgumentError, 'wrong number of arguments(2 for 1)' unless n.nil?
      return exp!(x)
    end
    
    zero = BigDecimal.new('0')
    return zero if x.infinite? == -1
    return x if x.nan? || x.infinite?
    
    if n.nil?
      n = [x.prec[0], x.exponent].max
    elsif n <= 0
      raise ArgumentError, "zero or negative argument"
    end
    
    one = BigDecimal.new('1')
    x1 = one
    y = one
    d = one
    m = zero

    z = nil
    maxnr.times do
      last = y
      x1 *= x
      m += one
      d *= m
      y += x1 / d
      z = y.truncate(n)
      break if y == last || (y.prec[0] >= n && z == last.truncate(n))
    end
    z
  end

  module_function :sqrt, :sqrt!
  module_function :log, :log!
  module_function :sin, :sin!
  module_function :cos, :cos!
  module_function :exp, :exp!

end

#sqrt(BigDecimal.new(ARGV[0]))

------=_NextPart_000_00E9_01C34578.F9AA5E70--