Module: Decimal::AuxiliarFunctions

Included in:
Decimal
Defined in:
lib/decimal/decimal.rb

Overview

:nodoc:

Class Attribute Summary collapse

Class Method Summary collapse

Class Attribute Details

.log10_digits ⇒ Object

Returns the value of attribute log10_digits.



3859
3860
3861
 
# File 'lib/decimal/decimal.rb', line 3859

def log10_digits
  @log10_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



3499
3500
3501
3502
3503
3504
3505
3506
3507
3508
3509
 
# File 'lib/decimal/decimal.rb', line 3499

def _convert(x, error=true)
  case x
  when Decimal
    x
  when *Decimal.context.coercible_types
    Decimal.new(x)
  else
    raise TypeError, "Unable to convert #{x.class} to Decimal" if error
    nil
  end
end

._dexp(c, e, p) ⇒ Object

Compute an approximation to exp(c*10**e), with p decimal places of precision. Returns integers d, f such that:

10**(p-1) <= d <= 10**p, and
(d-1)*10**f < exp(c*10**e) < (d+1)*10**f

In other words, d*10**f is an approximation to exp(c*10**e) with p digits of precision, and with an error in d of at most 1. This is almost, but not quite, the same as the error being < 1ulp: when d

10**(p-1) the error could be up to 10 ulp.



3613
3614
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627
3628
3629
3630
3631
3632
3633
3634
3635
3636
 
# File 'lib/decimal/decimal.rb', line 3613

def _dexp(c, e, p)
    # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision
    p += 2

    # compute log(10) with extra precision = adjusted exponent of c*10**e
    extra = [0, e + c.to_s.length - 1].max
    q = p + extra

    # compute quotient c*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q),
    # rounding down
    shift = e+q
    if shift >= 0
        cshift = c*10**shift
    else
        cshift = c/10**-shift
    end
    quot, rem = cshift.divmod(_log10_digits(q))

    # reduce remainder back to original precision
    rem = _div_nearest(rem, 10**extra)

    # error in result of _iexp < 120;  error after division < 0.62
    return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
end

._div_nearest(a, b) ⇒ Object

Closest integer to a/b, a and b positive integers; rounds to even in the case of a tie.



3640
3641
3642
3643
 
# File 'lib/decimal/decimal.rb', line 3640

def _div_nearest(a, b)
  q, r = a.divmod(b)
  q + (((2*r + (q&1)) > b) ? 1 : 0)
end

._dlog(c, e, p) ⇒ Object

Given integers c, e and p with c > 0, compute an integer approximation to 10**p * log(c*10**e), with an absolute error of at most 1. Assumes that c*10**e is not exactly 1.



3770
3771
3772
3773
3774
3775
3776
3777
3778
3779
3780
3781
3782
3783
3784
3785
3786
3787
3788
3789
3790
3791
3792
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
 
# File 'lib/decimal/decimal.rb', line 3770

def _dlog(c, e, p)

    # Increase precision by 2. The precision increase is compensated
    # for at the end with a division by 100.
    p += 2

    # rewrite c*10**e as d*10**f with either f >= 0 and 1 <= d <= 10,
    # or f <= 0 and 0.1 <= d <= 1.  Then we can compute 10**p * log(c*10**e)
    # as 10**p * log(d) + 10**p*f * log(10).
    l = c.to_s.length
    f = e+l - ((e+l >= 1) ? 1 : 0)

    # compute approximation to 10**p*log(d), with error < 27
    if p > 0
        k = e+p-f
        if k >= 0
            c *= 10**k
        else
            c = _div_nearest(c, 10**-k)  # error of <= 0.5 in c
        end

        # _ilog magnifies existing error in c by a factor of at most 10
        log_d = _ilog(c, 10**p) # error < 5 + 22 = 27
    else
        # p <= 0: just approximate the whole thing by 0; error < 2.31
        log_d = 0
    end

    # compute approximation to f*10**p*log(10), with error < 11.
    if f
        extra = f.abs.to_s.length - 1
        if p + extra >= 0
            # error in f * _log10_digits(p+extra) < |f| * 1 = |f|
            # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
            f_log_ten = _div_nearest(f*_log10_digits(p+extra), 10**extra)
        else
            f_log_ten = 0
        end
    else
        f_log_ten = 0
    end

    # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1
    return _div_nearest(f_log_ten + log_d, 100)
end

._dlog10(c, e, p) ⇒ Object

Given integers c, e and p with c > 0, p >= 0, compute an integer approximation to 10**p * log10(c*10**e), with an absolute error of at most 1. Assumes that c*10**e is not exactly 1.



3726
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
3737
3738
3739
3740
3741
3742
3743
3744
3745
3746
3747
3748
3749
3750
3751
3752
3753
3754
3755
3756
 
# File 'lib/decimal/decimal.rb', line 3726

def _dlog10(c, e, p)
   # increase precision by 2; compensate for this by dividing
  # final result by 100
  p += 2

  # write c*10**e as d*10**f with either:
  #   f >= 0 and 1 <= d <= 10, or
  #   f <= 0 and 0.1 <= d <= 1.
  # Thus for c*10**e close to 1, f = 0
  l = c.to_s.length
  f = e+l - ((e+l >= 1) ? 1 : 0)

  if p > 0
    m = 10**p
    k = e+p-f
    if k >= 0
      c *= 10**k
    else
      c = _div_nearest(c, 10**-k)
    end
    log_d = _ilog(c, m) # error < 5 + 22 = 27
    log_10 = _log10_digits(p) # error < 1
    log_d = _div_nearest(log_d*m, log_10)
    log_tenpower = f*m # exact
  else
    log_d = 0  # error < 2.31
    log_tenpower = _div_nearest(f, 10**-p) # error < 0.5
  end

  return _div_nearest(log_tenpower+log_d, 100)
end

._dpower(xc, xe, yc, ye, p) ⇒ Object

Given integers xc, xe, yc and ye representing Decimals x = xc*10**xe and y = yc*10**ye, compute x**y. Returns a pair of integers (c, e) such that:

10**(p-1) <= c <= 10**p, and
(c-1)*10**e < x**y < (c+1)*10**e

in other words, c*10**e is an approximation to x**y with p digits of precision, and with an error in c of at most 1. (This is almost, but not quite, the same as the error being < 1ulp: when c

10**(p-1) we can only guarantee error < 10ulp.)

We assume that: x is positive and not equal to 1, and y is nonzero.



3571
3572
3573
3574
3575
3576
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
3595
3596
3597
3598
3599
3600
3601
 
# File 'lib/decimal/decimal.rb', line 3571

def _dpower(xc, xe, yc, ye, p)
  # Find b such that 10**(b-1) <= |y| <= 10**b
  b = yc.abs.to_s.length + ye

  # log(x) = lxc*10**(-p-b-1), to p+b+1 places after the decimal point
  lxc = _dlog(xc, xe, p+b+1)

  # compute product y*log(x) = yc*lxc*10**(-p-b-1+ye) = pc*10**(-p-1)
  shift = ye-b
  if shift >= 0
      pc = lxc*yc*10**shift
  else
      pc = _div_nearest(lxc*yc, 10**-shift)
  end

  if pc == 0
      # we prefer a result that isn't exactly 1; this makes it
      # easier to compute a correctly rounded result in __pow__
      if (xc.to_s.length + xe >= 1) == (yc > 0) # if x**y > 1:
          coeff, exp = 10**(p-1)+1, 1-p
      else
          coeff, exp = 10**p-1, -p
      end
  else
      coeff, exp = _dexp(pc, -(p+1), p+1)
      coeff = _div_nearest(coeff, 10)
      exp += 1
  end

  return coeff, exp
end

._iexp(x, m, l = 8) ⇒ Object

Given integers x and M, M > 0, such that x/M is small in absolute value, compute an integer approximation to M*exp(x/M). For 0 <= x/M <= 2.4, the absolute error in the result is bounded by 60 (and is usually much smaller).



3820
3821
3822
3823
3824
3825
3826
3827
3828
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
 
# File 'lib/decimal/decimal.rb', line 3820

def _iexp(x, m, l=8)

    # Algorithm: to compute exp(z) for a real number z, first divide z
    # by a suitable power R of 2 so that |z/2**R| < 2**-L.  Then
    # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor
    # series
    #
    #     expm1(x) = x + x**2/2! + x**3/3! + ...
    #
    # Now use the identity
    #
    #     expm1(2x) = expm1(x)*(expm1(x)+2)
    #
    # R times to compute the sequence expm1(z/2**R),
    # expm1(z/2**(R-1)), ... , exp(z/2), exp(z).

    # Find R such that x/2**R/M <= 2**-L
    r = _nbits((x<<l)/m)

    # Taylor series.  (2**L)**T > M
    t = -(-10*m.to_s.length/(3*l)).to_i
    y = _div_nearest(x, t)
    mshift = m<<r
    (1...t).to_a.reverse.each do |i|
        y = _div_nearest(x*(mshift + y), mshift * i)
    end

    # Expansion
    (0...r).to_a.reverse.each do |k|
        mshift = m<<(k+2)
        y = _div_nearest(y*(y+mshift), mshift)
    end

    return m+y
end

._ilog(x, m, l = 8) ⇒ Object

Integer approximation to M*log(x/M), with absolute error boundable in terms only of x/M.

Given positive integers x and M, return an integer approximation to M * log(x/M). For L = 8 and 0.1 <= x/M <= 10 the difference between the approximation and the exact result is at most 22. For L = 8 and 1.0 <= x/M <= 10.0 the difference is at most 15. In both cases these are upper bounds on the error; it will usually be much smaller.



3679
3680
3681
3682
3683
3684
3685
3686
3687
3688
3689
3690
3691
3692
3693
3694
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
 
# File 'lib/decimal/decimal.rb', line 3679

def _ilog(x, m, l = 8)
  # The basic algorithm is the following: let log1p be the function
  # log1p(x) = log(1+x).  Then log(x/M) = log1p((x-M)/M).  We use
  # the reduction
  #
  #    log1p(y) = 2*log1p(y/(1+sqrt(1+y)))
  #
  # repeatedly until the argument to log1p is small (< 2**-L in
  # absolute value).  For small y we can use the Taylor series
  # expansion
  #
  #    log1p(y) ~ y - y**2/2 + y**3/3 - ... - (-y)**T/T
  #
  # truncating at T such that y**T is small enough.  The whole
  # computation is carried out in a form of fixed-point arithmetic,
  # with a real number z being represented by an integer
  # approximation to z*M.  To avoid loss of precision, the y below
  # is actually an integer approximation to 2**R*y*M, where R is the
  # number of reductions performed so far.

  y = x-m
  # argument reduction; R = number of reductions performed
  r = 0
  # while (r <= l && y.abs << l-r >= m ||
  #        r > l and y.abs>> r-l >= m)
  while (((r <= l) && ((y.abs << (l-r)) >= m)) ||
         ((r > l) && ((y.abs>>(r-l)) >= m)))
      y = _div_nearest((m*y) << 1,
                       m + _sqrt_nearest(m*(m+_rshift_nearest(y, r)), m))
      r += 1
  end

  # Taylor series with T terms
  t = -(-10*m.to_s.length/(3*l)).to_i
  yshift = _rshift_nearest(y, r)
  w = _div_nearest(m, t)
  # (1...t).reverse_each do |k| # Ruby 1.9
  (1...t).to_a.reverse.each do |k|
     w = _div_nearest(m, k) - _div_nearest(yshift*w, m)
  end

  return _div_nearest(w*y, m)
end

._log10_digits(p) ⇒ Object

Given an integer p >= 0, return floor(10**p)*log(10).

Raises:

  • (ArgumentError) 


3863
3864
3865
3866
3867
3868
3869
3870
3871
3872
3873
3874
3875
3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886
 
# File 'lib/decimal/decimal.rb', line 3863

def _log10_digits(p)
  # digits are stored as a string, for quick conversion to
  # integer in the case that we've already computed enough
  # digits; the stored digits should always be correct
  # (truncated, not rounded to nearest).
  raise ArgumentError, "p should be nonnegative" if p<0
  if p >= AuxiliarFunctions.log10_digits.length
      digits = nil
      # compute p+3, p+6, p+9, ... digits; continue until at
      # least one of the extra digits is nonzero
      extra = 3
      loop do
        # compute p+extra digits, correct to within 1ulp
        m = 10**(p+extra+2)
        digits = _div_nearest(_ilog(10*m, m), 100).to_s
        break if digits[-extra..-1] != '0'*extra
        extra += 3
      end
      # keep all reliable digits so far; remove trailing zeros
      # and next nonzero digit
      AuxiliarFunctions.log10_digits = digits.sub(/0*$/,'')[0...-1]
  end
  return (AuxiliarFunctions.log10_digits[0...p+1]).to_i
end

._log10_lb(c, correction = { '1'=> 100, '2'=> 70, '3'=> 53, '4'=> 40, '5'=> 31, '6'=> 23, '7'=> 16, '8'=> 10, '9'=> 5}) ⇒ Object

Compute a lower bound for 100*log10© for a positive integer c.

Raises:

  • (ArgumentError) 


3759
3760
3761
3762
3763
3764
3765
 
# File 'lib/decimal/decimal.rb', line 3759

def _log10_lb(c, correction = {
        '1'=> 100, '2'=> 70, '3'=> 53, '4'=> 40, '5'=> 31,
        '6'=> 23, '7'=> 16, '8'=> 10, '9'=> 5})
    raise ArgumentError, "The argument to _log10_lb should be nonnegative." if c <= 0
    str_c = c.to_s
    return 100*str_c.length - correction[str_c[0,1]]
end

._nbits(n, correction = { #:nodoc: '0'=> 4, '1'=> 3, '2'=> 2, '3'=> 2, '4'=> 1, '5'=> 1, '6'=> 1, '7'=> 1, '8'=> 0, '9'=> 0, 'a'=> 0, 'b'=> 0, 'c'=> 0, 'd'=> 0, 'e'=> 0, 'f'=> 0}) ⇒ Object

Number of bits in binary representation of the positive integer n, or 0 if n == 0. – This function from Tim Peters was taken from here: mail.python.org/pipermail/python-list/1999-July/007758.html The correction being in the function definition is for speed, and the whole function is not resolved with math.log because of avoiding the use of floats. ++

Raises:

  • (TypeError) 


3549
3550
3551
3552
3553
3554
3555
3556
3557
 
# File 'lib/decimal/decimal.rb', line 3549

def _nbits(n, correction = {                      #:nodoc:
        '0'=> 4, '1'=> 3, '2'=> 2, '3'=> 2,
        '4'=> 1, '5'=> 1, '6'=> 1, '7'=> 1,
        '8'=> 0, '9'=> 0, 'a'=> 0, 'b'=> 0,
        'c'=> 0, 'd'=> 0, 'e'=> 0, 'f'=> 0})
  raise  TypeError, "The argument to _nbits should be nonnegative." if n < 0
  hex_n = "%x" % n
  4*hex_n.length - correction[hex_n[0,1]]
end

._normalize(op1, op2, prec = 0) ⇒ Object

Normalizes op1, op2 to have the same exp and length of coefficient. Used for addition.



3521
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533
3534
3535
3536
3537
3538
3539
 
# File 'lib/decimal/decimal.rb', line 3521

def _normalize(op1, op2, prec=0)
  if op1.integral_exponent < op2.integral_exponent
    swap = true
    tmp,other = op2,op1
  else
    swap = false
    tmp,other = op1,op2
  end
  tmp_len = tmp.number_of_digits
  other_len = other.number_of_digits
  exp = tmp.integral_exponent + [-1, tmp_len - prec - 2].min
  if (other_len+other.integral_exponent-1 < exp) && prec>0
    other = Decimal.new([other.sign, 1, exp])
  end
  tmp = Decimal.new(tmp.sign,
                    Decimal.int_mult_radix_power(tmp.integral_significand, tmp.integral_exponent-other.integral_exponent),
                    other.integral_exponent)
  return swap ? [other, tmp] : [tmp, other]
end

._parser(txt) ⇒ Object

Parse numeric text literals (internal use)



3512
3513
3514
3515
3516
3517
3518
 
# File 'lib/decimal/decimal.rb', line 3512

def _parser(txt)
  md = /^\s*([-+])?(?:(?:(\d+)(?:\.(\d*))?|\.(\d+))(?:[eE]([-+]?\d+))?|Inf(?:inity)?|(s)?NaN(\d*))\s*$/i.match(txt)
  if md
    OpenStruct.new :sign=>md[1], :int=>md[2], :frac=>md[3], :onlyfrac=>md[4], :exp=>md[5],
                   :signal=>md[6], :diag=>md[7]
  end
end

._rshift_nearest(x, shift) ⇒ Object

Given an integer x and a nonnegative integer shift, return closest integer to x / 2**shift; use round-to-even in case of a tie.



3664
3665
3666
3667
3668
 
# File 'lib/decimal/decimal.rb', line 3664

def _rshift_nearest(x, shift)
    b, q = (1 << shift), (x >> shift)
    return q + (((2*(x & (b-1)) + (q&1)) > b) ? 1 : 0)
    #return q + (2*(x & (b-1)) + (((q&1) > b) ? 1 : 0))
end

._sqrt_nearest(n, a) ⇒ Object

Closest integer to the square root of the positive integer n. a is an initial approximation to the square root. Any positive integer will do for a, but the closer a is to the square root of n the faster convergence will be.



3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
3660
 
# File 'lib/decimal/decimal.rb', line 3649

def _sqrt_nearest(n, a)

    if n <= 0 or a <= 0
        raise ArgumentError, "Both arguments to _sqrt_nearest should be positive."
    end

    b=0
    while a != b
        b, a = a, a--n/a>>1 # ??
    end
    return a
end

.dexp(c, e, p) ⇒ Object

Compute an approximation to exp(c*10**e), with p decimal places of precision.

Returns integers d, f such that:

10**(p-1) <= d <= 10**p, and
(d-1)*10**f < exp(c*10**e) < (d+1)*10**f

In other words, d*10**f is an approximation to exp(c*10**e) with p digits of precision, and with an error in d of at most 1. This is almost, but not quite, the same as the error being < 1ulp: when d

10**(p-1) the error could be up to 10 ulp.



3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923
 
# File 'lib/decimal/decimal.rb', line 3900

def dexp(c, e, p)
  # we'll call iexp with M = 10**(p+2), giving p+3 digits of precision
  p += 2

  # compute log(10) with extra precision = adjusted exponent of c*10**e
  extra = [0, e + c.to_s.length - 1].max
  q = p + extra

  # compute quotient c*10**e/(log(10)) = c*10**(e+q)/(log(10)*10**q),
  # rounding down
  shift = e+q
  if shift >= 0
      cshift = c*10**shift
  else
      cshift = c/10**-shift
  end
  quot, rem = cshift.divmod(_log10_digits(q))

  # reduce remainder back to original precision
  rem = _div_nearest(rem, 10**extra)

  # error in result of _iexp < 120;  error after division < 0.62
  return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
end