Parent

Included Modules

Class Index [+]

Quicksearch

Flt::DecNum

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.

Constants

DefaultContext

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

BasicContext
ExtendedContext

Public Class Methods

int_div_radix_power(x,n) click to toggle source

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

    # File lib/flt/dec_num.rb, line 30
30:     def int_div_radix_power(x,n)
31:       n < 0 ? (x * (10**(-n))) : (x / (10**n))
32:     end
int_mult_radix_power(x,n) click to toggle source

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

    # File lib/flt/dec_num.rb, line 24
24:     def int_mult_radix_power(x,n)
25:       n < 0 ? (x / (10**(-n))) : (x * (10**n))
26:     end
int_radix_power(n) click to toggle source

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

    # File lib/flt/dec_num.rb, line 18
18:     def int_radix_power(n)
19:       10**n
20:     end
new(*args) click to toggle source

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
     # File lib/flt/dec_num.rb, line 115
115:   def initialize(*args)
116:     super(*args)
117:   end
radix() click to toggle source

Numerical base of DecNum.

    # File lib/flt/dec_num.rb, line 13
13:     def radix
14:       10
15:     end

Public Instance Methods

_ln_exp_bound() click to toggle source

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
789:   def _ln_exp_bound
790:     # for 0.1 <= x <= 10 we use the inequalities 1-1/x <= ln(x) <= x-1
791:     #
792:     # The original Python cod used lexical order (having converted to strings) for (num < den))
793:     # so the results would be different e.g. for num = 9m den=200; Can this happen? What is the correct way?
794: 
795:     adj = self.exponent + number_of_digits - 1
796:     if adj >= 1
797:       # argument >= 10; we use 23/10 = 2.3 as a lower bound for ln(10)
798:       return _number_of_digits(adj*23/10) - 1
799:     end
800:     if adj <= 2
801:       # argument <= 0.1
802:       return _number_of_digits((1-adj)*23/10) - 1
803:     end
804:     c = self.coefficient
805:     e = self.exponent
806:     if adj == 0
807:       # 1 < self < 10
808:       num = c-(10**-e)
809:       den = c
810:       return _number_of_digits(num) - _number_of_digits(den) - ((num < den) ? 1 : 0)
811:     end
812:     # adj == -1, 0.1 <= self < 1
813:     return e + _number_of_digits(10**-e - c) - 1
814:   end
_log10_exp_bound() click to toggle source

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
759:   def _log10_exp_bound
760:     # For x >= 10 or x < 0.1 we only need a bound on the integer
761:     # part of log10(self), and this comes directly from the
762:     # exponent of x.  For 0.1 <= x <= 10 we use the inequalities
763:     # 1-1/x <= log(x) <= x-1. If x > 1 we have |log10(x)| >
764:     # (1-1/x)/2.31 > 0.  If x < 1 then |log10(x)| > (1-x)/2.31 > 0
765:     #
766:     # The original Python cod used lexical order (having converted to strings) for (num < den) and (num < 231)
767:     # so the results would be different e.g. for num = 9; Can this happen? What is the correct way?
768: 
769:     adj = self.exponent + number_of_digits - 1
770:     return _number_of_digits(adj) - 1 if adj >= 1 # self >= 10
771:     return _number_of_digits(1-adj)-1 if adj <= 2 # self < 0.1
772: 
773:     c = self.coefficient
774:     e = self.exponent
775:     if adj == 0
776:       # 1 < self < 10
777:       num = (c - DecNum.int_radix_power(-e))
778:       den = (231*c)
779:       return _number_of_digits(num) - _number_of_digits(den) - ((num < den) ? 1 : 0) + 2
780:     end
781:     # adj == -1, 0.1 <= self < 1
782:     num = (DecNum.int_radix_power(-e)-c)
783:     return _number_of_digits(num.to_i) + e - ((num < 231) ? 1 : 0) - 1
784:   end
_power_exact(other, p) click to toggle source

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
558:   def _power_exact(other, p)
559: 
560:     # In the comments below, we write x for the value of self and
561:     # y for the value of other.  Write x = xc*10**xe and y =
562:     # yc*10**ye.
563: 
564:     # The main purpose of this method is to identify the *failure*
565:     # of x**y to be exactly representable with as little effort as
566:     # possible.  So we look for cheap and easy tests that
567:     # eliminate the possibility of x**y being exact.  Only if all
568:     # these tests are passed do we go on to actually compute x**y.
569: 
570:     # Here's the main idea.  First normalize both x and y.  We
571:     # express y as a rational m/n, with m and n relatively prime
572:     # and n>0.  Then for x**y to be exactly representable (at
573:     # *any* precision), xc must be the nth power of a positive
574:     # integer and xe must be divisible by n.  If m is negative
575:     # then additionally xc must be a power of either 2 or 5, hence
576:     # a power of 2**n or 5**n.
577:     #
578:     # There's a limit to how small |y| can be: if y=m/n as above
579:     # then:
580:     #
581:     #  (1) if xc != 1 then for the result to be representable we
582:     #      need xc**(1/n) >= 2, and hence also xc**|y| >= 2.  So
583:     #      if |y| <= 1/nbits(xc) then xc < 2**nbits(xc) <=
584:     #      2**(1/|y|), hence xc**|y| < 2 and the result is not
585:     #      representable.
586:     #
587:     #  (2) if xe != 0, |xe|*(1/n) >= 1, so |xe|*|y| >= 1.  Hence if
588:     #      |y| < 1/|xe| then the result is not representable.
589:     #
590:     # Note that since x is not equal to 1, at least one of (1) and
591:     # (2) must apply.  Now |y| < 1/nbits(xc) iff |yc|*nbits(xc) <
592:     # 10**-ye iff len(str(|yc|*nbits(xc)) <= -ye.
593:     #
594:     # There's also a limit to how large y can be, at least if it's
595:     # positive: the normalized result will have coefficient xc**y,
596:     # so if it's representable then xc**y < 10**p, and y <
597:     # p/log10(xc).  Hence if y*log10(xc) >= p then the result is
598:     # not exactly representable.
599: 
600:     # if len(str(abs(yc*xe)) <= -ye then abs(yc*xe) < 10**-ye,
601:     # so |y| < 1/xe and the result is not representable.
602:     # Similarly, len(str(abs(yc)*xc_bits)) <= -ye implies |y|
603:     # < 1/nbits(xc).
604: 
605:     xc = self.coefficient
606:     xe = self.exponent
607:     while (xc % DecNum.radix) == 0
608:       xc /= DecNum.radix
609:       xe += 1
610:     end
611: 
612:     yc = other.coefficient
613:     ye = other.exponent
614:     while (yc % DecNum.radix) == 0
615:       yc /= DecNum.radix
616:       ye += 1
617:     end
618: 
619:     # case where xc == 1: result is 10**(xe*y), with xe*y
620:     # required to be an integer
621:     if xc == 1
622:       if ye >= 0
623:         exponent = xe*yc*DecNum.int_radix_power(ye)
624:       else
625:         exponent, remainder = (xe*yc).divmod(DecNum.int_radix_power(-ye))
626:         return nil if remainder!=0
627:       end
628:       exponent = -exponent if other.sign == 1
629:       # if other is a nonnegative integer, use ideal exponent
630:       if other.integral? and (other.sign == 1)
631:         ideal_exponent = self.exponent*other.to_i
632:         zeros = [exponent-ideal_exponent, p-1].min
633:       else
634:         zeros = 0
635:       end
636:       return Num(1, DecNum.int_radix_power(zeros), exponent-zeros)
637:     end
638: 
639:     # case where y is negative: xc must be either a power
640:     # of 2 or a power of 5.
641:     if other.sign == 1
642:       last_digit = (xc % 10)
643:       if [2,4,6,8].include?(last_digit)
644:         # quick test for power of 2
645:         return nil if xc & -xc != xc
646:         # now xc is a power of 2; e is its exponent
647:         e = _nbits(xc)-1
648:         # find e*y and xe*y; both must be integers
649:         if ye >= 0
650:           y_as_int = yc*DecNum.int_radix_power(ye)
651:           e = e*y_as_int
652:           xe = xe*y_as_int
653:         else
654:           ten_pow = DecNum.int_radix_power(-ye)
655:           e, remainder = (e*yc).divmod(ten_pow)
656:           return nil if remainder!=0
657:           xe, remainder = (xe*yc).divmod(ten_pow)
658:           return nil if remainder!=0
659:         end
660: 
661:         return nil if e*65 >= p*93 # 93/65 > log(10)/log(5)
662:         xc = 5**e
663:       elsif last_digit == 5
664:         # e >= log_5(xc) if xc is a power of 5; we have
665:         # equality all the way up to xc=5**2658
666:         e = _nbits(xc)*28/65
667:         xc, remainder = (5**e).divmod(xc)
668:         return nil if remainder!=0
669:         while (xc % 5) == 0
670:           xc /= 5
671:           e -= 1
672:         end
673:         if ye >= 0
674:           y_as_integer = DecNum.int_mult_radix_power(yc,ye)
675:           e = e*y_as_integer
676:           xe = xe*y_as_integer
677:         else
678:           ten_pow = DecNum.int_radix_power(-ye)
679:           e, remainder = (e*yc).divmod(ten_pow)
680:           return nil if remainder
681:           xe, remainder = (xe*yc).divmod(ten_pow)
682:           return nil if remainder
683:         end
684:         return nil if e*3 >= p*10 # 10/3 > log(10)/log(2)
685:         xc = 2**e
686:       else
687:         return nil
688:       end
689: 
690:       return nil if xc >= DecNum.int_radix_power(p)
691:       xe = -e-xe
692:       return Num(1, xc, xe)
693: 
694:     end
695: 
696:     # now y is positive; find m and n such that y = m/n
697:     if ye >= 0
698:       m, n = yc*10**ye, 1
699:     else
700:       return nil if (xe != 0) and (_number_of_digits((yc*xe).abs) <= -ye)
701:       xc_bits = _nbits(xc)
702:       return nil if (xc != 1) and (_number_of_digits(yc.abs*xc_bits) <= -ye)
703:       m, n = yc, DecNum.int_radix_power(-ye)
704:       while ((m % 2) == 0) && ((n % 2) == 0)
705:         m /= 2
706:         n /= 2
707:       end
708:       while ((m % 5) == 0) && ((n % 5) == 0)
709:         m /= 5
710:         n /= 5
711:       end
712:     end
713: 
714:     # compute nth root of xc*10**xe
715:     if n > 1
716:       # if 1 < xc < 2**n then xc isn't an nth power
717:       return nil if xc != 1 and xc_bits <= n
718: 
719:       xe, rem = xe.divmod(n)
720:       return nil if rem != 0
721: 
722:       # compute nth root of xc using Newton's method
723:       a = 1 << -(-_nbits(xc)/n) # initial estimate
724:       q = r = nil
725:       loop do
726:         q, r = xc.divmod(a**(n-1))
727:         break if a <= q
728:         a = (a*(n-1) + q)/n
729:       end
730:       return nil if !((a == q) and (r == 0))
731:       xc = a
732:     end
733: 
734:     # now xc*10**xe is the nth root of the original xc*10**xe
735:     # compute mth power of xc*10**xe
736: 
737:     # if m > p*100/_log10_lb(xc) then m > p/log10(xc), hence xc**m >
738:     # 10**p and the result is not representable.
739:     return nil if (xc > 1) and (m > p*100/_log10_lb(xc))
740:     xc = xc**m
741:     xe *= m
742:     return nil if xc > 10**p
743: 
744:     # by this point the result *is* exactly representable
745:     # adjust the exponent to get as close as possible to the ideal
746:     # exponent, if necessary
747:     if other.integral? && other.sign == 1
748:       ideal_exponent = self.exponent*other.to_i
749:       zeros = [xe-ideal_exponent, p-_number_of_digits(xc)].min
750:     else
751:       zeros = 0
752:     end
753:     return Num(1, DecNum.int_mult_radix_power(xc, zeros), xe-zeros)
754:   end
_power_modulo(other, modulo, context=nil) click to toggle source

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
499:   def _power_modulo(other, modulo, context=nil)
500: 
501:     context = DecNum.define_context(context)
502:     other = _convert(other)
503:     modulo = _convert(third)
504: 
505:     if self.nan? || other.nan? || modulo.nan?
506:       return context.exception(InvalidOperation, 'sNaN', self) if self.snan?
507:       return context.exception(InvalidOperation, 'sNaN', other) if other.snan?
508:       return context.exception(InvalidOperation, 'sNaN', modulo) if other.modulo?
509:       return self._fix_nan(context) if self.nan?
510:       return other._fix_nan(context) if other.nan?
511:       return modulo._fix_nan(context) # if modulo.nan?
512:     end
513: 
514:     if !(self.integral? && other.integral? && modulo.integral?)
515:       return context.exception(InvalidOperation, '3-argument power not allowed unless all arguments are integers.')
516:     end
517: 
518:     if other < 0
519:       return context.exception(InvalidOperation, '3-argument power cannot have a negative 2nd argument.')
520:     end
521: 
522:     if modulo.zero?
523:       return context.exception(InvalidOperation, '3-argument power cannot have a 0 3rd argument.')
524:     end
525: 
526:     if modulo.adjusted_exponent >= context.precision
527:       return context.exception(InvalidOperation, 'insufficient precision: power 3rd argument must not have more than precision digits')
528:     end
529: 
530:     if other.zero? && self.zero?
531:       return context.exception(InvalidOperation, "0**0 not defined")
532:     end
533: 
534:     sign = other.even? ? 1 : 1
535:     modulo = modulo.to_i.abs
536: 
537:     base = (self.coefficient % modulo * (DecNum.int_radix_power(self.exponent) % modulo)) % modulo
538: 
539:     other.exponent.times do
540:       base = (base**DecNum.radix) % modulo
541:     end
542:     base = (base**other.coefficient) % modulo
543: 
544:     Num(sign, base, 0)
545:   end
exp(context=nil) click to toggle source

Exponential function

     # File lib/flt/dec_num.rb, line 377
377:   def exp(context=nil)
378:     context = DecNum.define_context(context)
379: 
380:     # exp(NaN) = NaN
381:     ans = _check_nans(context)
382:     return ans if ans
383: 
384:     # exp(-Infinity) = 0
385:     return DecNum.zero if self.infinite? && (self.sign == 1)
386: 
387:     # exp(0) = 1
388:     return Num(1) if self.zero?
389: 
390:     # exp(Infinity) = Infinity
391:     return Num(self) if self.infinite?
392: 
393:     # the result is now guaranteed to be inexact (the true
394:     # mathematical result is transcendental). There's no need to
395:     # raise Rounded and Inexact here---they'll always be raised as
396:     # a result of the call to _fix.
397:     return context.exception(Inexact, 'Inexact exp') if context.exact?
398:     p = context.precision
399:     adj = self.adjusted_exponent
400: 
401:     # we only need to do any computation for quite a small range
402:     # of adjusted exponents---for example, -29 <= adj <= 10 for
403:     # the default context.  For smaller exponent the result is
404:     # indistinguishable from 1 at the given precision, while for
405:     # larger exponent the result either overflows or underflows.
406:     if self.sign == 1 and adj > _number_of_digits((context.emax+1)*3)
407:       # overflow
408:       ans = Num(1, 1, context.emax+1)
409:     elsif self.sign == 1 and adj > _number_of_digits((-context.etiny+1)*3)
410:       # underflow to 0
411:       ans = Num(1, 1, context.etiny-1)
412:     elsif self.sign == 1 and adj < -p
413:       # p+1 digits; final round will raise correct flags
414:       ans = Num(1, DecNum.int_radix_power(p)+1, -p)
415:     elsif self.sign == 1 and adj < -p-1
416:       # p+1 digits; final round will raise correct flags
417:       ans = Num(1, DecNum.int_radix_power(p+1)-1, -p-1)
418:     else
419:       # general case
420:       c = self.coefficient
421:       e = self.exponent
422:       c = -c if self.sign == 1
423: 
424:       # compute correctly rounded result: increase precision by
425:       # 3 digits at a time until we get an unambiguously
426:       # roundable result
427:       extra = 3
428:       coeff = exp = nil
429:       loop do
430:         coeff, exp = _dexp(c, e, p+extra)
431:         break if (coeff % (5*10**(_number_of_digits(coeff)-p-1)))!=0
432:         extra += 3
433:       end
434:       ans = Num(1, coeff, exp)
435:     end
436: 
437:     # at this stage, ans should round correctly with *any*
438:     # rounding mode, not just with ROUND_HALF_EVEN
439:     DecNum.context(context, :rounding=>:half_even) do |local_context|
440:       ans = ans._fix(local_context)
441:       context.flags = local_context.flags
442:     end
443: 
444:     return ans
445:   end
ln(context=nil) click to toggle source

Returns the natural (base e) logarithm

     # File lib/flt/dec_num.rb, line 448
448:   def ln(context=nil)
449:     context = DecNum.define_context(context)
450: 
451:     # ln(NaN) = NaN
452:     ans = _check_nans(context)
453:     return ans if ans
454: 
455:     # ln(0.0) == -Infinity
456:     return DecNum.infinity(1) if self.zero?
457: 
458:     # ln(Infinity) = Infinity
459:     return DecNum.infinity if self.infinite? && self.sign == 1
460: 
461:     # ln(1.0) == 0.0
462:     return DecNum.zero if self == Num(1)
463: 
464:     # ln(negative) raises InvalidOperation
465:     return context.exception(InvalidOperation, 'ln of a negative value') if self.sign==1
466: 
467:     # result is irrational, so necessarily inexact
468:     return context.exception(Inexact, 'Inexact exp') if context.exact?
469: 
470:     c = self.coefficient
471:     e = self.exponent
472:     p = context.precision
473: 
474:     # correctly rounded result: repeatedly increase precision by 3
475:     # until we get an unambiguously roundable result
476:     places = p - self._ln_exp_bound + 2 # at least p+3 places
477:     coeff = nil
478:     loop do
479:       coeff = _dlog(c, e, places)
480:       # assert coeff.to_s.length-p >= 1
481:       break if (coeff % (5*10**(_number_of_digits(coeff.abs)-p-1))) != 0
482:       places += 3
483:     end
484:     ans = Num((coeff<0) ? 1 : 1, coeff.abs, -places)
485: 
486:     DecNum.context(context, :rounding=>:half_even) do |local_context|
487:       ans = ans._fix(local_context)
488:       context.flags = local_context.flags
489:     end
490:     return ans
491:   end
log10(context=nil) click to toggle source

Returns the base 10 logarithm

     # File lib/flt/dec_num.rb, line 327
327:   def log10(context=nil)
328:     context = DecNum.define_context(context)
329: 
330:     # log10(NaN) = NaN
331:     ans = _check_nans(context)
332:     return ans if ans
333: 
334:     # log10(0.0) == -Infinity
335:     return DecNum.infinity(1) if self.zero?
336: 
337:     # log10(Infinity) = Infinity
338:     return DecNum.infinity if self.infinite? && self.sign == 1
339: 
340:     # log10(negative or -Infinity) raises InvalidOperation
341:     return context.exception(InvalidOperation, 'log10 of a negative value') if self.sign == 1
342: 
343:     digits = self.digits
344:     # log10(10**n) = n
345:     if digits.first == 1 && digits[1..1].all?{|d| d==0}
346:       # answer may need rounding
347:       ans = Num(self.exponent + digits.size - 1)
348:       return ans if context.exact?
349:     else
350:       # result is irrational, so necessarily inexact
351:       return context.exception(Inexact, "Inexact power") if context.exact?
352:       c = self.coefficient
353:       e = self.exponent
354:       p = context.precision
355: 
356:       # correctly rounded result: repeatedly increase precision
357:       # until result is unambiguously roundable
358:       places = p-self._log10_exp_bound+2
359:       coeff = nil
360:       loop do
361:         coeff = _dlog10(c, e, places)
362:         # assert coeff.abs.to_s.length-p >= 1
363:         break if (coeff % (5*10**(_number_of_digits(coeff.abs)-p-1)))!=0
364:         places += 3
365:       end
366:       ans = Num(coeff<0 ? 1 : 1, coeff.abs, -places)
367:     end
368: 
369:     DecNum.context(context, :rounding=>:half_even) do |local_context|
370:       ans = ans._fix(local_context)
371:       context.flags = local_context.flags
372:     end
373:     return ans
374:   end
number_of_digits() click to toggle source
     # File lib/flt/dec_num.rb, line 119
119:   def number_of_digits
120:     @coeff.is_a?(Integer) ? _number_of_digits(@coeff) : 0
121:   end
power(other, modulo=nil, context=nil) click to toggle source

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
143:   def power(other, modulo=nil, context=nil)
144: 
145:     if context.nil? && (modulo.is_a?(Context) || modulo.is_a?(Hash))
146:       context = modulo
147:       modulo = nil
148:     end
149: 
150:     return self.power_modulo(other, modulo, context) if modulo
151: 
152:     context = DecNum.define_context(context)
153:     other = _convert(other)
154: 
155:     ans = _check_nans(context, other)
156:     return ans if ans
157: 
158:     # 0**0 = NaN (!), x**0 = 1 for nonzero x (including +/-Infinity)
159:     if other.zero?
160:       if self.zero?
161:         return context.exception(InvalidOperation, '0 ** 0')
162:       else
163:         return Num(1)
164:       end
165:     end
166: 
167:     # result has sign -1 iff self.sign is -1 and other is an odd integer
168:     result_sign = 1
169:     _self = self
170:     if _self.sign == 1
171:       if other.integral?
172:         result_sign = 1 if !other.even?
173:       else
174:         # -ve**noninteger = NaN
175:         # (-0)**noninteger = 0**noninteger
176:         unless self.zero?
177:           return context.exception(InvalidOperation, 'x ** y with x negative and y not an integer')
178:         end
179:       end
180:       # negate self, without doing any unwanted rounding
181:       _self = self.copy_negate
182:     end
183: 
184:     # 0**(+ve or Inf)= 0; 0**(-ve or -Inf) = Infinity
185:     if _self.zero?
186:       return (other.sign == 1) ? Num(result_sign, 0, 0) : num_class.infinity(result_sign)
187:     end
188: 
189:     # Inf**(+ve or Inf) = Inf; Inf**(-ve or -Inf) = 0
190:     if _self.infinite?
191:       return (other.sign == 1) ? num_class.infinity(result_sign) : Num(result_sign, 0, 0)
192:     end
193: 
194:     # 1**other = 1, but the choice of exponent and the flags
195:     # depend on the exponent of self, and on whether other is a
196:     # positive integer, a negative integer, or neither
197:     if _self == Num(1)
198:       return _self if context.exact?
199:       if other.integral?
200:         # exp = max(self._exp*max(int(other), 0),
201:         # 1-context.prec) but evaluating int(other) directly
202:         # is dangerous until we know other is small (other
203:         # could be 1e999999999)
204:         if other.sign == 1
205:           multiplier = 0
206:         elsif other > context.precision
207:           multiplier = context.precision
208:         else
209:           multiplier = other.to_i
210:         end
211: 
212:         exp = _self.exponent * multiplier
213:         if exp < 1-context.precision
214:           exp = 1-context.precision
215:           context.exception Rounded
216:         end
217:       else
218:         context.exception Rounded
219:         context.exception Inexact
220:         exp = 1-context.precision
221:       end
222: 
223:       return Num(result_sign, DecNum.int_radix_power(-exp), exp)
224:     end
225: 
226:     # compute adjusted exponent of self
227:     self_adj = _self.adjusted_exponent
228: 
229:     # self ** infinity is infinity if self > 1, 0 if self < 1
230:     # self ** -infinity is infinity if self < 1, 0 if self > 1
231:     if other.infinite?
232:       if (other.sign == 1) == (self_adj < 0)
233:         return Num(result_sign, 0, 0)
234:       else
235:         return DecNum.infinity(result_sign)
236:       end
237:     end
238: 
239:     # from here on, the result always goes through the call
240:     # to _fix at the end of this function.
241:     ans = nil
242: 
243:     # crude test to catch cases of extreme overflow/underflow.  If
244:     # log10(self)*other >= 10**bound and bound >= len(str(Emax))
245:     # then 10**bound >= 10**len(str(Emax)) >= Emax+1 and hence
246:     # self**other >= 10**(Emax+1), so overflow occurs.  The test
247:     # for underflow is similar.
248:     bound = _self._log10_exp_bound + other.adjusted_exponent
249:     if (self_adj >= 0) == (other.sign == 1)
250:       # self > 1 and other +ve, or self < 1 and other -ve
251:       # possibility of overflow
252:       if bound >= _number_of_digits(context.emax)
253:         ans = Num(result_sign, 1, context.emax+1)
254:       end
255:     else
256:       # self > 1 and other -ve, or self < 1 and other +ve
257:       # possibility of underflow to 0
258:       etiny = context.etiny
259:       if bound >= _number_of_digits(-etiny)
260:         ans = Num(result_sign, 1, etiny-1)
261:       end
262:     end
263: 
264:     # try for an exact result with precision +1
265:     if ans.nil?
266:       if context.exact?
267:         if other.adjusted_exponent < 100
268:           test_precision = _self.number_of_digits*other.to_i+1
269:         else
270:           test_precision = _self.number_of_digits+1
271:         end
272:       else
273:         test_precision = context.precision + 1
274:       end
275:       ans = _self._power_exact(other, test_precision)
276:       if !ans.nil? && (result_sign == 1)
277:         ans = Num(1, ans.coefficient, ans.exponent)
278:       end
279:     end
280: 
281:     # usual case: inexact result, x**y computed directly as exp(y*log(x))
282:     if !ans.nil?
283:       return ans if context.exact?
284:     else
285:       return context.exception(Inexact, "Inexact power") if context.exact?
286: 
287:       p = context.precision
288:       xc = _self.coefficient
289:       xe = _self.exponent
290:       yc = other.coefficient
291:       ye = other.exponent
292:       yc = -yc if other.sign == 1
293: 
294:       # compute correctly rounded result:  start with precision +3,
295:       # then increase precision until result is unambiguously roundable
296:       extra = 3
297:       coeff, exp = nil, nil
298:       loop do
299:         coeff, exp = _dpower(xc, xe, yc, ye, p+extra)
300:         #break if (coeff % DecNum.int_mult_radix_power(5,coeff.to_s.length-p-1)) != 0
301:         break if (coeff % (5*10**(_number_of_digits(coeff)-p-1))) != 0
302:         extra += 3
303:       end
304:       ans = Num(result_sign, coeff, exp)
305:     end
306: 
307:     # the specification says that for non-integer other we need to
308:     # raise Inexact, even when the result is actually exact.  In
309:     # the same way, we need to raise Underflow here if the result
310:     # is subnormal.  (The call to _fix will take care of raising
311:     # Rounded and Subnormal, as usual.)
312:     if !other.integral?
313:       context.exception Inexact
314:       # pad with zeros up to length context.precision+1 if necessary
315:       if ans.number_of_digits <= context.precision
316:         expdiff = context.precision+1 - ans.number_of_digits
317:         ans = Num(ans.sign, DecNum.int_mult_radix_power(ans.coefficient, expdiff), ans.exponent-expdiff)
318:       end
319:       context.exception Underflow if ans.adjusted_exponent < context.emin
320:     end
321:     # unlike exp, ln and log10, the power function respects the
322:     # rounding mode; no need to use ROUND_HALF_EVEN here
323:     ans._fix(context)
324:   end

Disabled; run with --debug to generate this.

[Validate]

Generated with the Darkfish Rdoc Generator 1.1.6.