summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authornobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2016-07-12 14:13:46 +0000
committernobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2016-07-12 14:13:46 +0000
commit745a2aac690ef026ffd8a3c1ca82ddc61d6babe6 (patch)
tree30b7be7a49a2c78b48488a0ceb58fd1530000ce5
parent89ed4f41a3cc169893ee50765e7e996a63bbd2e9 (diff)
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
-rw-r--r--ChangeLog5
-rw-r--r--complex.c34
-rw-r--r--internal.h4
-rw-r--r--lib/mathn.rb9
-rw-r--r--math.c15
5 files changed, 45 insertions, 22 deletions
diff --git a/ChangeLog b/ChangeLog
index 8bc642c760f..66c160b4d18 100644
--- a/ChangeLog
+++ b/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
diff --git a/complex.c b/complex.c
index 9419ea98186..fc6ecbb0e09 100644
--- a/complex.c
+++ b/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 {
- 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)));
- }
- }
+ return rb_complex_sqrt(x);
}
#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
diff --git a/internal.h b/internal.h
index 1a75532e790..b47430cd7da 100644
--- a/internal.h
+++ b/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);
diff --git a/lib/mathn.rb b/lib/mathn.rb
index d6f2da02100..9cb9401e37b 100644
--- a/lib/mathn.rb
+++ b/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)
- 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
+ sqrt!(a)
elsif a.respond_to?(:nan?) and a.nan?
a
elsif a >= 0
diff --git a/math.c b/math.c
index 750a520c9b4..8d36389a288 100644
--- a/math.c
+++ b/math.c
@@ -588,8 +588,23 @@ 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");