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.



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

def log_radix_digits
  @log_radix_digits
end

Class Method Details

._convert(x, error = true) ⇒ Object

Convert a numeric value to decimal (internal use)



4331
4332
4333
4334
4335
4336
4337
4338
4339
4340
4341
 
# File 'lib/flt/num.rb', line 4331

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.



4633
4634
4635
4636
 
# File 'lib/flt/num.rb', line 4633

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.



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

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



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

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.



4563
4564
4565
4566
4567
4568
4569
4570
4571
4572
4573
4574
4575
4576
4577
4578
4579
4580
4581
4582
4583
4584
4585
4586
4587
4588
4589
4590
4591
4592
4593
4594
4595
4596
4597
4598
4599
4600
4601
4602
4603
4604
 
# File 'lib/flt/num.rb', line 4563

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.



4468
4469
4470
4471
4472
4473
4474
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
4510
4511
4512
 
# File 'lib/flt/num.rb', line 4468

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) 


4650
4651
4652
4653
4654
4655
4656
4657
4658
4659
4660
4661
4662
4663
4664
4665
4666
4667
4668
4669
4670
4671
4672
4673
4674
4675
4676
4677
 
# File 'lib/flt/num.rb', line 4650

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 



4714
4715
4716
4717
4718
4719
4720
4721
4722
4723
 
# File 'lib/flt/num.rb', line 4714

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 



4703
4704
4705
4706
4707
4708
4709
4710
4711
4712
 
# File 'lib/flt/num.rb', line 4703

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.



4362
4363
4364
4365
4366
4367
4368
4369
4370
4371
4372
4373
4374
4375
4376
4377
4378
4379
4380
 
# File 'lib/flt/num.rb', line 4362

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 



4725
4726
4727
 
# File 'lib/flt/num.rb', line 4725

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

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

Parse numeric text literals (internal use)



4344
4345
4346
4347
4348
4349
4350
4351
4352
4353
4354
4355
4356
4357
4358
4359
 
# File 'lib/flt/num.rb', line 4344

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.



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
4420
4421
4422
4423
4424
 
# File 'lib/flt/num.rb', line 4394

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.



4625
4626
4627
4628
4629
 
# File 'lib/flt/num.rb', line 4625

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.



4610
4611
4612
4613
4614
4615
4616
4617
4618
4619
4620
4621
 
# File 'lib/flt/num.rb', line 4610

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) 


4697
4698
4699
4700
4701
 
# File 'lib/flt/num.rb', line 4697

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) 


4685
4686
4687
4688
4689
 
# File 'lib/flt/num.rb', line 4685

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