Class: Flt::DecNum

Inherits:

Object
Numeric
Num
Flt::DecNum

Includes:: AuxiliarFunctions

Defined in:: lib/flt/dec_num.rb,
lib/flt/trigonometry.rb

Overview

DecNum arbitrary precision floating point number. This implementation of DecNum is based on the Decimal module of Python, written by Eric Price, Facundo Batista, Raymond Hettinger, Aahz and Tim Peters.

Defined Under Namespace

Modules: AuxiliarFunctions, CMath, Math, Trigonometry Classes: Context

Constant Summary collapse

DefaultContext = the DefaultContext is the base for new contexts; it can be changed.

DecNum::Context.new(
:exact=>false, :precision=>28, :rounding=>:half_even,
:emin=> -999999999, :emax=>+999999999,
:flags=>[],
:traps=>[DivisionByZero, Overflow, InvalidOperation],
:ignored_flags=>[],
:capitals=>true,
:clamp=>true)

BasicContext =

DecNum::Context.new(DefaultContext,
:precision=>9, :rounding=>:half_up,
:traps=>[DivisionByZero, Overflow, InvalidOperation, Clamped, Underflow],
:flags=>[])

ExtendedContext =

DecNum::Context.new(DefaultContext,
:precision=>9, :rounding=>:half_even,
:traps=>[], :flags=>[], :clamp=>false)

Class Method Summary collapse

.int_div_radix_power(x, n) ⇒ Object

Divide by an integral power of the base: x/(radix**n) for x,n integer; returns an integer.
.int_mult_radix_power(x, n) ⇒ Object

Multiply by an integral power of the base: x*(radix**n) for x,n integer; returns an integer.
.int_radix_power(n) ⇒ Object

Integral power of the base: radix**n for integer n; returns an integer.
.radix ⇒ Object

Numerical base of DecNum.

Instance Method Summary collapse

#_ln_exp_bound ⇒ Object

Compute a lower bound for the adjusted exponent of self.ln().
#_log10_exp_bound ⇒ Object

Compute a lower bound for the adjusted exponent of self.log10() In other words, find r such that self.log10() >= 10**r.
#_power_exact(other, p) ⇒ Object

Attempt to compute self**other exactly Given Decimals self and other and an integer p, attempt to compute an exact result for the power self**other, with p digits of precision.
#_power_modulo(other, modulo, context = nil) ⇒ Object

Power-modulo: self._power_modulo(other, modulo) == (self**other) % modulo This is equivalent to Python’s 3-argument version of pow().
#exp(context = nil) ⇒ Object

Exponential function.
#initialize(*args) ⇒ DecNum constructor

A DecNum value can be defined by: * A String containing a text representation of the number * An Integer * A Rational * A Value of a type for which conversion is defined in the context.
#ln(context = nil) ⇒ Object

Returns the natural (base e) logarithm.
#log10(context = nil) ⇒ Object

Returns the base 10 logarithm.
#number_of_digits ⇒ Object
#power(other, modulo = nil, context = nil) ⇒ Object

Raises to the power of x, to modulo if given.

Methods included from AuxiliarFunctions

_dexp, _div_nearest, _dlog, _dlog10, _dpower, _iexp, _ilog, _log10_digits, _log10_lb, _number_of_digits, _rshift_nearest, _sqrt_nearest, dexp

Methods inherited from Num

#%, #*, #**, #+, #+@, #-, #-@, #/, #<, #<=, #<=>, #==, #>, #>=, Context, Flags, Num, [], #_abs, #_check_nans, #_fix, #_fix_nan, #_neg, #_pos, #_rescale, #_watched_rescale, #abs, #add, #adjusted_exponent, base_coercible_types, base_conversions, ccontext, #ceil, #coefficient, #coerce, #compare, context, context=, #convert_to, #copy_abs, #copy_negate, #copy_sign, define_context, #digits, #div, #divide, #divide_int, #divmod, #divrem, #eql?, #even?, #exponent, #finite?, #floor, #fma, #fraction_part, #fractional_exponent, #hash, #infinite?, infinity, #inspect, #integer_part, #integral?, #integral_exponent, #integral_significand, local_context, #log, #log2, #logb, math, #minus, #modulo, #multiply, nan, #nan?, #next_minus, #next_plus, #next_toward, #nonzero?, #normal?, #normalize, num_class, #num_class, #number_class, #odd?, one_half, #plus, #qnan?, #quantize, #rationalize, #reduce, #reduced_exponent, #remainder, #remainder_near, #rescale, #round, #same_quantum?, #scaleb, #scientific_exponent, set_context, #sign, #snan?, #special?, #split, #sqrt, #subnormal?, #subtract, #to_f, #to_i, #to_int_scale, #to_integral_exact, #to_integral_value, #to_r, #to_s, #truncate, #ulp, zero, #zero?

Methods included from Support

FlagValues, Flags, adjust_digits, simplified_round_mode

Methods included from Num::AuxiliarFunctions

_convert, _div_nearest, _exp, _iexp, _ilog, _log, _log_radix_digits, _log_radix_lb, _log_radix_mult, _normalize, _number_of_digits, _parser, _power, _rshift_nearest, _sqrt_nearest, log10_lb, log2_lb

Methods included from Support::AuxiliarFunctions

_nbits, _ndigits, detect_float_rounding

Methods inherited from Numeric

#even?, #odd?, #sign

Constructor Details

#initialize(*args) ⇒ `DecNum`

A DecNum value can be defined by:

A String containing a text representation of the number
An Integer
A Rational
A Value of a type for which conversion is defined in the context.
Another DecNum.
A sign, coefficient and exponent (either as separate arguments, as an array or as a Hash with symbolic keys), or a signed coefficient and an exponent. This is the internal representation of Num, as returned by Num#split. The sign is +1 for plus and -1 for minus; the coefficient and exponent are integers, except for special values which are defined by :inf, :nan or :snan for the exponent.

An optional Context can be passed after the value-definint argument to override the current context and options can be passed in a last hash argument; alternatively context options can be overriden by options of the hash argument.

When the number is defined by a numeric literal (a String), it can be followed by a symbol that specifies the mode used to convert the literal to a floating-point value:

:free is currently the default for all cases. The precision of the input literal (including trailing zeros) is preserved and the precision of the context is ignored. When the literal is in base 10, (which is the case by default), the literal is preserved exactly. Otherwise, all significative digits that can be derived from the literal are generanted, significative meaning here that if the digit is changed and the value converted back to a literal of the same base and precision, the original literal will not be obtained.
:short is a variation of :free in which only the minimun number of digits that are necessary to produce the original literal when the value is converted back with the same original precision.
:fixed will round and normalize the value to the precision specified by the context (normalize meaning that exaclty the number of digits specified by the precision will be generated, even if the original literal has fewer digits.) This may fail returning NaN (and raising Inexact) if the context precision is :exact, but not if the floating-point radix is a multiple of the input base.

Options that can be passed for construction from literal:

:base is the numeric base of the input, 10 by default.

The Flt.DecNum() constructor admits the same parameters and can be used as a shortcut for DecNum creation. Examples:

DecNum('0.1000')                                  # -> 0.1000
DecNum('0.12345')                                 # -> 0.12345
DecNum('1.2345E-1')                               # -> 0.12345
DecNum('0.1000', :short)                          # -> 0.1
DecNum('0.1000',:fixed, :precision=>20)           # -> 0.10000000000000000000
DecNum('0.12345',:fixed, :precision=>20)          # -> 0.12345000000000000000
DecNum('0.100110E3', :base=>2)                    # -> 4.8
DecNum('0.1E-5', :free, :base=>2)                 # -> 0.016
DecNum('0.1E-5', :short, :base=>2)                # -> 0.02
DecNum('0.1E-5', :fixed, :base=>2, :exact=>true)  # -> 0.015625
DecNum('0.1E-5', :fixed, :base=>2)                # -> 0.01562500000000000000000000000



115
116
117

# File 'lib/flt/dec_num.rb', line 115

def initialize(*args)
  super(*args)
end

Class Method Details

.int_div_radix_power(x, n) ⇒ `Object`

Divide by an integral power of the base: x/(radix**n) for x,n integer; returns an integer.



30
31
32

# File 'lib/flt/dec_num.rb', line 30

def int_div_radix_power(x,n)
  n < 0 ? (x * (10**(-n))) : (x / (10**n))
end

.int_mult_radix_power(x, n) ⇒ `Object`

Multiply by an integral power of the base: x*(radix**n) for x,n integer; returns an integer.



24
25
26

# File 'lib/flt/dec_num.rb', line 24

def int_mult_radix_power(x,n)
  n < 0 ? (x / (10**(-n))) : (x * (10**n))
end

.int_radix_power(n) ⇒ `Object`

Integral power of the base: radix**n for integer n; returns an integer.



18
19
20

# File 'lib/flt/dec_num.rb', line 18

def int_radix_power(n)
  10**n
end

.radix ⇒ `Object`

Numerical base of DecNum.



13
14
15

# File 'lib/flt/dec_num.rb', line 13

def radix
  10
end

Instance Method Details

#_ln_exp_bound ⇒ `Object`

Compute a lower bound for the adjusted exponent of self.ln(). In other words, compute r such that self.ln() >= 10**r. Assumes that self is finite and positive and that self != 1.

# File 'lib/flt/dec_num.rb', line 789

def _ln_exp_bound
  # for 0.1 <= x <= 10 we use the inequalities 1-1/x <= ln(x) <= x-1
  #
  # The original Python cod used lexical order (having converted to strings) for (num < den))
  # so the results would be different e.g. for num = 9m den=200; Can this happen? What is the correct way?

  adj = self.exponent + number_of_digits - 1
  if adj >= 1
    # argument >= 10; we use 23/10 = 2.3 as a lower bound for ln(10)
    return _number_of_digits(adj*23/10) - 1
  end
  if adj <= -2
    # argument <= 0.1
    return _number_of_digits((-1-adj)*23/10) - 1
  end
  c = self.coefficient
  e = self.exponent
  if adj == 0
    # 1 < self < 10
    num = c-(10**-e)
    den = c
    return _number_of_digits(num) - _number_of_digits(den) - ((num < den) ? 1 : 0)
  end
  # adj == -1, 0.1 <= self < 1
  return e + _number_of_digits(10**-e - c) - 1
end

#_log10_exp_bound ⇒ `Object`

Compute a lower bound for the adjusted exponent of self.log10() In other words, find r such that self.log10() >= 10**r. Assumes that self is finite and positive and that self != 1.

# File 'lib/flt/dec_num.rb', line 759

def _log10_exp_bound
  # For x >= 10 or x < 0.1 we only need a bound on the integer
  # part of log10(self), and this comes directly from the
  # exponent of x.  For 0.1 <= x <= 10 we use the inequalities
  # 1-1/x <= log(x) <= x-1. If x > 1 we have |log10(x)| >
  # (1-1/x)/2.31 > 0.  If x < 1 then |log10(x)| > (1-x)/2.31 > 0
  #
  # The original Python cod used lexical order (having converted to strings) for (num < den) and (num < 231)
  # so the results would be different e.g. for num = 9; Can this happen? What is the correct way?

  adj = self.exponent + number_of_digits - 1
  return _number_of_digits(adj) - 1 if adj >= 1 # self >= 10
  return _number_of_digits(-1-adj)-1 if adj <= -2 # self < 0.1

  c = self.coefficient
  e = self.exponent
  if adj == 0
    # 1 < self < 10
    num = (c - DecNum.int_radix_power(-e))
    den = (231*c)
    return _number_of_digits(num) - _number_of_digits(den) - ((num < den) ? 1 : 0) + 2
  end
  # adj == -1, 0.1 <= self < 1
  num = (DecNum.int_radix_power(-e)-c)
  return _number_of_digits(num.to_i) + e - ((num < 231) ? 1 : 0) - 1
end

#_power_exact(other, p) ⇒ `Object`

Attempt to compute self**other exactly Given Decimals self and other and an integer p, attempt to compute an exact result for the power self**other, with p digits of precision. Return nil if self**other is not exactly representable in p digits.

Assumes that elimination of special cases has already been performed: self and other must both be nonspecial; self must be positive and not numerically equal to 1; other must be nonzero. For efficiency, other.exponent should not be too large, so that 10**other.exponent.abs is a feasible calculation.

# File 'lib/flt/dec_num.rb', line 558

def _power_exact(other, p)

  # In the comments below, we write x for the value of self and
  # y for the value of other.  Write x = xc*10**xe and y =
  # yc*10**ye.

  # The main purpose of this method is to identify the *failure*
  # of x**y to be exactly representable with as little effort as
  # possible.  So we look for cheap and easy tests that
  # eliminate the possibility of x**y being exact.  Only if all
  # these tests are passed do we go on to actually compute x**y.

  # Here's the main idea.  First normalize both x and y.  We
  # express y as a rational m/n, with m and n relatively prime
  # and n>0.  Then for x**y to be exactly representable (at
  # *any* precision), xc must be the nth power of a positive
  # integer and xe must be divisible by n.  If m is negative
  # then additionally xc must be a power of either 2 or 5, hence
  # a power of 2**n or 5**n.
  #
  # There's a limit to how small |y| can be: if y=m/n as above
  # then:
  #
  #  (1) if xc != 1 then for the result to be representable we
  #      need xc**(1/n) >= 2, and hence also xc**|y| >= 2.  So
  #      if |y| <= 1/nbits(xc) then xc < 2**nbits(xc) <=
  #      2**(1/|y|), hence xc**|y| < 2 and the result is not
  #      representable.
  #
  #  (2) if xe != 0, |xe|*(1/n) >= 1, so |xe|*|y| >= 1.  Hence if
  #      |y| < 1/|xe| then the result is not representable.
  #
  # Note that since x is not equal to 1, at least one of (1) and
  # (2) must apply.  Now |y| < 1/nbits(xc) iff |yc|*nbits(xc) <
  # 10**-ye iff len(str(|yc|*nbits(xc)) <= -ye.
  #
  # There's also a limit to how large y can be, at least if it's
  # positive: the normalized result will have coefficient xc**y,
  # so if it's representable then xc**y < 10**p, and y <
  # p/log10(xc).  Hence if y*log10(xc) >= p then the result is
  # not exactly representable.

  # if len(str(abs(yc*xe)) <= -ye then abs(yc*xe) < 10**-ye,
  # so |y| < 1/xe and the result is not representable.
  # Similarly, len(str(abs(yc)*xc_bits)) <= -ye implies |y|
  # < 1/nbits(xc).

  xc = self.coefficient
  xe = self.exponent
  while (xc % DecNum.radix) == 0
    xc /= DecNum.radix
    xe += 1
  end

  yc = other.coefficient
  ye = other.exponent
  while (yc % DecNum.radix) == 0
    yc /= DecNum.radix
    ye += 1
  end

  # case where xc == 1: result is 10**(xe*y), with xe*y
  # required to be an integer
  if xc == 1
    if ye >= 0
      exponent = xe*yc*DecNum.int_radix_power(ye)
    else
      exponent, remainder = (xe*yc).divmod(DecNum.int_radix_power(-ye))
      return nil if remainder!=0
    end
    exponent = -exponent if other.sign == -1
    # if other is a nonnegative integer, use ideal exponent
    if other.integral? and (other.sign == +1)
      ideal_exponent = self.exponent*other.to_i
      zeros = [exponent-ideal_exponent, p-1].min
    else
      zeros = 0
    end
    return Num(+1, DecNum.int_radix_power(zeros), exponent-zeros)
  end

  # case where y is negative: xc must be either a power
  # of 2 or a power of 5.
  if other.sign == -1
    last_digit = (xc % 10)
    if [2,4,6,8].include?(last_digit)
      # quick test for power of 2
      return nil if xc & -xc != xc
      # now xc is a power of 2; e is its exponent
      e = _nbits(xc)-1
      # find e*y and xe*y; both must be integers
      if ye >= 0
        y_as_int = yc*DecNum.int_radix_power(ye)
        e = e*y_as_int
        xe = xe*y_as_int
      else
        ten_pow = DecNum.int_radix_power(-ye)
        e, remainder = (e*yc).divmod(ten_pow)
        return nil if remainder!=0
        xe, remainder = (xe*yc).divmod(ten_pow)
        return nil if remainder!=0
      end

      return nil if e*65 >= p*93 # 93/65 > log(10)/log(5)
      xc = 5**e
    elsif last_digit == 5
      # e >= log_5(xc) if xc is a power of 5; we have
      # equality all the way up to xc=5**2658
      e = _nbits(xc)*28/65
      xc, remainder = (5**e).divmod(xc)
      return nil if remainder!=0
      while (xc % 5) == 0
        xc /= 5
        e -= 1
      end
      if ye >= 0
        y_as_integer = DecNum.int_mult_radix_power(yc,ye)
        e = e*y_as_integer
        xe = xe*y_as_integer
      else
        ten_pow = DecNum.int_radix_power(-ye)
        e, remainder = (e*yc).divmod(ten_pow)
        return nil if remainder
        xe, remainder = (xe*yc).divmod(ten_pow)
        return nil if remainder
      end
      return nil if e*3 >= p*10 # 10/3 > log(10)/log(2)
      xc = 2**e
    else
      return nil
    end

    return nil if xc >= DecNum.int_radix_power(p)
    xe = -e-xe
    return Num(+1, xc, xe)

  end

  # now y is positive; find m and n such that y = m/n
  if ye >= 0
    m, n = yc*10**ye, 1
  else
    return nil if (xe != 0) and (_number_of_digits((yc*xe).abs) <= -ye)
    xc_bits = _nbits(xc)
    return nil if (xc != 1) and (_number_of_digits(yc.abs*xc_bits) <= -ye)
    m, n = yc, DecNum.int_radix_power(-ye)
    while ((m % 2) == 0) && ((n % 2) == 0)
      m /= 2
      n /= 2
    end
    while ((m % 5) == 0) && ((n % 5) == 0)
      m /= 5
      n /= 5
    end
  end

  # compute nth root of xc*10**xe
  if n > 1
    # if 1 < xc < 2**n then xc isn't an nth power
    return nil if xc != 1 and xc_bits <= n

    xe, rem = xe.divmod(n)
    return nil if rem != 0

    # compute nth root of xc using Newton's method
    a = 1 << -(-_nbits(xc)/n) # initial estimate
    q = r = nil
    loop do
      q, r = xc.divmod(a**(n-1))
      break if a <= q
      a = (a*(n-1) + q)/n
    end
    return nil if !((a == q) and (r == 0))
    xc = a
  end

  # now xc*10**xe is the nth root of the original xc*10**xe
  # compute mth power of xc*10**xe

  # if m > p*100/_log10_lb(xc) then m > p/log10(xc), hence xc**m >
  # 10**p and the result is not representable.
  return nil if (xc > 1) and (m > p*100/_log10_lb(xc))
  xc = xc**m
  xe *= m
  return nil if xc > 10**p

  # by this point the result *is* exactly representable
  # adjust the exponent to get as close as possible to the ideal
  # exponent, if necessary
  if other.integral? && other.sign == +1
    ideal_exponent = self.exponent*other.to_i
    zeros = [xe-ideal_exponent, p-_number_of_digits(xc)].min
  else
    zeros = 0
  end
  return Num(+1, DecNum.int_mult_radix_power(xc, zeros), xe-zeros)
end

#_power_modulo(other, modulo, context = nil) ⇒ `Object`

Power-modulo: self._power_modulo(other, modulo) == (self**other) % modulo This is equivalent to Python’s 3-argument version of pow()

# File 'lib/flt/dec_num.rb', line 499

def _power_modulo(other, modulo, context=nil)

  context = DecNum.define_context(context)
  other = _convert(other)
  modulo = _convert(third)

  if self.nan? || other.nan? || modulo.nan?
    return context.exception(InvalidOperation, 'sNaN', self) if self.snan?
    return context.exception(InvalidOperation, 'sNaN', other) if other.snan?
    return context.exception(InvalidOperation, 'sNaN', modulo) if other.modulo?
    return self._fix_nan(context) if self.nan?
    return other._fix_nan(context) if other.nan?
    return modulo._fix_nan(context) # if modulo.nan?
  end

  if !(self.integral? && other.integral? && modulo.integral?)
    return context.exception(InvalidOperation, '3-argument power not allowed unless all arguments are integers.')
  end

  if other < 0
    return context.exception(InvalidOperation, '3-argument power cannot have a negative 2nd argument.')
  end

  if modulo.zero?
    return context.exception(InvalidOperation, '3-argument power cannot have a 0 3rd argument.')
  end

  if modulo.adjusted_exponent >= context.precision
    return context.exception(InvalidOperation, 'insufficient precision: power 3rd argument must not have more than precision digits')
  end

  if other.zero? && self.zero?
    return context.exception(InvalidOperation, "0**0 not defined")
  end

  sign = other.even? ? +1 : -1
  modulo = modulo.to_i.abs

  base = (self.coefficient % modulo * (DecNum.int_radix_power(self.exponent) % modulo)) % modulo

  other.exponent.times do
    base = (base**DecNum.radix) % modulo
  end
  base = (base**other.coefficient) % modulo

  Num(sign, base, 0)
end

#exp(context = nil) ⇒ `Object`

Exponential function

# File 'lib/flt/dec_num.rb', line 377

def exp(context=nil)
  context = DecNum.define_context(context)

  # exp(NaN) = NaN
  ans = _check_nans(context)
  return ans if ans

  # exp(-Infinity) = 0
  return DecNum.zero if self.infinite? && (self.sign == -1)

  # exp(0) = 1
  return Num(1) if self.zero?

  # exp(Infinity) = Infinity
  return Num(self) if self.infinite?

  # the result is now guaranteed to be inexact (the true
  # mathematical result is transcendental). There's no need to
  # raise Rounded and Inexact here---they'll always be raised as
  # a result of the call to _fix.
  return context.exception(Inexact, 'Inexact exp') if context.exact?
  p = context.precision
  adj = self.adjusted_exponent

  # we only need to do any computation for quite a small range
  # of adjusted exponents---for example, -29 <= adj <= 10 for
  # the default context.  For smaller exponent the result is
  # indistinguishable from 1 at the given precision, while for
  # larger exponent the result either overflows or underflows.
  if self.sign == +1 and adj > _number_of_digits((context.emax+1)*3)
    # overflow
    ans = Num(+1, 1, context.emax+1)
  elsif self.sign == -1 and adj > _number_of_digits((-context.etiny+1)*3)
    # underflow to 0
    ans = Num(+1, 1, context.etiny-1)
  elsif self.sign == +1 and adj < -p
    # p+1 digits; final round will raise correct flags
    ans = Num(+1, DecNum.int_radix_power(p)+1, -p)
  elsif self.sign == -1 and adj < -p-1
    # p+1 digits; final round will raise correct flags
    ans = Num(+1, DecNum.int_radix_power(p+1)-1, -p-1)
  else
    # general case
    c = self.coefficient
    e = self.exponent
    c = -c if self.sign == -1

    # compute correctly rounded result: increase precision by
    # 3 digits at a time until we get an unambiguously
    # roundable result
    extra = 3
    coeff = exp = nil
    loop do
      coeff, exp = _dexp(c, e, p+extra)
      break if (coeff % (5*10**(_number_of_digits(coeff)-p-1)))!=0
      extra += 3
    end
    ans = Num(+1, coeff, exp)
  end

  # at this stage, ans should round correctly with *any*
  # rounding mode, not just with ROUND_HALF_EVEN
  DecNum.context(context, :rounding=>:half_even) do |local_context|
    ans = ans._fix(local_context)
    context.flags = local_context.flags
  end

  return ans
end

#ln(context = nil) ⇒ `Object`

Returns the natural (base e) logarithm

# File 'lib/flt/dec_num.rb', line 448

def ln(context=nil)
  context = DecNum.define_context(context)

  # ln(NaN) = NaN
  ans = _check_nans(context)
  return ans if ans

  # ln(0.0) == -Infinity
  return DecNum.infinity(-1) if self.zero?

  # ln(Infinity) = Infinity
  return DecNum.infinity if self.infinite? && self.sign == +1

  # ln(1.0) == 0.0
  return DecNum.zero if self == Num(1)

  # ln(negative) raises InvalidOperation
  return context.exception(InvalidOperation, 'ln of a negative value') if self.sign==-1

  # result is irrational, so necessarily inexact
  return context.exception(Inexact, 'Inexact exp') if context.exact?

  c = self.coefficient
  e = self.exponent
  p = context.precision

  # correctly rounded result: repeatedly increase precision by 3
  # until we get an unambiguously roundable result
  places = p - self._ln_exp_bound + 2 # at least p+3 places
  coeff = nil
  loop do
    coeff = _dlog(c, e, places)
    # assert coeff.to_s.length-p >= 1
    break if (coeff % (5*10**(_number_of_digits(coeff.abs)-p-1))) != 0
    places += 3
  end
  ans = Num((coeff<0) ? -1 : +1, coeff.abs, -places)

  DecNum.context(context, :rounding=>:half_even) do |local_context|
    ans = ans._fix(local_context)
    context.flags = local_context.flags
  end
  return ans
end

#log10(context = nil) ⇒ `Object`

Returns the base 10 logarithm

# File 'lib/flt/dec_num.rb', line 327

def log10(context=nil)
  context = DecNum.define_context(context)

  # log10(NaN) = NaN
  ans = _check_nans(context)
  return ans if ans

  # log10(0.0) == -Infinity
  return DecNum.infinity(-1) if self.zero?

  # log10(Infinity) = Infinity
  return DecNum.infinity if self.infinite? && self.sign == +1

  # log10(negative or -Infinity) raises InvalidOperation
  return context.exception(InvalidOperation, 'log10 of a negative value') if self.sign == -1

  digits = self.digits
  # log10(10**n) = n
  if digits.first == 1 && digits[1..-1].all?{|d| d==0}
    # answer may need rounding
    ans = Num(self.exponent + digits.size - 1)
    return ans if context.exact?
  else
    # result is irrational, so necessarily inexact
    return context.exception(Inexact, "Inexact power") if context.exact?
    c = self.coefficient
    e = self.exponent
    p = context.precision

    # correctly rounded result: repeatedly increase precision
    # until result is unambiguously roundable
    places = p-self._log10_exp_bound+2
    coeff = nil
    loop do
      coeff = _dlog10(c, e, places)
      # assert coeff.abs.to_s.length-p >= 1
      break if (coeff % (5*10**(_number_of_digits(coeff.abs)-p-1)))!=0
      places += 3
    end
    ans = Num(coeff<0 ? -1 : +1, coeff.abs, -places)
  end

  DecNum.context(context, :rounding=>:half_even) do |local_context|
    ans = ans._fix(local_context)
    context.flags = local_context.flags
  end
  return ans
end

#number_of_digits ⇒ `Object`



119
120
121

# File 'lib/flt/dec_num.rb', line 119

def number_of_digits
  @coeff.is_a?(Integer) ? _number_of_digits(@coeff) : 0
end

#power(other, modulo = nil, context = nil) ⇒ `Object`

Raises to the power of x, to modulo if given.

With two arguments, compute self**other. If self is negative then other must be integral. The result will be inexact unless other is integral and the result is finite and can be expressed exactly in ‘precision’ digits.

With three arguments, compute (self**other) % modulo. For the three argument form, the following restrictions on the arguments hold:

- all three arguments must be integral
- other must be nonnegative
- at least one of self or other must be nonzero
- modulo must be nonzero and have at most 'precision' digits

The result of a.power(b, modulo) is identical to the result that would be obtained by computing (a**b) % modulo with unbounded precision, but is computed more efficiently. It is always exact.

# File 'lib/flt/dec_num.rb', line 143

def power(other, modulo=nil, context=nil)

  if context.nil? && (modulo.is_a?(Context) || modulo.is_a?(Hash))
    context = modulo
    modulo = nil
  end

  return self.power_modulo(other, modulo, context) if modulo

  context = DecNum.define_context(context)
  other = _convert(other)

  ans = _check_nans(context, other)
  return ans if ans

  # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
  if other.zero?
    if self.zero?
      return context.exception(InvalidOperation, '0 ** 0')
    else
      return Num(1)
    end
  end

  # result has sign -1 iff self.sign is -1 and other is an odd integer
  result_sign = +1
  _self = self
  if _self.sign == -1
    if other.integral?
      result_sign = -1 if !other.even?
    else
      # -ve**noninteger = NaN
      # (-0)**noninteger = 0**noninteger
      unless self.zero?
        return context.exception(InvalidOperation, 'x ** y with x negative and y not an integer')
      end
    end
    # negate self, without doing any unwanted rounding
    _self = self.copy_negate
  end

  # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
  if _self.zero?
    return (other.sign == +1) ? Num(result_sign, 0, 0) : num_class.infinity(result_sign)
  end

  # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
  if _self.infinite?
    return (other.sign == +1) ? num_class.infinity(result_sign) : Num(result_sign, 0, 0)
  end

  # 1**other = 1, but the choice of exponent and the flags
  # depend on the exponent of self, and on whether other is a
  # positive integer, a negative integer, or neither
  if _self == Num(1)
    return _self if context.exact?
    if other.integral?
      # exp = max(self._exp*max(int(other), 0),
      # 1-context.prec) but evaluating int(other) directly
      # is dangerous until we know other is small (other
      # could be 1e999999999)
      if other.sign == -1
        multiplier = 0
      elsif other > context.precision
        multiplier = context.precision
      else
        multiplier = other.to_i
      end

      exp = _self.exponent * multiplier
      if exp < 1-context.precision
        exp = 1-context.precision
        context.exception Rounded
      end
    else
      context.exception Rounded
      context.exception Inexact
      exp = 1-context.precision
    end

    return Num(result_sign, DecNum.int_radix_power(-exp), exp)
  end

  # compute adjusted exponent of self
  self_adj = _self.adjusted_exponent

  # self ** infinity is infinity if self > 1, 0 if self < 1
  # self ** -infinity is infinity if self < 1, 0 if self > 1
  if other.infinite?
    if (other.sign == +1) == (self_adj < 0)
      return Num(result_sign, 0, 0)
    else
      return DecNum.infinity(result_sign)
    end
  end

  # from here on, the result always goes through the call
  # to _fix at the end of this function.
  ans = nil

  # crude test to catch cases of extreme overflow/underflow.  If
  # log10(self)*other >= 10**bound and bound >= len(str(Emax))
  # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
  # self**other >= 10**(Emax+1), so overflow occurs.  The test
  # for underflow is similar.
  bound = _self._log10_exp_bound + other.adjusted_exponent
  if (self_adj >= 0) == (other.sign == +1)
    # self > 1 and other +ve, or self < 1 and other -ve
    # possibility of overflow
    if bound >= _number_of_digits(context.emax)
      ans = Num(result_sign, 1, context.emax+1)
    end
  else
    # self > 1 and other -ve, or self < 1 and other +ve
    # possibility of underflow to 0
    etiny = context.etiny
    if bound >= _number_of_digits(-etiny)
      ans = Num(result_sign, 1, etiny-1)
    end
  end

  # try for an exact result with precision +1
  if ans.nil?
    if context.exact?
      if other.adjusted_exponent < 100
        test_precision = _self.number_of_digits*other.to_i+1
      else
        test_precision = _self.number_of_digits+1
      end
    else
      test_precision = context.precision + 1
    end
    ans = _self._power_exact(other, test_precision)
    if !ans.nil? && (result_sign == -1)
      ans = Num(-1, ans.coefficient, ans.exponent)
    end
  end

  # usual case: inexact result, x**y computed directly as exp(y*log(x))
  if !ans.nil?
    return ans if context.exact?
  else
    return context.exception(Inexact, "Inexact power") if context.exact?

    p = context.precision
    xc = _self.coefficient
    xe = _self.exponent
    yc = other.coefficient
    ye = other.exponent
    yc = -yc if other.sign == -1

    # compute correctly rounded result:  start with precision +3,
    # then increase precision until result is unambiguously roundable
    extra = 3
    coeff, exp = nil, nil
    loop do
      coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
      #break if (coeff % DecNum.int_mult_radix_power(5,coeff.to_s.length-p-1)) != 0
      break if (coeff % (5*10**(_number_of_digits(coeff)-p-1))) != 0
      extra += 3
    end
    ans = Num(result_sign, coeff, exp)
  end

  # the specification says that for non-integer other we need to
  # raise Inexact, even when the result is actually exact.  In
  # the same way, we need to raise Underflow here if the result
  # is subnormal.  (The call to _fix will take care of raising
  # Rounded and Subnormal, as usual.)
  if !other.integral?
    context.exception Inexact
    # pad with zeros up to length context.precision+1 if necessary
    if ans.number_of_digits <= context.precision
      expdiff = context.precision+1 - ans.number_of_digits
      ans = Num(ans.sign, DecNum.int_mult_radix_power(ans.coefficient, expdiff), ans.exponent-expdiff)
    end
    context.exception Underflow if ans.adjusted_exponent < context.emin
  end
  # unlike exp, ln and log10, the power function respects the
  # rounding mode; no need to use ROUND_HALF_EVEN here
  ans._fix(context)
end

Class: Flt::DecNum

Overview

Defined Under Namespace

Constant Summary collapse

Constants included from AuxiliarFunctions

Constants inherited from Num

Constants included from Num::AuxiliarFunctions

Constants included from Support::AuxiliarFunctions

Class Method Summary collapse

Instance Method Summary collapse

Methods included from AuxiliarFunctions

Methods inherited from Num

Methods included from Support

Methods included from Num::AuxiliarFunctions

Methods included from Support::AuxiliarFunctions

Methods inherited from Numeric

Constructor Details

#initialize(*args) ⇒ DecNum

Class Method Details

.int_div_radix_power(x, n) ⇒ Object

.int_mult_radix_power(x, n) ⇒ Object

.int_radix_power(n) ⇒ Object

.radix ⇒ Object

Instance Method Details

#_ln_exp_bound ⇒ Object

#_log10_exp_bound ⇒ Object

#_power_exact(other, p) ⇒ Object

#_power_modulo(other, modulo, context = nil) ⇒ Object

#exp(context = nil) ⇒ Object

#ln(context = nil) ⇒ Object

#log10(context = nil) ⇒ Object

#number_of_digits ⇒ Object

#power(other, modulo = nil, context = nil) ⇒ Object