math.c: Complex sqrt

* math.c (rb_math_sqrt): [EXPERIMENTAL] move Complex sqrt support
  from mathn.rb.

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@55646 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

ChangeLog

@@ -1,3 +1,8 @@
+Tue Jul 12 23:13:43 2016  Nobuyoshi Nakada  <nobu@ruby-lang.org>
+
+	* math.c (rb_math_sqrt): [EXPERIMENTAL] move Complex sqrt support
+	  from mathn.rb.
+
 Tue Jul 12 21:59:40 2016  Martin Duerst  <duerst@it.aoyama.ac.jp>
 
 	* revert r55642 (previous commit) because of test failure at

complex.c

@@ -535,6 +535,21 @@ m_sin(VALUE x)
 #if 0
 imp1(sqrt)
 
+VALUE
+rb_complex_sqrt(VALUE x)
+{
+    int pos;
+    VALUE a, re, im;
+    get_dat1(x);
+
+    pos = f_positive_p(dat->imag);
+    a = f_abs(x);
+    re = m_sqrt_bang(f_div(f_add(a, dat->real), TWO));
+    im = m_sqrt_bang(f_div(f_sub(a, dat->real), TWO));
+    if (!pos) im = f_negate(im);
+    return f_complex_new2(rb_cComplex, re, im);
+}
+
 static VALUE
 m_sqrt(VALUE x)
 {
@@ -543,18 +558,7 @@ m_sqrt(VALUE x)
 	    return m_sqrt_bang(x);
 	return f_complex_new2(rb_cComplex, ZERO, m_sqrt_bang(f_negate(x)));
     }
-    else {
+    return rb_complex_sqrt(x);
-	get_dat1(x);
-
-	if (f_negative_p(dat->imag))
-	    return f_conj(m_sqrt(f_conj(x)));
-	else {
-	    VALUE a = f_abs(x);
-	    return f_complex_new2(rb_cComplex,
-				  m_sqrt_bang(f_div(f_add(a, dat->real), TWO)),
-				  m_sqrt_bang(f_div(f_sub(a, dat->real), TWO)));
-	}
-    }
 }
 #endif
 
@@ -1415,6 +1419,12 @@ rb_complex_set_imag(VALUE cmp, VALUE i)
     return cmp;
 }
 
+VALUE
+rb_complex_abs(VALUE cmp)
+{
+    return nucomp_abs(cmp);
+}
+
 /*
  * call-seq:
  *    cmp.to_i  ->  integer

internal.h

@@ -903,6 +903,8 @@ VALUE rb_insns_name_array(void);
 /* complex.c */
 VALUE rb_complex_plus(VALUE, VALUE);
 VALUE rb_complex_mul(VALUE, VALUE);
+VALUE rb_complex_abs(VALUE x);
+VALUE rb_complex_sqrt(VALUE x);
 
 /* cont.c */
 VALUE rb_obj_is_fiber(VALUE);
@@ -1083,9 +1085,7 @@ VALUE rb_math_hypot(VALUE, VALUE);
 VALUE rb_math_log(int argc, const VALUE *argv);
 VALUE rb_math_sin(VALUE);
 VALUE rb_math_sinh(VALUE);
-#if 0
 VALUE rb_math_sqrt(VALUE);
-#endif
 
 /* newline.c */
 void Init_newline(void);

lib/mathn.rb

@@ -103,14 +103,7 @@ module Math
 
   def sqrt(a)
     if a.kind_of?(Complex)
-      abs = sqrt(a.real*a.real + a.imag*a.imag)
+      sqrt!(a)
-      x = sqrt((a.real + abs)/Rational(2))
-      y = sqrt((-a.real + abs)/Rational(2))
-      if a.imag >= 0
-        Complex(x, y)
-      else
-        Complex(x, -y)
-      end
     elsif a.respond_to?(:nan?) and a.nan?
       a
     elsif a >= 0

math.c

@@ -587,9 +587,24 @@ math_log10(VALUE obj, VALUE x)
 
 static VALUE
 math_sqrt(VALUE obj, VALUE x)
+{
+    return rb_math_sqrt(x);
+}
+
+VALUE
+rb_math_sqrt(VALUE x)
 {
     double d;
 
+    if (RB_TYPE_P(x, T_COMPLEX)) {
+	int neg = signbit(RCOMPLEX(x)->imag);
+	double re = Get_Double(RCOMPLEX(x)->real), im;
+	d = Get_Double(rb_complex_abs(x));
+	im = sqrt((d - re) / 2.0);
+	re = sqrt((d + re) / 2.0);
+	if (neg) im = -im;
+	return rb_complex_new(DBL2NUM(re), DBL2NUM(im));
+    }
     d = Get_Double(x);
     /* check for domain error */
     if (d < 0.0) domain_error("sqrt");