From bb6384761e6c32acd169faa337d77d7e8d1efe5e Mon Sep 17 00:00:00 2001 From: mrkn Date: Mon, 27 Jun 2011 16:26:09 +0000 Subject: * ext/bigdecimal/bigdecimal.c (BigMath_s_log): move BigMath.log from bigdecimal/math.rb. * ext/bigdecimal/lib/bigdecimal/math.rb: ditto. * test/bigdecimal/test_bigdecimal.rb: move test for BigMath.log from test/bigdecimal/test_bigmath.rb. * test/bigdecimal/test_bigmath.rb: ditto. git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@32258 b2dd03c8-39d4-4d8f-98ff-823fe69b080e --- ext/bigdecimal/bigdecimal.c | 191 ++++++++++++++++++++++++++++++++++ ext/bigdecimal/lib/bigdecimal/math.rb | 35 ------- 2 files changed, 191 insertions(+), 35 deletions(-) (limited to 'ext/bigdecimal') diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index 9e22dc55e1..bdeb12a5a2 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -69,6 +69,33 @@ static ID id_to_r; #define DBLE_FIG (DBL_DIG+1) /* figure of double */ #endif +#ifndef RBIGNUM_ZERO_P +# define RBIGNUM_ZERO_P(x) (RBIGNUM_LEN(x) == 0 || \ + (RBIGNUM_DIGITS(x)[0] == 0 && \ + (RBIGNUM_LEN(x) == 1 || bigzero_p(x)))) +#endif + +static inline int +bigzero_p(VALUE x) +{ + long i; + BDIGIT *ds = RBIGNUM_DIGITS(x); + + for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) { + if (ds[i]) return 0; + } + return 1; +} + +#ifndef RRATIONAL_ZERO_P +# define RRATIONAL_ZERO_P(x) (FIXNUM_P(RRATIONAL(x)->num) && \ + FIX2LONG(RRATIONAL(x)->num) == 0) +#endif + +#ifndef RRATIONAL_NEGATIVE_P +# define RRATIONAL_NEGATIVE_P(x) RTEST(rb_funcall((x), '<', 1, INT2FIX(0))) +#endif + /* * ================== Ruby Interface part ========================== */ @@ -2132,6 +2159,169 @@ BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec) } } +/* call-seq: + * BigMath.log(x, prec) + * + * Computes the natural logarithm of x to the specified number of digits of + * precision. + * + * If x is zero or negative, raises Math::DomainError. + * + * If x is positive infinite, returns Infinity. + * + * If x is NaN, returns NaN. + */ +static VALUE +BigMath_s_log(VALUE klass, VALUE x, VALUE vprec) +{ + ssize_t prec, n, i; + SIGNED_VALUE expo; + Real* vx = NULL; + VALUE argv[2], vn, one, two, w, x2, y, d; + int zero = 0; + int negative = 0; + int infinite = 0; + int nan = 0; + double flo; + long fix; + + if (TYPE(vprec) != T_FIXNUM && TYPE(vprec) != T_BIGNUM) { + rb_raise(rb_eArgError, "precision must be an Integer"); + } + + prec = NUM2SSIZET(vprec); + if (prec <= 0) { + rb_raise(rb_eArgError, "Zero or negative precision for exp"); + } + + /* TODO: the following switch statement is almostly the same as one in the + * BigDecimalCmp function. */ + switch (TYPE(x)) { + case T_DATA: + if (!is_kind_of_BigDecimal(x)) break; + vx = DATA_PTR(x); + zero = VpIsZero(vx); + negative = VpGetSign(vx) < 0; + infinite = VpIsPosInf(vx) || VpIsNegInf(vx); + nan = VpIsNaN(vx); + break; + + case T_FIXNUM: + fix = FIX2LONG(x); + zero = fix == 0; + negative = fix < 0; + goto get_vp_value; + + case T_BIGNUM: + zero = RBIGNUM_ZERO_P(x); + negative = RBIGNUM_NEGATIVE_P(x); +get_vp_value: + if (zero || negative) break; + vx = GetVpValue(x, 0); + break; + + case T_FLOAT: + flo = RFLOAT_VALUE(x); + zero = flo == 0; + negative = flo < 0; + infinite = isinf(flo); + nan = isnan(flo); + if (!zero && !negative && !infinite && !nan) { + vx = GetVpValueWithPrec(x, DBL_DIG+1, 1); + } + break; + + case T_RATIONAL: + zero = RRATIONAL_ZERO_P(x); + negative = RRATIONAL_NEGATIVE_P(x); + if (zero || negative) break; + vx = GetVpValueWithPrec(x, prec, 1); + break; + + case T_COMPLEX: + rb_raise(rb_eMathDomainError, + "Complex argument for BigMath.log"); + + default: + break; + } + if (infinite && !negative) { + Real* vy; + vy = VpCreateRbObject(prec, "#0"); + RB_GC_GUARD(vy->obj); + VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE); + return ToValue(vy); + } + else if (nan) { + Real* vy; + vy = VpCreateRbObject(prec, "#0"); + RB_GC_GUARD(vy->obj); + VpSetNaN(vy); + return ToValue(vy); + } + else if (zero || negative) { + rb_raise(rb_eMathDomainError, + "Zero or negative argument for log"); + } + else if (vx == NULL) { + rb_raise(rb_eArgError, "%s can't be coerced into BigDecimal", + rb_special_const_p(x) ? RSTRING_PTR(rb_inspect(x)) : rb_obj_classname(x)); + } + x = ToValue(vx); + + RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1")); + RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2")); + + n = prec + rmpd_double_figures(); + RB_GC_GUARD(vn) = SSIZET2NUM(n); + expo = VpExponent10(vx); + if (expo < 0 || expo >= 3) { + char buf[16]; + snprintf(buf, 16, "1E%ld", -expo); + x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn); + } + else { + expo = 0; + } + w = BigDecimal_sub(x, one); + argv[0] = BigDecimal_add(x, one); + argv[1] = vn; + x = BigDecimal_div2(2, argv, w); + RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn); + RB_GC_GUARD(y) = x; + RB_GC_GUARD(d) = y; + i = 1; + while (!VpIsZero((Real*)DATA_PTR(d))) { + SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y)); + SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d)); + ssize_t m = n - vabs(ey - ed); + if (m <= 0) { + break; + } + else if ((size_t)m < rmpd_double_figures()) { + m = rmpd_double_figures(); + } + + x = BigDecimal_mult2(x2, x, vn); + i += 2; + argv[0] = SSIZET2NUM(i); + argv[1] = SSIZET2NUM(m); + d = BigDecimal_div2(2, argv, x); + y = BigDecimal_add(y, d); + } + + y = BigDecimal_mult(y, two); + if (expo != 0) { + VALUE log10, vexpo, dy; + log10 = BigMath_s_log(klass, INT2FIX(10), vprec); + vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1)); + dy = BigDecimal_mult(log10, vexpo); + y = BigDecimal_add(y, dy); + } + + return y; +} + /* Document-class: BigDecimal * BigDecimal provides arbitrary-precision floating point decimal arithmetic. * @@ -2430,6 +2620,7 @@ Init_bigdecimal(void) /* mathematical functions */ rb_mBigMath = rb_define_module("BigMath"); rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2); + rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2); id_BigDecimal_exception_mode = rb_intern_const("BigDecimal.exception_mode"); id_BigDecimal_rounding_mode = rb_intern_const("BigDecimal.rounding_mode"); diff --git a/ext/bigdecimal/lib/bigdecimal/math.rb b/ext/bigdecimal/lib/bigdecimal/math.rb index dc557deef2..03c59bfccb 100644 --- a/ext/bigdecimal/lib/bigdecimal/math.rb +++ b/ext/bigdecimal/lib/bigdecimal/math.rb @@ -145,41 +145,6 @@ module BigMath y end - # Computes the natural logarithm of x to the specified number of digits - # of precision. - # - # Returns x if x is infinite or NaN. - # - def log(x, prec) - raise ArgumentError, "Zero or negative argument for log" if x <= 0 || prec <= 0 - return x if x.infinite? || x.nan? - one = BigDecimal("1") - two = BigDecimal("2") - n = prec + BigDecimal.double_fig - if (expo = x.exponent) < 0 || expo >= 3 - x = x.mult(BigDecimal("1E#{-expo}"), n) - else - expo = nil - end - x = (x - one).div(x + one,n) - x2 = x.mult(x,n) - y = x - d = y - i = one - while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) - m = BigDecimal.double_fig if m < BigDecimal.double_fig - x = x2.mult(x,n) - i += two - d = x.div(i,m) - y += d - end - y *= two - if expo - y += log(BigDecimal("10"),prec) * BigDecimal(expo.to_s) - end - y - end - # Computes the value of pi to the specified number of digits of precision. def PI(prec) raise ArgumentError, "Zero or negative argument for PI" if prec <= 0 -- cgit v1.2.3