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.



4601
4602
4603
 
# File 'lib/flt/num.rb', line 4601

def log_radix_digits
  @log_radix_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



4288
4289
4290
4291
4292
4293
4294
4295
4296
4297
4298
 
# File 'lib/flt/num.rb', line 4288

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.



4590
4591
4592
4593
 
# File 'lib/flt/num.rb', line 4590

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.



4394
4395
4396
4397
4398
4399
4400
4401
4402
4403
4404
4405
4406
4407
4408
4409
4410
4411
4412
4413
4414
4415
4416
4417
4418
4419
 
# File 'lib/flt/num.rb', line 4394

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



4475
4476
4477
4478
4479
4480
4481
4482
4483
4484
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
 
# File 'lib/flt/num.rb', line 4475

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.



4520
4521
4522
4523
4524
4525
4526
4527
4528
4529
4530
4531
4532
4533
4534
4535
4536
4537
4538
4539
4540
4541
4542
4543
4544
4545
4546
4547
4548
4549
4550
4551
4552
4553
4554
4555
4556
4557
4558
4559
4560
4561
 
# File 'lib/flt/num.rb', line 4520

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.



4425
4426
4427
4428
4429
4430
4431
4432
4433
4434
4435
4436
4437
4438
4439
4440
4441
4442
4443
4444
4445
4446
4447
4448
4449
4450
4451
4452
4453
4454
4455
4456
4457
4458
4459
4460
4461
4462
4463
4464
4465
4466
4467
4468
4469
 
# File 'lib/flt/num.rb', line 4425

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) 


4607
4608
4609
4610
4611
4612
4613
4614
4615
4616
4617
4618
4619
4620
4621
4622
4623
4624
4625
4626
4627
4628
4629
4630
4631
4632
4633
4634
 
# File 'lib/flt/num.rb', line 4607

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 



4671
4672
4673
4674
4675
4676
4677
4678
4679
4680
 
# File 'lib/flt/num.rb', line 4671

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 



4660
4661
4662
4663
4664
4665
4666
4667
4668
4669
 
# File 'lib/flt/num.rb', line 4660

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.



4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
4329
4330
4331
4332
4333
4334
4335
4336
4337
 
# File 'lib/flt/num.rb', line 4319

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 



4682
4683
4684
 
# File 'lib/flt/num.rb', line 4682

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

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

Parse numeric text literals (internal use)



4301
4302
4303
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
 
# File 'lib/flt/num.rb', line 4301

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.



4351
4352
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
 
# File 'lib/flt/num.rb', line 4351

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.



4582
4583
4584
4585
4586
 
# File 'lib/flt/num.rb', line 4582

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.



4567
4568
4569
4570
4571
4572
4573
4574
4575
4576
4577
4578
 
# File 'lib/flt/num.rb', line 4567

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) 


4654
4655
4656
4657
4658
 
# File 'lib/flt/num.rb', line 4654

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) 


4642
4643
4644
4645
4646
 
# File 'lib/flt/num.rb', line 4642

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