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)


256
257
258
259
260
261
262
263
264
265
266
267
268
# File 'lib/bigdecimal/math.rb', line 256

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


446
447
448
449
450
451
452
453
454
# File 'lib/bigdecimal/math.rb', line 446

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


230
231
232
233
234
235
236
237
238
239
240
241
242
243
# File 'lib/bigdecimal/math.rb', line 230

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.double_fig
  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"


424
425
426
427
428
429
430
431
432
433
# File 'lib/bigdecimal/math.rb', line 424

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.double_fig
  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"


281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
# File 'lib/bigdecimal/math.rb', line 281

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.double_fig
  pi = PI(n)
  x = -x if neg = x < 0
  return pi.div(neg ? -2 : 2, prec) if x.infinite?
  return pi.div(neg ? -4 : 4, prec) if x.round(n) == 1
  x = BigDecimal("1").div(x, n) if inv = x > 1
  x = (-1 + sqrt(1 + x.mult(x, n), n)).div(x, n) if dbl = x > 0.5
  y = x
  d = y
  t = x
  r = BigDecimal("3")
  x2 = x.mult(x,n)
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t = -t.mult(x2,n)
    d = t.div(r,m)
    y += d
    r += 2
  end
  y *= 2 if dbl
  y = pi / 2 - y if inv
  y = -y if neg
  y.mult(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"


319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
# File 'lib/bigdecimal/math.rb', line 319

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


467
468
469
470
471
472
473
474
475
476
477
# File 'lib/bigdecimal/math.rb', line 467

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.double_fig
  (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"


106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# File 'lib/bigdecimal/math.rb', line 106

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(x.to_f), 0)
  precs = [prec + BigDecimal.double_fig]
  precs << 2 + precs.last / 2 while precs.last > BigDecimal.double_fig
  precs.reverse_each 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"


192
193
194
195
196
197
198
# File 'lib/bigdecimal/math.rb', line 192

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?
  sign, x = _sin_periodic_reduction(x, prec + BigDecimal.double_fig, add_half_pi: true)
  sign * sin(x, 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"


379
380
381
382
383
384
385
386
387
388
# File 'lib/bigdecimal/math.rb', line 379

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.double_fig
  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"


944
945
946
947
# File 'lib/bigdecimal/math.rb', line 944

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

.erf(x, prec) ⇒ Object

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"


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

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 = x.to_f
    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)
    if erfc
      # Workaround for https://github.com/ruby/bigdecimal/issues/464
      return BigDecimal(1) if erfc.exponent < -prec - BigDecimal.double_fig

      return BigDecimal(1).sub(erfc, prec)
    end
  end

  prec2 = prec + BigDecimal.double_fig
  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"


638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
# File 'lib/bigdecimal/math.rb', line 638

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.double_fig), prec) if x < 0
  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)) extra digits
  log10 = 2.302585092994046
  high_prec = prec + BigDecimal.double_fig + (x.ceil**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.



332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
# File 'lib/bigdecimal.rb', line 332

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.double_fig + cnt
  x = x._decimal_shift(-cnt)

  # Calculation of exp(small_prec) is fast because calculation of x**n is fast
  # Calculation of exp(small_abs) converges fast.
  # exp(x) = exp(small_prec_part + small_abs_part) = exp(small_prec_part) * exp(small_abs_part)
  x_small_prec = x.round(Integer.sqrt(prec2))
  y = _exp_taylor(x_small_prec, prec2).mult(_exp_taylor(x.sub(x_small_prec, prec2), prec2), prec2)

  # 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"


559
560
561
562
563
564
565
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 559

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.double_fig
  elsif x < 1
    exp_prec = prec - x.exponent + BigDecimal.double_fig
  else
    exp_prec = prec
  end

  return BigDecimal(-1) if exp_prec <= 0

  exp = exp(x, exp_prec)

  if exp.exponent > prec + BigDecimal.double_fig
    # Workaroudn for https://github.com/ruby/bigdecimal/issues/464
    exp
  else
    exp.sub(1, prec)
  end
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]


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

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"


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

def gamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
  prec2 = prec + BigDecimal.double_fig
  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"


134
135
136
137
138
139
140
141
142
# File 'lib/bigdecimal/math.rb', line 134

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


885
886
887
888
# File 'lib/bigdecimal/math.rb', line 885

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]


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
798
799
# File 'lib/bigdecimal/math.rb', line 757

def lgamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
  prec2 = prec + BigDecimal.double_fig
  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.double_fig

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


255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
# File 'lib/bigdecimal.rb', line 255

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.double_fig
  BigDecimal.save_limit do
    BigDecimal.limit(0)
    if x > 10 || x < 0.1
      log10 = log(BigDecimal(10), prec2)
      exponent = x.exponent
      x = x._decimal_shift(-exponent)
      if x < 0.3
        x *= 10
        exponent -= 1
      end
      return (log10 * exponent).add(log(x, prec2), prec)
    end

    x_minus_one_exponent = (x - 1).exponent

    # log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
    sqrt_steps = [Integer.sqrt(prec2) + 3 * x_minus_one_exponent, 0].max

    lg2 = 0.3010299956639812
    sqrt_prec = prec2 + [-x_minus_one_exponent, 0].max + (sqrt_steps * lg2).ceil

    sqrt_steps.times do
      x = x.sqrt(sqrt_prec)
    end

    # Taylor series for log(x) around 1
    # log(x) = -log((1 + X) / (1 - X)) where X = (x - 1) / (x + 1)
    # log(x) = 2 * (X + X**3 / 3 + X**5 / 5 + X**7 / 7 + ...)
    x = (x - 1).div(x + 1, sqrt_prec)
    y = x
    x2 = x.mult(x, prec2)
    1.step do |i|
      n = prec2 + x.exponent - y.exponent + x2.exponent
      break if n <= 0 || x.zero?
      x = x.mult(x2.round(n - x2.exponent), n)
      y = y.add(x.div(2 * i + 1, n), prec2)
    end

    y.mult(2 ** (sqrt_steps + 1), prec)
  end
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"


522
523
524
525
526
527
528
529
530
531
532
533
# File 'lib/bigdecimal/math.rb', line 522

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.double_fig * 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.double_fig - (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)


543
544
545
546
547
548
549
# File 'lib/bigdecimal/math.rb', line 543

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"


494
495
496
497
498
499
500
501
502
503
504
505
# File 'lib/bigdecimal/math.rb', line 494

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.double_fig * 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.double_fig - (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"


899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
# File 'lib/bigdecimal/math.rb', line 899

def PI(prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI)
  n      = prec + BigDecimal.double_fig
  zero   = BigDecimal("0")
  one    = BigDecimal("1")
  two    = BigDecimal("2")

  m25    = BigDecimal("-0.04")
  m57121 = BigDecimal("-57121")

  pi     = zero

  d = one
  k = one
  t = BigDecimal("-80")
  while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t   = t*m25
    d   = t.div(k,m)
    k   = k+two
    pi  = pi + d
  end

  d = one
  k = one
  t = BigDecimal("956")
  while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t   = t.div(m57121,n)
    d   = t.div(k,m)
    pi  = pi + d
    k   = k+two
  end
  pi.mult(1, 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"


155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
# File 'lib/bigdecimal/math.rb', line 155

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.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  sign, x = _sin_periodic_reduction(x, n)
  x1   = x
  x2   = x.mult(x,n)
  y    = x
  d    = y
  i    = one
  z    = one
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    x1  = -x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = x1.div(z,m)
    y  += d
  end
  y = BigDecimal("1") if y > 1
  y.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"


356
357
358
359
360
361
362
363
364
365
366
# File 'lib/bigdecimal/math.rb', line 356

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.double_fig
  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"


214
215
216
217
# File 'lib/bigdecimal/math.rb', line 214

def tan(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :tan)
  sin(x, prec + BigDecimal.double_fig).div(cos(x, prec + BigDecimal.double_fig), 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"


401
402
403
404
405
406
407
408
409
410
411
# File 'lib/bigdecimal/math.rb', line 401

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.double_fig + [-x.exponent, 0].max
  e = exp(x, prec2)
  einv = BigDecimal(1).div(e, prec2)
  (e - einv).div(e + einv, prec)
end