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.



4086
4087
4088
 
# File 'lib/decimal/decimal.rb', line 4086

def log10_digits
  @log10_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



3726
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
 
# File 'lib/decimal/decimal.rb', line 3726

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.



3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863
 
# File 'lib/decimal/decimal.rb', line 3840

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.



3867
3868
3869
3870
 
# File 'lib/decimal/decimal.rb', line 3867

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.



3997
3998
3999
4000
4001
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
4030
4031
4032
4033
4034
4035
4036
4037
4038
4039
4040
4041
 
# File 'lib/decimal/decimal.rb', line 3997

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.



3953
3954
3955
3956
3957
3958
3959
3960
3961
3962
3963
3964
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979
3980
3981
3982
3983
 
# File 'lib/decimal/decimal.rb', line 3953

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.



3798
3799
3800
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
3828
 
# File 'lib/decimal/decimal.rb', line 3798

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



4047
4048
4049
4050
4051
4052
4053
4054
4055
4056
4057
4058
4059
4060
4061
4062
4063
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080
4081
 
# File 'lib/decimal/decimal.rb', line 4047

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.



3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923
3924
3925
3926
3927
3928
3929
3930
3931
3932
3933
3934
3935
3936
3937
3938
3939
3940
3941
3942
3943
3944
3945
3946
3947
3948
 
# File 'lib/decimal/decimal.rb', line 3906

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) 


4090
4091
4092
4093
4094
4095
4096
4097
4098
4099
4100
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112
4113
 
# File 'lib/decimal/decimal.rb', line 4090

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) 


3986
3987
3988
3989
3990
3991
3992
 
# File 'lib/decimal/decimal.rb', line 3986

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) 


3776
3777
3778
3779
3780
3781
3782
3783
3784
 
# File 'lib/decimal/decimal.rb', line 3776

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.



3748
3749
3750
3751
3752
3753
3754
3755
3756
3757
3758
3759
3760
3761
3762
3763
3764
3765
3766
 
# File 'lib/decimal/decimal.rb', line 3748

def _normalize(op1, op2, prec=0)
  if op1.exponent < op2.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.exponent + [-1, tmp_len - prec - 2].min
  if (other_len+other.exponent-1 < exp) && prec>0
    other = Decimal.new([other.sign, 1, exp])
  end
  tmp = Decimal.new(tmp.sign,
                    Decimal.int_mult_radix_power(tmp.coefficient, tmp.exponent-other.exponent),
                    other.exponent)
  return swap ? [other, tmp] : [tmp, other]
end

._parser(txt) ⇒ Object

Parse numeric text literals (internal use)



3739
3740
3741
3742
3743
3744
3745
 
# File 'lib/decimal/decimal.rb', line 3739

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.



3891
3892
3893
3894
3895
 
# File 'lib/decimal/decimal.rb', line 3891

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.



3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886
3887
 
# File 'lib/decimal/decimal.rb', line 3876

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.



4127
4128
4129
4130
4131
4132
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
4150
 
# File 'lib/decimal/decimal.rb', line 4127

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