[前][次][番号順一覧][スレッド一覧]

ruby-changes:52976

From: nobu <ko1@a...>
Date: Sat, 20 Oct 2018 11:49:26 +0900 (JST)
Subject: [ruby-changes:52976] nobu:r65190 (trunk): complex.c: small optimization of Complex#**

nobu	2018-10-20 11:49:18 +0900 (Sat, 20 Oct 2018)

  New Revision: 65190

  https://svn.ruby-lang.org/cgi-bin/viewvc.cgi?view=revision&revision=65190

  Log:
    complex.c: small optimization of Complex#**
    
    * complex.c (rb_complex_pow): calculate power of a Fixnum without
      allocating intermediate Complex objects, and avoid unexpected
      NaNs.

  Modified files:
    trunk/complex.c
    trunk/internal.h
    trunk/numeric.c
    trunk/rational.c
    trunk/test/ruby/test_complex.rb
Index: internal.h
===================================================================
--- internal.h	(revision 65189)
+++ internal.h	(revision 65190)
@@ -1174,6 +1174,7 @@ VALUE rb_complex_mul(VALUE, VALUE); https://github.com/ruby/ruby/blob/trunk/internal.h#L1174
 VALUE rb_complex_abs(VALUE x);
 VALUE rb_complex_sqrt(VALUE x);
 VALUE rb_dbl_complex_polar(double abs, double ang);
+VALUE rb_complex_pow(VALUE self, VALUE other);
 
 /* cont.c */
 VALUE rb_obj_is_fiber(VALUE);
@@ -1471,6 +1472,7 @@ VALUE rb_int_abs(VALUE num); https://github.com/ruby/ruby/blob/trunk/internal.h#L1472
 VALUE rb_int_odd_p(VALUE num);
 int rb_int_positive_p(VALUE num);
 int rb_int_negative_p(VALUE num);
+VALUE rb_num_pow(VALUE x, VALUE y);
 
 static inline VALUE
 rb_num_compare_with_zero(VALUE num, ID mid)
@@ -1735,6 +1737,7 @@ VALUE rb_rational_reciprocal(VALUE x); https://github.com/ruby/ruby/blob/trunk/internal.h#L1737
 VALUE rb_cstr_to_rat(const char *, int);
 VALUE rb_rational_abs(VALUE self);
 VALUE rb_rational_cmp(VALUE self, VALUE other);
+VALUE rb_rational_pow(VALUE self, VALUE other);
 VALUE rb_numeric_quo(VALUE x, VALUE y);
 
 /* re.c */
Index: rational.c
===================================================================
--- rational.c	(revision 65189)
+++ rational.c	(revision 65190)
@@ -1001,8 +1001,8 @@ f_odd_p(VALUE integer) https://github.com/ruby/ruby/blob/trunk/rational.c#L1001
  *    Rational(1, 2) ** 0               #=> (1/1)
  *    Rational(1, 2) ** 0.0             #=> 1.0
  */
-static VALUE
-nurat_expt(VALUE self, VALUE other)
+VALUE
+rb_rational_pow(VALUE self, VALUE other)
 {
     if (k_numeric_p(other) && k_exact_zero_p(other))
 	return f_rational_new_bang1(CLASS_OF(self), ONE);
@@ -1077,6 +1077,7 @@ nurat_expt(VALUE self, VALUE other) https://github.com/ruby/ruby/blob/trunk/rational.c#L1077
 	return rb_num_coerce_bin(self, other, rb_intern("**"));
     }
 }
+#define nurat_expt rb_rational_pow
 
 /*
  * call-seq:
Index: numeric.c
===================================================================
--- numeric.c	(revision 65189)
+++ numeric.c	(revision 65190)
@@ -4077,6 +4077,22 @@ rb_int_pow(VALUE x, VALUE y) https://github.com/ruby/ruby/blob/trunk/numeric.c#L4077
     return Qnil;
 }
 
+VALUE
+rb_num_pow(VALUE x, VALUE y)
+{
+    VALUE z = rb_int_pow(x, y);
+    if (!NIL_P(z)) return z;
+    if (RB_FLOAT_TYPE_P(x)) return rb_float_pow(x, y);
+    if (SPECIAL_CONST_P(x)) return Qnil;
+    switch (BUILTIN_TYPE(x)) {
+      case T_COMPLEX:
+        return rb_complex_pow(x, y);
+      case T_RATIONAL:
+        return rb_rational_pow(x, y);
+    }
+    return Qnil;
+}
+
 /*
  * Document-method: Integer#==
  * Document-method: Integer#===
Index: test/ruby/test_complex.rb
===================================================================
--- test/ruby/test_complex.rb	(revision 65189)
+++ test/ruby/test_complex.rb	(revision 65190)
@@ -407,6 +407,10 @@ class Complex_Test < Test::Unit::TestCas https://github.com/ruby/ruby/blob/trunk/test/ruby/test_complex.rb#L407
     r = c ** Rational(-2,3)
     assert_in_delta(0.432, r.real, 0.001)
     assert_in_delta(-0.393, r.imag, 0.001)
+
+    c = Complex(0.0, -888888888888888.0)**8888
+    assert_not_predicate(c.real, :nan?)
+    assert_not_predicate(c.imag, :nan?)
   end
 
   def test_cmp
Index: complex.c
===================================================================
--- complex.c	(revision 65189)
+++ complex.c	(revision 65190)
@@ -725,6 +725,19 @@ safe_mul(VALUE a, VALUE b, int az, int b https://github.com/ruby/ruby/blob/trunk/complex.c#L725
     return f_mul(a, b);
 }
 
+static void
+comp_mul(VALUE areal, VALUE aimag, VALUE breal, VALUE bimag, VALUE *real, VALUE *imag)
+{
+    int arzero = f_zero_p(areal);
+    int aizero = f_zero_p(aimag);
+    int brzero = f_zero_p(breal);
+    int bizero = f_zero_p(bimag);
+    *real = f_sub(safe_mul(areal, breal, arzero, brzero),
+                  safe_mul(aimag, bimag, aizero, bizero));
+    *imag = f_add(safe_mul(areal, bimag, arzero, bizero),
+                  safe_mul(aimag, breal, aizero, brzero));
+}
+
 /*
  * call-seq:
  *    cmp * numeric  ->  complex
@@ -742,19 +755,9 @@ rb_complex_mul(VALUE self, VALUE other) https://github.com/ruby/ruby/blob/trunk/complex.c#L755
 {
     if (RB_TYPE_P(other, T_COMPLEX)) {
 	VALUE real, imag;
-	VALUE areal, aimag, breal, bimag;
-	int arzero, aizero, brzero, bizero;
-
 	get_dat2(self, other);
 
-	arzero = f_zero_p(areal = adat->real);
-	aizero = f_zero_p(aimag = adat->imag);
-	brzero = f_zero_p(breal = bdat->real);
-	bizero = f_zero_p(bimag = bdat->imag);
-	real = f_sub(safe_mul(areal, breal, arzero, brzero),
-		     safe_mul(aimag, bimag, aizero, bizero));
-	imag = f_add(safe_mul(areal, bimag, arzero, bizero),
-		     safe_mul(aimag, breal, aizero, brzero));
+        comp_mul(adat->real, adat->imag, bdat->real, bdat->imag, &real, &imag);
 
 	return f_complex_new2(CLASS_OF(self), real, imag);
     }
@@ -867,8 +870,8 @@ f_reciprocal(VALUE x) https://github.com/ruby/ruby/blob/trunk/complex.c#L870
  *    Complex('i') ** 2              #=> (-1+0i)
  *    Complex(-8) ** Rational(1, 3)  #=> (1.0000000000000002+1.7320508075688772i)
  */
-static VALUE
-nucomp_expt(VALUE self, VALUE other)
+VALUE
+rb_complex_pow(VALUE self, VALUE other)
 {
     if (k_numeric_p(other) && k_exact_zero_p(other))
 	return f_complex_new_bang1(CLASS_OF(self), ONE);
@@ -898,38 +901,45 @@ nucomp_expt(VALUE self, VALUE other) https://github.com/ruby/ruby/blob/trunk/complex.c#L901
 	return f_complex_polar(CLASS_OF(self), nr, ntheta);
     }
     if (FIXNUM_P(other)) {
-	if (f_gt_p(other, ZERO)) {
-	    VALUE x, z;
-	    long n;
-
-	    x = self;
-	    z = x;
-	    n = FIX2LONG(other) - 1;
-
-	    while (n) {
-		long q, r;
-
-		while (1) {
-		    get_dat1(x);
-
-		    q = n / 2;
-		    r = n % 2;
-
-		    if (r)
-			break;
-
-		    x = nucomp_s_new_internal(CLASS_OF(self),
-				       f_sub(f_mul(dat->real, dat->real),
-					     f_mul(dat->imag, dat->imag)),
-				       f_mul(f_mul(TWO, dat->real), dat->imag));
-		    n = q;
-		}
-		z = f_mul(z, x);
-		n--;
-	    }
-	    return z;
+        long n = FIX2LONG(other);
+        if (n == 0) {
+	    return nucomp_s_new_internal(CLASS_OF(self), ONE, ZERO);
+        }
+        if (n < 0) {
+            self = f_reciprocal(self);
+            other = rb_int_uminus(other);
+            n = -n;
+        }
+        {
+            get_dat1(self);
+            VALUE xr = dat->real, xi = dat->imag, zr = xr, zi = xi;
+
+            if (f_zero_p(xi)) {
+                zr = rb_num_pow(zr, other);
+            }
+            else if (f_zero_p(xr)) {
+                zi = rb_num_pow(zi, other);
+                if (n & 2) zi = f_negate(zi);
+                if (!(n & 1)) {
+                    VALUE tmp = zr;
+                    zr = zi;
+                    zi = tmp;
+                }
+            }
+            else {
+                while (--n) {
+                    long q, r;
+
+                    for (; q = n / 2, r = n % 2, r == 0; n = q) {
+                        VALUE tmp = f_sub(f_mul(xr, xr), f_mul(xi, xi));
+                        xi = f_mul(f_mul(TWO, xr), xi);
+                        xr = tmp;
+                    }
+                    comp_mul(zr, zi, xr, xi, &zr, &zi);
+                }
+            }
+	    return nucomp_s_new_internal(CLASS_OF(self), zr, zi);
 	}
-	return f_expt(f_reciprocal(self), rb_int_uminus(other));
     }
     if (k_numeric_p(other) && f_real_p(other)) {
 	VALUE r, theta;
@@ -945,6 +955,7 @@ nucomp_expt(VALUE self, VALUE other) https://github.com/ruby/ruby/blob/trunk/complex.c#L955
     }
     return rb_num_coerce_bin(self, other, id_expt);
 }
+#define nucomp_expt rb_complex_pow
 
 /*
  * call-seq:

--
ML: ruby-changes@q...
Info: http://www.atdot.net/~ko1/quickml/

[前][次][番号順一覧][スレッド一覧]