Module: Distribution::MathExtension::Log

Defined in:
lib/distribution/math_extension/log_utilities.rb

Overview

Various logarithm shortcuts, adapted from GSL-1.9.

Constant Summary collapse

C1 = 
-1.quo(2)
C2 = 
1.quo(3)
C3 = 
-1.quo(4)
C4 = 
1.quo(5)
C5 = 
-1.quo(6)
C6 = 
1.quo(7)
C7 = 
-1.quo(8)
C8 = 
1.quo(9)
C9 = 
-1.quo(10)

Class Method Summary collapse

Class Method Details

.log1p(x) ⇒ Object

gsl_log1p from GSL-1.9 sys/log1p.c log for very small x



20
21
22
23
24
25
 
# File 'lib/distribution/math_extension/log_utilities.rb', line 20

def log1p(x)
  # in C, this is volatile double y.
  # Not sure how to reproduce that in Ruby.
  y = 1 + x
  Math.log(y) - ((y - 1) - x).quo(y) # cancel errors with IEEE arithmetic
end

.log_1plusx(x, with_error = false) ⇒ Object

\log(1+x) for x > -1 gsl_sf_log_1plusx_e



29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
# File 'lib/distribution/math_extension/log_utilities.rb', line 29

def log_1plusx(x, with_error = false)
  fail(ArgumentError, 'Range error: x must be > -1') if x <= -1

  if x.abs < Math::ROOT6_FLOAT_EPSILON
    result = x * (1.0 + x * (C1 + x * (C2 + x * (C3 + x * (C4 + x * begin
      C5 + x * (C6 + x * (C7 + x * (C8 + x * C9))) # formerly t = this
    end)))))
    return with_error ? [result, Float::EPSILON * result.abs] : result
  elsif x.abs < 0.5
    c = ChebyshevSeries.evaluate(:lopx, (8 * x + 1).quo(2 * x + 4), with_error)
    return with_error ? [x * c.first, x * c.last] : x * c
  else
    result = Math.log(1 + x)
    return with_error ? [result, Float::EPSILON * result.abs] : result
  end
end

.log_1plusx_minusx(x, with_error = false) ⇒ Object

\log(1+x)-x for x > -1 gsl_sf_log_1plusx_mx_e



48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
# File 'lib/distribution/math_extension/log_utilities.rb', line 48

def log_1plusx_minusx(x, with_error = false)
  fail(ArgumentError, 'Range error: x must be > -1') if x <= -1

  if x.abs < Math::ROOT5_FLOAT_EPSILON
    result = x * x * (C1 + x * (C2 + x * (C3 + x * (C4 + x * begin
      C5 + x * (C6 + x * (C7 + x * (C8 + x * C9))) # formerly t = this
    end))))
    return with_error ? [result, Float::EPSILON * result.abs] : result
  elsif x.abs < 0.5
    c = ChebyshevSeries.evaluate(:lopxmx, (8 * x + 1).quo(2 * x + 4), with_error)
    return with_error ? [x * x * c.first, x * x * c.last] : x * x * c
  else
    lterm = Math.log(1.0 + x)
    error = Float::EPSILON * (lterm.abs + x.abs) if with_error
    result = lterm - x
    return with_error ? [result, error] : result
  end
end