Module: Flt::DecNum::AuxiliarFunctions

Included in:
Flt::DecNum
Defined in:
lib/flt/dec_num.rb

Overview

:nodoc:

Constant Summary collapse

LOG10_LB_CORRECTION =

(1..9).map_hash{|i| 100 - (100*Math.log10(i)).floor}

{ # (1..9).map_hash{|i| 100 - (100*Math.log10(i)).floor}
'1'=> 100, '2'=> 70, '3'=> 53, '4'=> 40, '5'=> 31,
'6'=> 23, '7'=> 16, '8'=> 10, '9'=> 5}
NUMBER_OF_DIGITS_MAX_VALID_LOG = 
10**(Float::DIG-1)

Class Attribute Summary collapse

Class Method Summary collapse

Class Attribute Details

.log10_digits ⇒ Object

Returns the value of attribute log10_digits.



1122
1123
1124
 
# File 'lib/flt/dec_num.rb', line 1122

def log10_digits
  @log10_digits
end

Class Method Details

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



874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
 
# File 'lib/flt/dec_num.rb', line 874

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
    # TODO: without the .abs tests fail because c is negative: c should not be negative!!
    extra = [0, e + _number_of_digits(c.abs) - 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.



902
903
904
905
 
# File 'lib/flt/dec_num.rb', line 902

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.



1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
 
# File 'lib/flt/dec_num.rb', line 1033

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 = _number_of_digits(c)
    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 = _number_of_digits(f.abs) - 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.



988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
 
# File 'lib/flt/dec_num.rb', line 988

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 = _number_of_digits(c)
  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.



832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
 
# File 'lib/flt/dec_num.rb', line 832

def _dpower(xc, xe, yc, ye, p)
  # Find b such that 10**(b-1) <= |y| <= 10**b
  b = _number_of_digits(yc.abs) + 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 (_number_of_digits(xc) + 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).



1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
 
# File 'lib/flt/dec_num.rb', line 1083

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*_number_of_digits(m)/(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.



941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
 
# File 'lib/flt/dec_num.rb', line 941

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*_number_of_digits(m)/(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) 


1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
 
# File 'lib/flt/dec_num.rb', line 1126

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) ⇒ Object

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

Raises:

  • (ArgumentError) 


1021
1022
1023
1024
1025
 
# File 'lib/flt/dec_num.rb', line 1021

def _log10_lb(c)
    raise ArgumentError, "The argument to _log10_lb should be nonnegative." if c <= 0
    str_c = c.to_s
    return 100*str_c.length - LOG10_LB_CORRECTION[str_c[0,1]]
end

._number_of_digits(i) ⇒ Object

number of bits in a nonnegative integer

Raises:

  • (TypeError) 


1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
 
# File 'lib/flt/dec_num.rb', line 1189

def _number_of_digits(i)
  raise  TypeError, "The argument to _number_of_digits should be nonnegative." if i < 0
  if i.is_a?(Fixnum) || (i > NUMBER_OF_DIGITS_MAX_VALID_LOG)
    # for short integers this is faster
    # note that here we return 1 for 0
    i.to_s.length
  else
    (::Math.log10(i)+1).floor
  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.



926
927
928
929
930
 
# File 'lib/flt/dec_num.rb', line 926

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.



911
912
913
914
915
916
917
918
919
920
921
922
 
# File 'lib/flt/dec_num.rb', line 911

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.



1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
 
# File 'lib/flt/dec_num.rb', line 1163

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 + _number_of_digits(c) - 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 < 1s20;  error after division < 0.62
  return _div_nearest(_iexp(rem, 10**p), 1000), quot - p + 3
end