Module: Flt::Num::AuxiliarFunctions

Included in:
Flt::Num, Flt::Num
Defined in:
lib/flt/num.rb

Constant Summary collapse

EXP_INC = 
4
LOG_PREC_INC = 
4
LOG_RADIX_INC = 
2
LOG_RADIX_EXTRA = 
3
LOG2_MULT =

TODO: K=100? K=64? …

100
LOG2_LB_CORRECTION =

(1..15).map{|i| (LOG2_MULT*Math.log(16.0/i)/Math.log(2)).ceil}

[ # (1..15).map{|i| (LOG2_MULT*Math.log(16.0/i)/Math.log(2)).ceil}
  400, 300, 242, 200, 168, 142, 120, 100, 84, 68, 55, 42, 30, 20, 10
  # for LOG2_MULT=64: 256, 192, 155, 128, 108, 91, 77, 64, 54, 44, 35, 27, 20, 13, 6
]
LOG10_MULT = 
100
LOG10_LB_CORRECTION =

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

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

Class Attribute Summary collapse

Class Method Summary collapse

Class Attribute Details

.log_radix_digits ⇒ Object (readonly)

Returns the value of attribute log_radix_digits.



4479
4480
4481
 
# File 'lib/flt/num.rb', line 4479

def log_radix_digits
  @log_radix_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



4166
4167
4168
4169
4170
4171
4172
4173
4174
4175
4176
 
# File 'lib/flt/num.rb', line 4166

def _convert(x, error=true)
  case x
  when num_class
    x
  when *num_class.context.coercible_types
    num_class.new(x)
  else
    raise TypeError, "Unable to convert #{x.class} to #{num_class}" if error
    nil
  end
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.



4468
4469
4470
4471
 
# File 'lib/flt/num.rb', line 4468

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

._exp(c, e, p) ⇒ Object

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

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

In other words, d*radix**f is an approximation to exp(c*radix**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

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



4272
4273
4274
4275
4276
4277
4278
4279
4280
4281
4282
4283
4284
4285
4286
4287
4288
4289
4290
4291
4292
4293
4294
4295
4296
4297
 
# File 'lib/flt/num.rb', line 4272

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

    # compute log(radix) with extra precision = adjusted exponent of c*radix**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*radix**e/(log(radix)) = c*radix**(e+q)/(log(radix)*radix**q),
    # rounding down
    shift = e+q
    if shift >= 0
        cshift = c*num_class.int_radix_power(shift)
    else
        cshift = c/num_class.int_radix_power(-shift)
    end
    quot, rem = cshift.divmod(_log_radix_digits(q))

    # reduce remainder back to original precision
    rem = _div_nearest(rem, num_class.int_radix_power(extra))

    # for radix=10: error in result of _iexp < 120;  error after division < 0.62
    r = _div_nearest(_iexp(rem, num_class.int_radix_power(p)), num_class.int_radix_power(EXP_INC+1)), quot - p + (EXP_INC+1)
    return r
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 redix=10, and 0 <= x/M <= 2.4, the absolute error in the result is bounded by 60 (and is usually much smaller).



4353
4354
4355
4356
4357
4358
4359
4360
4361
4362
4363
4364
4365
4366
4367
4368
4369
4370
4371
4372
4373
4374
4375
4376
4377
4378
4379
4380
4381
4382
4383
4384
4385
4386
4387
 
# File 'lib/flt/num.rb', line 4353

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 = -(-num_class.radix*_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 radix=10, 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.



4398
4399
4400
4401
4402
4403
4404
4405
4406
4407
4408
4409
4410
4411
4412
4413
4414
4415
4416
4417
4418
4419
4420
4421
4422
4423
4424
4425
4426
4427
4428
4429
4430
4431
4432
4433
4434
4435
4436
4437
4438
4439
 
# File 'lib/flt/num.rb', line 4398

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

._log(c, e, p) ⇒ Object

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



4303
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
4317
4318
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
4329
4330
4331
4332
4333
4334
4335
4336
4337
4338
4339
4340
4341
4342
4343
4344
4345
4346
4347
 
# File 'lib/flt/num.rb', line 4303

def _log(c, e, p)
    # Increase precision by 2. The precision increase is compensated
    # for at the end with a division
    p += LOG_PREC_INC

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

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

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

    # compute approximation to f*radix**p*log(radix), with error < 11 for radix=10.
    if f
        extra = _number_of_digits(f.abs) - 1
        if p + extra >= 0
            # for radix=10:
            # error in f * _log10_digits(p+extra) < |f| * 1 = |f|
            # after division, error < |f|/10**extra + 0.5 < 10 + 0.5 < 11
            f_log_r = _div_nearest(f*_log_radix_digits(p+extra), num_class.int_radix_power(extra))
        else
            f_log_r = 0
        end
    else
        f_log_r = 0
    end

    # error in sum < 11+27 = 38; error after division < 0.38 + 0.5 < 1 for radix=10
    return _div_nearest(f_log_r + log_d, num_class.int_radix_power(LOG_PREC_INC)) # extra radix factor for base 2 ???
end

._log_radix_digits(p) ⇒ Object

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

Raises:

  • (ArgumentError) 


4485
4486
4487
4488
4489
4490
4491
4492
4493
4494
4495
4496
4497
4498
4499
4500
4501
4502
4503
4504
4505
4506
4507
4508
4509
4510
4511
4512
 
# File 'lib/flt/num.rb', line 4485

def _log_radix_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 bge correct
  # (truncated, not rounded to nearest).
  raise ArgumentError, "p should be nonnegative" if p<0
  stored_digits = (AuxiliarFunctions.log_radix_digits[num_class.radix] || "")
  if p >= stored_digits.length
      digits = nil
      # compute p+3, p+6, p+9, ... digits; continue until at
      # least one of the extra digits is nonzero
      extra = LOG_RADIX_EXTRA
      loop do
        # compute p+extra digits, correct to within 1ulp
        m = num_class.int_radix_power(p+extra+LOG_RADIX_INC)
        digits = _div_nearest(_ilog(num_class.radix*m, m), num_class.int_radix_power(LOG_RADIX_INC)).to_s(num_class.radix)
        break if digits[-extra..-1] != '0'*extra
        extra += LOG_RADIX_EXTRA
      end
      # if the radix < e (i.e. only for radix==2), we must prefix with a 0 because log(radix)<1
      # BUT THIS REDUCES PRECISION BY ONE? : may be avoid prefix and adjust scaling in the caller
      prefix = num_class.radix==2 ? '0' : ''
      # keep all reliable digits so far; remove trailing zeros
      # and next nonzero digit
      AuxiliarFunctions.log_radix_digits[num_class.radix] = prefix + digits.sub(/0*$/,'')[0...-1]
  end
  return (AuxiliarFunctions.log_radix_digits[num_class.radix][0..p]).to_i(num_class.radix)
end

._log_radix_lb(c) ⇒ Object 



4549
4550
4551
4552
4553
4554
4555
4556
4557
4558
 
# File 'lib/flt/num.rb', line 4549

def _log_radix_lb(c)
  case num_class.radix
  when 10
    log10_lb(c)
  when 2
    log2_lb(c)
  else
    raise ArgumentError, "_log_radix_lb not implemented for base #{num_class.radix}"
  end
end

._log_radix_mult ⇒ Object 



4538
4539
4540
4541
4542
4543
4544
4545
4546
4547
 
# File 'lib/flt/num.rb', line 4538

def _log_radix_mult
  case num_class.radix
  when 10
    LOG10_MULT
  when 2
    LOG2_MULT
  else
    raise ArgumentError, "_log_radix_mult not implemented for base #{num_class.radix}"
  end
end

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

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



4197
4198
4199
4200
4201
4202
4203
4204
4205
4206
4207
4208
4209
4210
4211
4212
4213
4214
4215
 
# File 'lib/flt/num.rb', line 4197

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 = num_class.new([other.sign, 1, exp])
  end
  tmp = Num(tmp.sign,
                    num_class.int_mult_radix_power(tmp.coefficient, tmp.exponent-other.exponent),
                    other.exponent)
  return swap ? [other, tmp] : [tmp, other]
end

._number_of_digits(v) ⇒ Object 



4560
4561
4562
 
# File 'lib/flt/num.rb', line 4560

def _number_of_digits(v)
  _ndigits(v, num_class.radix)
end

._parser(txt, options = {}) ⇒ Object

Parse numeric text literals (internal use)



4179
4180
4181
4182
4183
4184
4185
4186
4187
4188
4189
4190
4191
4192
4193
4194
 
# File 'lib/flt/num.rb', line 4179

def _parser(txt, options={})
  base = options[:base]
  md = /^\s*([-+])?(?:(?:(\d+)(?:\.(\d*))?|\.(\d+))(?:E([-+]?\d+))?|Inf(?:inity)?|(s)?NaN(\d*))\s*$/i.match(txt)
  if md
    base ||= 10
    OpenStruct.new :sign=>md[1], :int=>md[2], :frac=>md[3], :onlyfrac=>md[4], :exp=>md[5],
                   :signal=>md[6], :diag=>md[7], :base=>base
  else
    md = /^\s*([-+])?0x(?:(?:([\da-f]+)(?:\.([\da-f]*))?|\.([\da-f]+))(?:P([-+]?\d+))?)\s*$/i.match(txt)
    if md
      base = 16
      OpenStruct.new :sign=>md[1], :int=>md[2], :frac=>md[3], :onlyfrac=>md[4], :exp=>md[5],
                     :signal=>nil, :diag=>nil, :base=>base, :exp_base=>2
    end
  end
end

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

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

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

in other words, c*radix**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

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

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



4229
4230
4231
4232
4233
4234
4235
4236
4237
4238
4239
4240
4241
4242
4243
4244
4245
4246
4247
4248
4249
4250
4251
4252
4253
4254
4255
4256
4257
4258
4259
 
# File 'lib/flt/num.rb', line 4229

def _power(xc, xe, yc, ye, p)
  # Find b such that radix**(b-1) <= |y| <= radix**b
  b = _number_of_digits(yc.abs) + ye

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

  # compute product y*log(x) = yc*lxc*radix**(-p-b-1+ye) = pc*radix**(-p-1)
  shift = ye-b
  if shift >= 0
      pc = lxc*yc*num_class.int_radix_power(shift)
  else
      pc = _div_nearest(lxc*yc, num_class.int_radix_power(-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 = num_class.int_radix_power(p-1)+1, 1-p
      else
          coeff, exp = num_class.int_radix_power(p)-1, -p
      end
  else
      coeff, exp = _exp(pc, -(p+1), p+1)
      coeff = _div_nearest(coeff, num_class.radix)
      exp += 1
  end

  return coeff, exp
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.



4460
4461
4462
4463
4464
 
# File 'lib/flt/num.rb', line 4460

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.



4445
4446
4447
4448
4449
4450
4451
4452
4453
4454
4455
4456
 
# File 'lib/flt/num.rb', line 4445

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

.log10_lb(c) ⇒ Object

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

Raises:

  • (ArgumentError) 


4532
4533
4534
4535
4536
 
# File 'lib/flt/num.rb', line 4532

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

.log2_lb(c) ⇒ Object

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

Raises:

  • (ArgumentError) 


4520
4521
4522
4523
4524
 
# File 'lib/flt/num.rb', line 4520

def log2_lb(c)
    raise ArgumentError, "The argument to _log2_lb should be nonnegative." if c <= 0
    str_c = c.to_s(16)
    return LOG2_MULT*4*str_c.length - LOG2_LB_CORRECTION[str_c[0,1].to_i(16)-1]
end