Module: BigMath

Defined in:
lib/bigdecimal.rb,
lib/bigdecimal/math.rb

Overview

– Contents:

sqrt(x, prec)
cbrt(x, prec)
hypot(x, y, prec)
sin (x, prec)
cos (x, prec)
tan (x, prec)
asin(x, prec)
acos(x, prec)
atan(x, prec)
atan2(y, x, prec)
sinh (x, prec)
cosh (x, prec)
tanh (x, prec)
asinh(x, prec)
acosh(x, prec)
atanh(x, prec)
log2 (x, prec)
log10(x, prec)
log1p(x, prec)
expm1(x, prec)
erf (x, prec)
erfc(x, prec)
gamma(x, prec)
lgamma(x, prec)
frexp(x)
ldexp(x, exponent)
PI  (prec)
E   (prec) == exp(1.0,prec)

where:

x, y ... BigDecimal number to be computed.
prec ... Number of digits to be obtained.

++

Provides mathematical functions.

Example:

require "bigdecimal/math"

include BigMath

a = BigDecimal((PI(49)/2).to_s)
puts sin(a,100) # => 0.9999999999...9999999986e0

Class Method Summary collapse

Class Method Details

.acos(x, prec) ⇒ Object

call-seq:

acos(decimal, numeric) -> BigDecimal

Computes the arccosine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.acos(BigDecimal('0.5'), 32).to_s
#=> "0.10471975511965977461542144610932e1"

Raises:

  • (Math::DomainError) 


270
271
272
273
274
275
276
277
278
279
280
281
282
 
# File 'lib/bigdecimal/math.rb', line 270

def acos(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :acos)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acos)
  raise Math::DomainError, "Out of domain argument for acos" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  return (PI(prec2) / 2).sub(asin(x, prec2), prec) if x < 0
  return PI(prec2).div(2, prec) if x.zero?

  sin = (1 - x**2).sqrt(prec2)
  atan(sin.div(x, prec2), prec)
end

.acosh(x, prec) ⇒ Object

call-seq:

acosh(decimal, numeric) -> BigDecimal

Computes the inverse hyperbolic cosine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.acosh(BigDecimal('2'), 32).to_s
#=> "0.1316957896924816708625046347308e1"

Raises:

  • (Math::DomainError) 


453
454
455
456
457
458
459
460
461
 
# File 'lib/bigdecimal/math.rb', line 453

def acosh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :acosh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acosh)
  raise Math::DomainError, "Out of domain argument for acosh" if x < 1
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  return BigDecimal::Internal.nan_computation_result if x.nan?

  log(x + sqrt(x**2 - 1, prec + BigDecimal::Internal::EXTRA_PREC), prec)
end

.asin(x, prec) ⇒ Object

call-seq:

asin(decimal, numeric) -> BigDecimal

Computes the arcsine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.asin(BigDecimal('0.5'), 32).to_s
#=> "0.52359877559829887307710723054658e0"

Raises:

  • (Math::DomainError) 


244
245
246
247
248
249
250
251
252
253
254
255
256
257
 
# File 'lib/bigdecimal/math.rb', line 244

def asin(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :asin)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asin)
  raise Math::DomainError, "Out of domain argument for asin" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  cos = (1 - x**2).sqrt(prec2)
  if cos.zero?
    PI(prec2).div(x > 0 ? 2 : -2, prec)
  else
    atan(x.div(cos, prec2), prec)
  end
end

.asinh(x, prec) ⇒ Object

call-seq:

asinh(decimal, numeric) -> BigDecimal

Computes the inverse hyperbolic sine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.asinh(BigDecimal('1'), 32).to_s
#=> "0.88137358701954302523260932497979e0"



431
432
433
434
435
436
437
438
439
440
 
# File 'lib/bigdecimal/math.rb', line 431

def asinh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :asinh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asinh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
  return -asinh(-x, prec) if x < 0

  sqrt_prec = prec + [-x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC
  log(x + sqrt(x**2 + 1, sqrt_prec), prec)
end

.atan(x, prec) ⇒ Object

call-seq:

atan(decimal, numeric) -> BigDecimal

Computes the arctangent of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.atan(BigDecimal('-1'), 32).to_s
#=> "-0.78539816339744830961566084581988e0"



295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
 
# File 'lib/bigdecimal/math.rb', line 295

def atan(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  n = prec + BigDecimal::Internal::EXTRA_PREC
  return PI(n).div(2 * x.infinite?, prec) if x.infinite?

  x = -x if neg = x < 0
  x = BigDecimal(1).div(x, n) if inv = x < -1 || x > 1

  # Solve tan(y) - x = 0 with Newton's method
  # Repeat: y -= (tan(y) - x) * cos(y)**2
  y = BigDecimal(Math.atan(BigDecimal::Internal.fast_to_f(x)), 0)
  BigDecimal::Internal.newton_loop(n) do |p|
    s = sin(y, p)
    c = (1 - s * s).sqrt(p)
    y = y.sub(c * (s.sub(c * x.mult(1, p), p)), p)
  end
  y = PI(n) / 2 - y if inv
  y.mult(neg ? -1 : 1, prec)
end

.atan2(y, x, prec) ⇒ Object

call-seq:

atan2(decimal, decimal, numeric) -> BigDecimal

Computes the arctangent of y and x to the specified number of digits of precision, numeric.

BigMath.atan2(BigDecimal('-1'), BigDecimal('1'), 32).to_s
#=> "-0.78539816339744830961566084581988e0"



326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
 
# File 'lib/bigdecimal/math.rb', line 326

def atan2(y, x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan2)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan2)
  y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :atan2)
  return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan?

  if x.infinite? || y.infinite?
    one = BigDecimal(1)
    zero = BigDecimal(0)
    x = x.infinite? ? (x > 0 ? one : -one) : zero
    y = y.infinite? ? (y > 0 ? one : -one) : y.sign * zero
  end

  return x.sign >= 0 ? BigDecimal(0) : y.sign * PI(prec) if y.zero?

  y = -y if neg = y < 0
  xlarge = y.abs < x.abs
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x > 0
    v = xlarge ? atan(y.div(x, prec2), prec) : PI(prec2) / 2 - atan(x.div(y, prec2), prec2)
  else
    v = xlarge ? PI(prec2) - atan(-y.div(x, prec2), prec2) : PI(prec2) / 2 + atan(x.div(-y, prec2), prec2)
  end
  v.mult(neg ? -1 : 1, prec)
end

.atanh(x, prec) ⇒ Object

call-seq:

atanh(decimal, numeric) -> BigDecimal

Computes the inverse hyperbolic tangent of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.atanh(BigDecimal('0.5'), 32).to_s
#=> "0.54930614433405484569762261846126e0"

Raises:

  • (Math::DomainError) 


474
475
476
477
478
479
480
481
482
483
484
 
# File 'lib/bigdecimal/math.rb', line 474

def atanh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atanh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atanh)
  raise Math::DomainError, "Out of domain argument for atanh" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x == 1
  return -BigDecimal::Internal.infinity_computation_result if x == -1

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  (log(x + 1, prec2) - log(1 - x, prec2)).div(2, prec)
end

.cbrt(x, prec) ⇒ Object

call-seq:

cbrt(decimal, numeric) -> BigDecimal

Computes the cube root of decimal to the specified number of digits of precision, numeric.

BigMath.cbrt(BigDecimal('2'), 32).to_s
#=> "0.12599210498948731647672106072782e1"



137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
 
# File 'lib/bigdecimal/math.rb', line 137

def cbrt(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cbrt)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cbrt)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
  return BigDecimal(0) if x.zero?

  x = -x if neg = x < 0
  ex = x.exponent / 3
  x = x._decimal_shift(-3 * ex)
  y = BigDecimal(Math.cbrt(BigDecimal::Internal.fast_to_f(x)), 0)
  BigDecimal::Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p|
    y = (2 * y + x.div(y, p).div(y, p)).div(3, p)
  end
  y._decimal_shift(ex).mult(neg ? -1 : 1, prec)
end

.cos(x, prec) ⇒ Object

call-seq:

cos(decimal, numeric) -> BigDecimal

Computes the cosine of decimal to the specified number of digits of precision, numeric.

If decimal is Infinity or NaN, returns NaN.

BigMath.cos(BigMath.PI(16), 32).to_s
#=> "-0.99999999999999999999999999999997e0"



204
205
206
207
208
209
210
211
 
# File 'lib/bigdecimal/math.rb', line 204

def cos(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cos)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cos)
  return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
  n = prec + BigDecimal::Internal::EXTRA_PREC
  sign, x = _sin_periodic_reduction(x, n, add_half_pi: true)
  _sin_around_zero(x, n).mult(sign, prec)
end

.cosh(x, prec) ⇒ Object

call-seq:

cosh(decimal, numeric) -> BigDecimal

Computes the hyperbolic cosine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.cosh(BigDecimal('1'), 32).to_s
#=> "0.15430806348152437784779056207571e1"



386
387
388
389
390
391
392
393
394
395
 
# File 'lib/bigdecimal/math.rb', line 386

def cosh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cosh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cosh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite?

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  e = exp(x, prec2)
  (e + BigDecimal(1).div(e, prec2)).div(2, prec)
end

.E(prec) ⇒ Object

call-seq:

E(numeric) -> BigDecimal

Computes e (the base of natural logarithms) to the specified number of digits of precision, numeric.

BigMath.E(32).to_s
#=> "0.27182818284590452353602874713527e1"



923
924
925
926
 
# File 'lib/bigdecimal/math.rb', line 923

def E(prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :E)
  exp(1, prec)
end

.erf(x, prec) ⇒ Object

call-seq:

erf(decimal, numeric) -> BigDecimal

Computes the error function of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.erf(BigDecimal('1'), 32).to_s
#=> "0.84270079294971486934122063508261e0"



598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
 
# File 'lib/bigdecimal/math.rb', line 598

def erf(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :erf)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erf)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(x.infinite?) if x.infinite?
  return BigDecimal(0) if x == 0
  return -erf(-x, prec) if x < 0
  return BigDecimal(1) if x > 5000000000 # erf(5000000000) > 1 - 1e-10000000000000000000

  if x > 8
    xf = BigDecimal::Internal.fast_to_f(x)
    log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5)
    erfc_prec = [prec + log10_erfc.ceil, 1].max
    erfc = _erfc_asymptotic(x, erfc_prec)
    return BigDecimal(1).sub(erfc, prec) if erfc
  end

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  x_smallprec = x.mult(1, Integer.sqrt(prec2) / 2)
  # Taylor series of x with small precision is fast
  erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec2)
  # Taylor series converges quickly for small x
  _erf_taylor(x - x_smallprec, x_smallprec, erf1, prec2).mult(1, prec)
end

.erfc(x, prec) ⇒ Object

call-seq:

erfc(decimal, numeric) -> BigDecimal

Computes the complementary error function of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.erfc(BigDecimal('10'), 32).to_s
#=> "0.20884875837625447570007862949578e-44"



634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
 
# File 'lib/bigdecimal/math.rb', line 634

def erfc(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :erfc)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erfc)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(1 - x.infinite?) if x.infinite?
  return BigDecimal(1).sub(erf(x, prec + BigDecimal::Internal::EXTRA_PREC), prec) if x < 0.5
  return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)

  if x >= 8
    y = _erfc_asymptotic(x, prec)
    return y.mult(1, prec) if y
  end

  # erfc(x) = 1 - erf(x) < exp(-x**2)/x/sqrt(pi)
  # Precision of erf(x) needs about log10(exp(-x**2)/x/sqrt(pi)) extra digits
  log10 = 2.302585092994046
  xf = BigDecimal::Internal.fast_to_f(x)
  high_prec = prec + BigDecimal::Internal::EXTRA_PREC + ((xf**2 + Math.log(xf) + Math.log(Math::PI)/2) / log10).ceil
  BigDecimal(1).sub(erf(x, high_prec), prec)
end

.exp(x, prec) ⇒ Object

call-seq:

BigMath.exp(decimal, numeric)    -> BigDecimal

Computes the value of e (the base of natural logarithms) raised to the power of decimal, to the specified number of digits of precision.

If decimal is infinity, returns Infinity.

If decimal is NaN, returns NaN.



358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
 
# File 'lib/bigdecimal.rb', line 358

def exp(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :exp)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :exp)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return x.positive? ? BigDecimal::Internal.infinity_computation_result : BigDecimal(0) if x.infinite?
  return BigDecimal(1) if x.zero?

  # exp(x * 10**cnt) = exp(x)**(10**cnt)
  cnt = x < -1 || x > 1 ? x.exponent : 0
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC + cnt
  x = x._decimal_shift(-cnt)

  # Decimal form of bit-burst algorithm
  # Calculate exp(x.xxxxxxxxxxxxxxxx) as
  # exp(x.xx) * exp(0.00xx) * exp(0.0000xxxx) * exp(0.00000000xxxxxxxx)
  x = x.mult(1, prec2)
  n = 2
  y = BigDecimal(1)
  BigDecimal.save_limit do
    BigDecimal.limit(0)
    while x != 0 do
      partial_x = x.truncate(n)
      x -= partial_x
      y = y.mult(_exp_binary_splitting(partial_x, prec2), prec2)
      n *= 2
    end
  end

  # calculate exp(x * 10**cnt) from exp(x)
  # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10
  cnt.times do
    y2 = y.mult(y, prec2)
    y5 = y2.mult(y2, prec2).mult(y, prec2)
    y = y5.mult(y5, prec2)
  end

  y.mult(1, prec)
end

.expm1(x, prec) ⇒ Object

call-seq:

BigMath.expm1(decimal, numeric)    -> BigDecimal

Computes exp(decimal) - 1 to the specified number of digits of precision, numeric.

BigMath.expm1(BigDecimal('0.1'), 32).to_s
#=> "0.10517091807564762481170782649025e0"



566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
 
# File 'lib/bigdecimal/math.rb', line 566

def expm1(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :expm1)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :expm1)
  return BigDecimal(-1) if x.infinite? == -1

  exp_prec = prec
  if x < -1
    # log10(exp(x)) = x * log10(e)
    lg_e = 0.4342944819032518
    exp_prec = prec + (lg_e * x).ceil + BigDecimal::Internal::EXTRA_PREC
  elsif x < 1
    exp_prec = prec - x.exponent + BigDecimal::Internal::EXTRA_PREC
  else
    exp_prec = prec
  end

  return BigDecimal(-1) if exp_prec <= 0

  exp(x, exp_prec).sub(1, prec)
end

.frexp(x) ⇒ Object

call-seq:

frexp(x) -> [BigDecimal, Integer]

Decomposes x into a normalized fraction and an integral power of ten.

BigMath.frexp(BigDecimal(123.456))
#=> [0.123456e0, 3]



866
867
868
869
870
871
872
 
# File 'lib/bigdecimal/math.rb', line 866

def frexp(x)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :frexp)
  return [x, 0] unless x.finite?

  exponent = x.exponent
  [x._decimal_shift(-exponent), exponent]
end

.gamma(x, prec) ⇒ Object

call-seq:

BigMath.gamma(decimal, numeric)    -> BigDecimal

Computes the gamma function of decimal to the specified number of digits of precision, numeric.

BigMath.gamma(BigDecimal('0.5'), 32).to_s
#=> "0.17724538509055160272981674833411e1"



727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
 
# File 'lib/bigdecimal/math.rb', line 727

def gamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x < 0.5
    raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero?

    # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
    pi = PI(prec2)
    sin = _sinpix(x, pi, prec2)
    return pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec)
  elsif x.frac.zero? && x < 1000 * prec
    return _gamma_positive_integer(x, prec2).mult(1, prec)
  end

  a, sum = _gamma_spouge_sum_part(x, prec2)
  (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec)
end

.hypot(x, y, prec) ⇒ Object

call-seq:

hypot(x, y, numeric) -> BigDecimal

Returns sqrt(x**2 + y**2) to the specified number of digits of precision, numeric.

BigMath.hypot(BigDecimal('1'), BigDecimal('2'), 32).to_s
#=> "0.22360679774997896964091736687313e1"



163
164
165
166
167
168
169
170
171
 
# File 'lib/bigdecimal/math.rb', line 163

def hypot(x, y, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :hypot)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :hypot)
  y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :hypot)
  return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? || y.infinite?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  sqrt(x.mult(x, prec2) + y.mult(y, prec2), prec)
end

.ldexp(x, exponent) ⇒ Object

call-seq:

ldexp(fraction, exponent) -> BigDecimal

Inverse of frexp. Returns the value of fraction * 10**exponent.

BigMath.ldexp(BigDecimal("0.123456e0"), 3)
#=> 0.123456e3



883
884
885
886
 
# File 'lib/bigdecimal/math.rb', line 883

def ldexp(x, exponent)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :ldexp)
  x.finite? ? x._decimal_shift(exponent) : x
end

.lgamma(x, prec) ⇒ Object

call-seq:

BigMath.lgamma(decimal, numeric)    -> [BigDecimal, Integer]

Computes the natural logarithm of the absolute value of the gamma function of decimal to the specified number of digits of precision, numeric and its sign.

BigMath.lgamma(BigDecimal('0.5'), 32)
#=> [0.57236494292470008707171367567653e0, 1]



755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
 
# File 'lib/bigdecimal/math.rb', line 755

def lgamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x < 0.5
    return [BigDecimal::INFINITY, 1] if x.frac.zero?

    loop do
      # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
      pi = PI(prec2)
      sin = _sinpix(x, pi, prec2)
      log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec)
      return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC

      # Retry with higher precision if loss of significance is too large
      prec2 = prec2 * 3 / 2
    end
  elsif x.frac.zero? && x < 1000 * prec
    log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), prec)
    [log_gamma, 1]
  else
    # if x is close to 1 or 2, increase precision to reduce loss of significance
    diff1_exponent = (x - 1).exponent
    diff2_exponent = (x - 2).exponent
    extremely_near_one = diff1_exponent < -prec2
    extremely_near_two = diff2_exponent < -prec2

    if extremely_near_one || extremely_near_two
      # If x is extreamely close to base = 1 or 2, linear interpolation is accurate enough.
      # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ...
      # And we can ignore (x - base) ** 2 and higher order terms.
      base = extremely_near_one ? 1 : 2
      d = BigDecimal(1)._decimal_shift(1 - prec2)
      log_gamma_d, sign = lgamma(base + d, prec2)
      return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign]
    end

    prec2 += [-diff1_exponent, -diff2_exponent, 0].max
    a, sum = _gamma_spouge_sum_part(x, prec2)
    log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec)
    [log_gamma, 1]
  end
end

.log(x, prec) ⇒ Object

call-seq:

BigMath.log(decimal, numeric)    -> BigDecimal

Computes the natural logarithm of decimal to the specified number of digits of precision, numeric.

If decimal is zero or negative, raises Math::DomainError.

If decimal is positive infinity, returns Infinity.

If decimal is NaN, returns NaN.

Raises:

  • (Math::DomainError) 


305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
 
# File 'lib/bigdecimal.rb', line 305

def log(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log)
  raise Math::DomainError, 'Complex argument for BigMath.log' if Complex === x

  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  raise Math::DomainError, 'Negative argument for log' if x < 0
  return -BigDecimal::Internal.infinity_computation_result if x.zero?
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  return BigDecimal(0) if x == 1

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC

  # Reduce x to near 1
  if x > 1.01 || x < 0.99
    # log(x) = log(x/exp(logx_approx)) + logx_approx
    logx_approx = BigDecimal(BigDecimal::Internal.float_log(x), 0)
    x = x.div(exp(logx_approx, prec2), prec2)
  else
    logx_approx = BigDecimal(0)
  end

  # Solve exp(y) - x = 0 with Newton's method
  # Repeat: y -= (exp(y) - x) / exp(y)
  y = BigDecimal(BigDecimal::Internal.float_log(x), 0)
  exp_additional_prec = [-(x - 1).exponent, 0].max
  BigDecimal::Internal.newton_loop(prec2) do |p|
    expy = exp(y, p + exp_additional_prec)
    y = y.sub(expy.sub(x, p).div(expy, p), p)
  end
  y.add(logx_approx, prec)
end

.log10(x, prec) ⇒ Object

call-seq:

BigMath.log10(decimal, numeric)    -> BigDecimal

Computes the base 10 logarithm of decimal to the specified number of digits of precision, numeric.

If decimal is zero or negative, raises Math::DomainError.

If decimal is positive infinity, returns Infinity.

If decimal is NaN, returns NaN.

BigMath.log10(BigDecimal('3'), 32).to_s
#=> "0.47712125471966243729502790325512e0"



529
530
531
532
533
534
535
536
537
538
539
540
 
# File 'lib/bigdecimal/math.rb', line 529

def log10(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log10)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log10)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
  v = log(x, prec2).div(log(BigDecimal(10), prec2), prec2)
  # Perform half-up rounding to calculate log10(10**n)==n correctly in every rounding mode
  v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
  v.mult(1, prec)
end

.log1p(x, prec) ⇒ Object

call-seq:

BigMath.log1p(decimal, numeric)    -> BigDecimal

Computes log(1 + decimal) to the specified number of digits of precision, numeric.

BigMath.log1p(BigDecimal('0.1'), 32).to_s
#=> "0.95310179804324860043952123280765e-1"

Raises:

  • (Math::DomainError) 


550
551
552
553
554
555
556
 
# File 'lib/bigdecimal/math.rb', line 550

def log1p(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log1p)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log1p)
  raise Math::DomainError, 'Out of domain argument for log1p' if x < -1

  return log(x + 1, prec)
end

.log2(x, prec) ⇒ Object

call-seq:

BigMath.log2(decimal, numeric)    -> BigDecimal

Computes the base 2 logarithm of decimal to the specified number of digits of precision, numeric.

If decimal is zero or negative, raises Math::DomainError.

If decimal is positive infinity, returns Infinity.

If decimal is NaN, returns NaN.

BigMath.log2(BigDecimal('3'), 32).to_s
#=> "0.15849625007211561814537389439478e1"



501
502
503
504
505
506
507
508
509
510
511
512
 
# File 'lib/bigdecimal/math.rb', line 501

def log2(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log2)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log2)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
  v = log(x, prec2).div(log(BigDecimal(2), prec2), prec2)
  # Perform half-up rounding to calculate log2(2**n)==n correctly in every rounding mode
  v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
  v.mult(1, prec)
end

.PI(prec) ⇒ Object

call-seq:

PI(numeric) -> BigDecimal

Computes the value of pi to the specified number of digits of precision, numeric.

BigMath.PI(32).to_s
#=> "0.31415926535897932384626433832795e1"



897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
 
# File 'lib/bigdecimal/math.rb', line 897

def PI(prec)
  # Gauss–Legendre algorithm
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI)
  n = prec + BigDecimal::Internal::EXTRA_PREC
  a = BigDecimal(1)
  b = BigDecimal(0.5, 0).sqrt(n)
  s = BigDecimal(0.25, 0)
  t = 1
  while a != b && (a - b).exponent > 1 - n
    c = (a - b).div(2, n)
    a, b = (a + b).div(2, n), (a * b).sqrt(n)
    s = s.sub(c * c * t, n)
    t *= 2
  end
  (a * b).div(s, prec)
end

.sin(x, prec) ⇒ Object

call-seq:

sin(decimal, numeric) -> BigDecimal

Computes the sine of decimal to the specified number of digits of precision, numeric.

If decimal is Infinity or NaN, returns NaN.

BigMath.sin(BigMath.PI(5)/4, 32).to_s
#=> "0.70710807985947359435812921837984e0"



184
185
186
187
188
189
190
191
 
# File 'lib/bigdecimal/math.rb', line 184

def sin(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sin)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sin)
  return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
  n = prec + BigDecimal::Internal::EXTRA_PREC
  sign, x = _sin_periodic_reduction(x, n)
  _sin_around_zero(x, n).mult(sign, prec)
end

.sinh(x, prec) ⇒ Object

call-seq:

sinh(decimal, numeric) -> BigDecimal

Computes the hyperbolic sine of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.sinh(BigDecimal('1'), 32).to_s
#=> "0.11752011936438014568823818505956e1"



363
364
365
366
367
368
369
370
371
372
373
 
# File 'lib/bigdecimal/math.rb', line 363

def sinh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sinh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sinh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  prec2 -= x.exponent if x.exponent < 0
  e = exp(x, prec2)
  (e - BigDecimal(1).div(e, prec2)).div(2, prec)
end

.sqrt(x, prec) ⇒ Object

call-seq:

sqrt(decimal, numeric) -> BigDecimal

Computes the square root of decimal to the specified number of digits of precision, numeric.

BigMath.sqrt(BigDecimal('2'), 32).to_s
#=> "0.14142135623730950488016887242097e1"



64
65
66
67
68
 
# File 'lib/bigdecimal/math.rb', line 64

def sqrt(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sqrt)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sqrt)
  x.sqrt(prec)
end

.tan(x, prec) ⇒ Object

call-seq:

tan(decimal, numeric) -> BigDecimal

Computes the tangent of decimal to the specified number of digits of precision, numeric.

If decimal is Infinity or NaN, returns NaN.

BigMath.tan(BigDecimal("0.0"), 4).to_s
#=> "0.0"

BigMath.tan(BigMath.PI(24) / 4, 32).to_s
#=> "0.99999999999999999999999830836025e0"



227
228
229
230
231
 
# File 'lib/bigdecimal/math.rb', line 227

def tan(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :tan)
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  sin(x, prec2).div(cos(x, prec2), prec)
end

.tanh(x, prec) ⇒ Object

call-seq:

tanh(decimal, numeric) -> BigDecimal

Computes the hyperbolic tangent of decimal to the specified number of digits of precision, numeric.

If decimal is NaN, returns NaN.

BigMath.tanh(BigDecimal('1'), 32).to_s
#=> "0.76159415595576488811945828260479e0"



408
409
410
411
412
413
414
415
416
417
418
 
# File 'lib/bigdecimal/math.rb', line 408

def tanh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :tanh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :tanh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(x.infinite?) if x.infinite?

  prec2 = prec + BigDecimal::Internal::EXTRA_PREC + [-x.exponent, 0].max
  e = exp(x, prec2)
  einv = BigDecimal(1).div(e, prec2)
  (e - einv).div(e + einv, prec)
end