[ruby-list:45859] Re: patch for Complex#sqrt in lib/cmath.rb of ruby-1.9.1-p0
From:
NISHIMATSU Takeshi <t_nissie@...>
Date:
2009-02-08 03:25:11 UTC
List:
ruby-list #45859
反応がなくてさみしい西松です。
--- I wrote in [ruby-list:45852]:
> 西松と申します。
> [ruby-list:44275] で触れた -0.0（符号付きゼロ）のためのパッチです。
> このときにちゃんとチェックせず、申し訳ありませんでした。
>
> % diff -u cmath.rb.original cmath.rb
> --- cmath.rb.original 2009-02-04 19:21:19.946793851 +0900
> +++ cmath.rb 2009-02-04 19:28:45.593345258 +0900
> @@ -63,7 +63,9 @@
> sqrt!(z)
> end
> else
> - if z.imag < 0
> + if (1.0/z.imag).infinite? == -1 and z.real < 0
> + Complex(0.0,-sqrt(-z.real))
> + elsif z.imag < 0
> sqrt(z.conjugate).conjugate
> else
> r = z.abs
このパッチでは、システムが「division by zero」をトラップ
するのかしないのかを考えていないので、ダメかもしれません。
参考: http://en.wikipedia.org/wiki/-0_(number)
C99のcopysign()と同様のモノがRubyにあれば
if x==0
if 1.copysgin(x)>0
# x is +0.0 or 0
else
# x=-0.0
end if
end if
などと書けるのに… missing/copysign.cが必要なんですよね。
ちなみにFortranにはsign(X,Y)があります。
もしかしたら誰かが必要とするかもしれないので、「division
by zero」版のsqrtとcprojとそれらのテストを最後に添付して
おきます。お騒がせしました。
--- cmath.ORIGINAL.rb 2008-12-13 10:03:25.000000000 +0900
+++ cmath.rb 2009-02-08 01:27:24.000000000 +0900
@@ -63,7 +63,9 @@
sqrt!(z)
end
else
- if z.imag < 0
+ if (1.0/z.imag).infinite? == -1 and z.real < 0
+ Complex(0.0,-sqrt(-z.real))
+ elsif z.imag < 0
sqrt(z.conjugate).conjugate
else
r = z.abs
@@ -73,6 +75,18 @@
end
end
+ def cproj(z)
+ if z.real.infinite? or z.image.infinite?
+ if z.imag>0.0 or (1.0/z.imag).infinite? == 1
+ Complex(1.0/0.0, 0.0)
+ else
+ Complex(1.0/0.0,-0.0)
+ end
+ else
+ z
+ end
+ end
+
def sin(z)
if z.real?
sin!(z)
@@ -190,6 +204,8 @@
module_function :sqrt!
module_function :sqrt
+ module_function :cproj
+
module_function :sin!
module_function :sin
module_function :cos!
--- test_complex.ORIGINAL.rb 2009-02-07 23:02:12.000000000 +0900
+++ test_complex.rb 2009-02-08 10:20:29.000000000 +0900
@@ -962,6 +962,43 @@
assert_in_delta(0.0, c.real, 0.001)
assert_in_delta(3.0, c.imag, 0.001)
+ c = Math.sqrt(Complex(-9,0.0))
+ assert_equal(0.0, c.real)
+ assert_equal( 3.0, c.imag)
+ c = Math.sqrt(Complex(-9,-0.0))
+ assert_equal(0.0, c.real)
+ assert_equal(-3.0, c.imag)
+
+ # normal
+ c = Complex(1.0, 2.0)
+ assert_equal(c,Math.cproj(c))
+ # upper half-plane
+ c = Complex( 1.0/0.0, 2)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal( 1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex(-1.0/0.0, 2)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal( 1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex( 1.0/0.0, 0.0)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal( 1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex(-0.0, 1.0/0.0)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal( 1,(1.0/Math.cproj(c).imag).infinite?)
+ # lower half-plane
+ c = Complex( 1.0/0.0, -2)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal(-1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex(-1.0/0.0, -2)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal(-1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex( 1.0/0.0, -0.0)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal(-1,(1.0/Math.cproj(c).imag).infinite?)
+ c = Complex( 2.0,-1.0/0.0)
+ assert_equal( 1,Math.cproj(c).real.infinite?)
+ assert_equal(-1,(1.0/Math.cproj(c).imag).infinite?)
+
c = Math.exp(Complex(1, 2))
assert_in_delta(-1.131, c.real, 0.001)
assert_in_delta(2.471, c.imag, 0.001)
--
西松タケシ
love && peace && free_software
http://loto.sourceforge.net/feram/ Fast MD program for perovskite-type ferroelectrics