Class: Flt::DecNum

Inherits:
Num show all
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

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)

Constants included from AuxiliarFunctions

AuxiliarFunctions::LOG10_LB_CORRECTION, AuxiliarFunctions::NUMBER_OF_DIGITS_MAX_VALID_LOG

Constants inherited from Num

Num::EXCEPTIONS, Num::ROUND_05UP, Num::ROUND_CEILING, Num::ROUND_DOWN, Num::ROUND_FLOOR, Num::ROUND_HALF_DOWN, Num::ROUND_HALF_EVEN, Num::ROUND_HALF_UP, Num::ROUND_UP

Constants included from Num::AuxiliarFunctions

Num::AuxiliarFunctions::EXP_INC, Num::AuxiliarFunctions::LOG10_LB_CORRECTION, Num::AuxiliarFunctions::LOG10_MULT, Num::AuxiliarFunctions::LOG2_LB_CORRECTION, Num::AuxiliarFunctions::LOG2_MULT, Num::AuxiliarFunctions::LOG_PREC_INC, Num::AuxiliarFunctions::LOG_RADIX_EXTRA, Num::AuxiliarFunctions::LOG_RADIX_INC

Constants included from Support::AuxiliarFunctions

Support::AuxiliarFunctions::NBITS_BLOCK, Support::AuxiliarFunctions::NBITS_LIMIT, Support::AuxiliarFunctions::NDIGITS_BLOCK, Support::AuxiliarFunctions::NDIGITS_LIMIT

Class Method Summary (collapse)

Instance Method Summary (collapse)

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, #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 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

#FlagValues, #Flags, #simplified_round_mode

Methods included from Support::AuxiliarFunctions

#_nbits, #_ndigits, #detect_float_rounding

Methods inherited from Numeric

#even?, #odd?, #sign

Constructor Details

- (DecNum) initialize(*args)

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.1000
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

+ (Object) int_div_radix_power(x, n)

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

+ (Object) int_mult_radix_power(x, n)

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

+ (Object) int_radix_power(n)

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

+ (Object) radix

Numerical base of DecNum.



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

def radix
  10
end

Instance Method Details

- (Object) _ln_exp_bound

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.



789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
# 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

- (Object) _log10_exp_bound

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.



759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
# 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

- (Object) _power_exact(other, p)

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.



558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
# 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

- (Object) _power_modulo(other, modulo, context = nil)

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



499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
# 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

- (Object) exp(context = nil)

Exponential function



377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
# 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

- (Object) ln(context = nil)

Returns the natural (base e) logarithm



448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
# 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

- (Object) log10(context = nil)

Returns the base 10 logarithm



327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
# 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

- (Object) number_of_digits



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

- (Object) power(other, modulo = nil, context = nil)

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.



143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
# 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