diff options
| -rw-r--r-- | ChangeLog | 5 | ||||
| -rw-r--r-- | complex.c | 34 | ||||
| -rw-r--r-- | internal.h | 4 | ||||
| -rw-r--r-- | lib/mathn.rb | 9 | ||||
| -rw-r--r-- | math.c | 15 |
5 files changed, 45 insertions, 22 deletions
@@ -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 @@ -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 1a75532e79..b47430cd7d 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 d6f2da0210..9cb9401e37 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 @@ -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"); |