ruby-dev

> もっとも、こうした話は微妙であることは確かなので、
> ComplexFloat をことさら強く推すというわけではなく、
> Complex の中で Float のパフォーマンスがよくなってくれさえすれば
> Complex でもいいかな、という気もしています。

最初のメールでも言ったことですが、そういったことは可能だと思います。

Complex での浮動小数点数のあつかいですが、村田が勧める C99 試案のコー
ドを試しにつかうパッチを添付します。

少しテストした結果からは、大方の予測通りだと思いますが、ComplexFloat 
が速くなっているのは、メソッド呼び出しの軽減によるところが大きいと思わ
れます。Complex も同様の傾向で速度が向上しています。データの持ち方を変
えるなどすれば、もっと速くなるのかなと思います。

ついでに python も少し試してみましたが、ruby の Complex よりも速そうで
すが、それほど極端に速いわけではないようです。perl の Math::Complex は、
ruby 1.8 と比べても、かなり遅いようです。いづれも、内容は見ておらず、
きちんとした性能の評価ではなく、単純な計算を繰り返したりしてみただけで
す。

田中昌宏さんたちが評価してくれるならば、これを取り込んでもよいと思いま
す。ただ、たとえ、C コンパイラが複素数を最適化しようともデータを ruby 
や python で転がしている限り、殆ど速度は変らないということです。おそら
く、計算方法としては、よく練られているのだと思いますから、そういう面で
は意味があるのだと思います。

Attachments (1)

diff (7.47 KB, text/x-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:

Thread

Prev Next

In This Thread

Prev Next