From: Tadayoshi Funaba Date: 2008-05-14T23:24:57+09:00 Subject: [ruby-dev:34706] Re: int/int -> rational ----Next_Part(Wed_May_14_23_23_46_2008_944)-- Content-Type: Text/Plain; charset=iso-2022-jp Content-Transfer-Encoding: 7bit > もっとも、こうした話は微妙であることは確かなので、 > ComplexFloat をことさら強く推すというわけではなく、 > Complex の中で Float のパフォーマンスがよくなってくれさえすれば > Complex でもいいかな、という気もしています。 最初のメールでも言ったことですが、そういったことは可能だと思います。 Complex での浮動小数点数のあつかいですが、村田が勧める C99 試案のコー ドを試しにつかうパッチを添付します。 少しテストした結果からは、大方の予測通りだと思いますが、ComplexFloat が速くなっているのは、メソッド呼び出しの軽減によるところが大きいと思わ れます。Complex も同様の傾向で速度が向上しています。データの持ち方を変 えるなどすれば、もっと速くなるのかなと思います。 ついでに python も少し試してみましたが、ruby の Complex よりも速そうで すが、それほど極端に速いわけではないようです。perl の Math::Complex は、 ruby 1.8 と比べても、かなり遅いようです。いづれも、内容は見ておらず、 きちんとした性能の評価ではなく、単純な計算を繰り返したりしてみただけで す。 田中昌宏さんたちが評価してくれるならば、これを取り込んでもよいと思いま す。ただ、たとえ、C コンパイラが複素数を最適化しようともデータを ruby や python で転がしている限り、殆ど速度は変らないということです。おそら く、計算方法としては、よく練られているのだと思いますから、そういう面で は意味があるのだと思います。 ----Next_Part(Wed_May_14_23_23_46_2008_944)-- Content-Type: Text/Plain; charset=us-ascii Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename=diff Index: complex.c =================================================================== --- complex.c (revision 16344) +++ complex.c (working copy) @@ -202,7 +202,32 @@ fun1(numerator) fun1(polar) fun1(scalar_p) + +#if 0 fun1(to_f) +#else +inline static VALUE +f_to_f(VALUE x) +{ + VALUE r; + switch (TYPE(x)) { + case T_FIXNUM: + r = rb_float_new((double)FIX2LONG(x)); + break; + case T_BIGNUM: + r = rb_float_new(rb_big2dbl(x)); + break; + case T_FLOAT: + r = x; + break; + default: + r = rb_funcall(x, id_to_f, 0); + break; + } + return r; +} +#endif + fun1(to_i) fun1(to_r) fun1(to_s) @@ -585,7 +610,59 @@ return dat->image; } +#if defined(__STDC_IEC_559__) && defined(INFINITY) && defined(NAN) +#define C99COMP +#endif + +inline static double +c99_creal(VALUE x) +{ + switch (TYPE(x)) { + case T_FIXNUM: + return (double)FIX2LONG(x); + case T_BIGNUM: + return rb_big2dbl(x); + case T_FLOAT: + return RFLOAT_VALUE(x); + case T_RATIONAL: + return RFLOAT_VALUE(f_to_f(x)); + case T_COMPLEX: + { + get_dat1(x); + return RFLOAT_VALUE(f_to_f(dat->real)); + } + } +} + +inline static double +c99_cimag(VALUE x) +{ + switch (TYPE(x)) { + case T_FIXNUM: + case T_BIGNUM: + case T_FLOAT: + case T_RATIONAL: + return 0.0; + case T_COMPLEX: + { + get_dat1(x); + return RFLOAT_VALUE(f_to_f(dat->image)); + } + } +} + static VALUE +c99_add(VALUE z, VALUE w) +{ + double a, b, c, d; + a = c99_creal(z); b = c99_cimag(z); + c = c99_creal(w); d = c99_cimag(w); + return f_complex_new_bang2(CLASS_OF(z), + rb_float_new(a + c), + rb_float_new(b + d)); +} + +static VALUE nucomp_add(VALUE self, VALUE other) { switch (TYPE(other)) { @@ -605,6 +682,13 @@ get_dat2(self, other); +#ifdef C99COMP + if ((TYPE(adat->real) == T_FLOAT || + TYPE(bdat->real) == T_FLOAT) && + (TYPE(adat->image) == T_FLOAT || + TYPE(bdat->image) == T_FLOAT)) + return c99_add(self, other); +#endif real = f_add(adat->real, bdat->real); image = f_add(adat->image, bdat->image); @@ -616,6 +700,17 @@ } static VALUE +c99_sub(VALUE z, VALUE w) +{ + double a, b, c, d; + a = c99_creal(z); b = c99_cimag(z); + c = c99_creal(w); d = c99_cimag(w); + return f_complex_new_bang2(CLASS_OF(z), + rb_float_new(a - c), + rb_float_new(b - d)); +} + +static VALUE nucomp_sub(VALUE self, VALUE other) { switch (TYPE(other)) { @@ -635,6 +730,13 @@ get_dat2(self, other); +#ifdef C99COMP + if ((TYPE(adat->real) == T_FLOAT || + TYPE(bdat->real) == T_FLOAT) && + (TYPE(adat->image) == T_FLOAT || + TYPE(bdat->image) == T_FLOAT)) + return c99_sub(self, other); +#endif real = f_sub(adat->real, bdat->real); image = f_sub(adat->image, bdat->image); @@ -646,12 +748,66 @@ } static VALUE +c99_mul(VALUE z, VALUE w) +{ + double a, b, c, d, ac, bd, ad, bc, x, y; + a = c99_creal(z); b = c99_cimag(z); + c = c99_creal(w); d = c99_cimag(w); + ac = a * c; bd = b * d; + ad = a * d; bc = b * c; + x = ac - bd; + y = ad + bc; + /* Recover infinities that computed as NaN+iNaN ... */ + if (isnan(x) && isnan(y)) { + int recalc = 0; + if ( isinf(a) || isinf(b) ) { /* z is infinite */ + /* "Box" the infinity ... */ + a = copysign(isinf(a) ? 1.0 : 0.0, a); + b = copysign(isinf(b) ? 1.0 : 0.0, b); + /* Change NaNs in the other factor to 0 ... */ + if (isnan(c)) c = copysign(0.0, c); + if (isnan(d)) d = copysign(0.0, d); + recalc = 1; + } + if ( isinf(c) || isinf(d) ) { /* w is infinite */ + /* "Box" the infinity ... */ + c = copysign(isinf(c) ? 1.0 : 0.0, c); + d = copysign(isinf(d) ? 1.0 : 0.0, d); + /* Change NaNs in the other factor to 0 ... */ + if (isnan(a)) a = copysign(0.0, a); + if (isnan(b)) b = copysign(0.0, b); + recalc = 1; + } + if (!recalc) { + /* *Recover infinities from overflow cases ... */ + if (isinf(ac) || isinf(bd) || + isinf(ad) || isinf(bc)) { + /* Change all NaNs to 0 ... */ + if (isnan(a)) a = copysign(0.0, a); + if (isnan(b)) b = copysign(0.0, b); + if (isnan(c)) c = copysign(0.0, c); + if (isnan(d)) d = copysign(0.0, d); + recalc = 1; + } + } + if (recalc) { + x = INFINITY * ( a * c - b * d ); + y = INFINITY * ( a * d + b * c ); + } + } + return f_complex_new_bang2(CLASS_OF(z), rb_float_new(x), rb_float_new(y)); +} + +static VALUE nucomp_mul(VALUE self, VALUE other) { switch (TYPE(other)) { + case T_FLOAT: +#ifdef C99COMP + return c99_mul(self, other); +#endif case T_FIXNUM: case T_BIGNUM: - case T_FLOAT: case T_RATIONAL: { get_dat1(self); @@ -666,6 +822,13 @@ get_dat2(self, other); +#ifdef C99COMP + if (TYPE(adat->real) == T_FLOAT || + TYPE(adat->image) == T_FLOAT || + TYPE(bdat->real) == T_FLOAT || + TYPE(bdat->image) == T_FLOAT) + return c99_mul(self, other); +#endif real = f_sub(f_mul(adat->real, bdat->real), f_mul(adat->image, bdat->image)); image = f_add(f_mul(adat->real, bdat->image), @@ -678,13 +841,62 @@ } } + static VALUE +c99_div(VALUE z, VALUE w) +{ + double a, b, c, d, logbw, denom, x, y; + int ilogbw = 0; + a = c99_creal(z); b = c99_cimag(z); + c = c99_creal(w); d = c99_cimag(w); + logbw = logb(fmax(fabs(c), fabs(d))); + if (isfinite(logbw)) { + ilogbw = (int)logbw; + c = scalbn(c, -ilogbw); + d = scalbn(d, -ilogbw); + } + denom = c * c + d * d; + x = scalbn((a * c + b * d) / denom, -ilogbw); + y = scalbn((b * c - a * d) / denom, -ilogbw); + /* + * Recover infinities and zeros that computed + * as NaN+iNaN; the only cases are non-zero/zero, + * infinite/finite, and finite/infinite, ... + */ + if (isnan(x) && isnan(y)) { + if ((denom == 0.0) && + (!isnan(a) || !isnan(b))) { + x = copysign(INFINITY, c) * a; + y = copysign(INFINITY, c) * b; + } + else if ((isinf(a) || isinf(b)) && + isfinite(c) && isfinite(d)) { + a = copysign(isinf(a) ? 1.0 : 0.0, a); + b = copysign(isinf(b) ? 1.0 : 0.0, b); + x = INFINITY * ( a * c + b * d ); + y = INFINITY * ( b * c - a * d ); + } + else if (isinf(logbw) && + isfinite(a) && isfinite(b)) { + c = copysign(isinf(c) ? 1.0 : 0.0, c); + d = copysign(isinf(d) ? 1.0 : 0.0, d); + x = 0.0 * ( a * c + b * d ); + y = 0.0 * ( b * c - a * d ); + } + } + return f_complex_new_bang2(CLASS_OF(z), rb_float_new(x), rb_float_new(y)); +} + +static VALUE nucomp_div(VALUE self, VALUE other) { switch (TYPE(other)) { + case T_FLOAT: +#ifdef C99COMP + return c99_div(self, other); +#endif case T_FIXNUM: case T_BIGNUM: - case T_FLOAT: case T_RATIONAL: { get_dat1(self); @@ -697,16 +909,13 @@ { get_dat2(self, other); - if (TYPE(adat->real) == T_FLOAT || - TYPE(adat->image) == T_FLOAT || - TYPE(bdat->real) == T_FLOAT || - TYPE(bdat->image) == T_FLOAT) { - VALUE magn = m_hypot(bdat->real, bdat->image); - VALUE tmp = f_complex_new_bang2(CLASS_OF(self), - f_div(bdat->real, magn), - f_div(bdat->image, magn)); - return f_div(f_mul(self, f_conjugate(tmp)), magn); - } +#ifdef C99COMP + if (TYPE(adat->real) == T_FLOAT || + TYPE(adat->image) == T_FLOAT || + TYPE(bdat->real) == T_FLOAT || + TYPE(bdat->image) == T_FLOAT) + return c99_div(self, other); +#endif return f_div(f_mul(self, f_conjugate(other)), f_abs2(other)); } default: ----Next_Part(Wed_May_14_23_23_46_2008_944)----