Module: BigMath

Defined in:

lib/bigdecimal/math.rb,
bigdecimal.c

Overview

– Contents:

sqrt(x, prec)
sin (x, prec)
cos (x, prec)
atan(x, prec)  Note: |x|<1, x=0.9999 may not converge.
PI  (prec)
E   (prec) == exp(1.0,prec)

where:

x    ... BigDecimal number to be computed.
         |x| must be small enough to get convergence.
prec ... Number of digits to be obtained.

++

Provides mathematical functions.

Example:

require "bigdecimal/math"

include BigMath

a = BigDecimal((PI(100)/2).to_s)
puts sin(a,100) # => 0.10000000000000000000......E1

Class Method Summary

• .atan(x, prec) ⇒ Object

  call-seq: atan(decimal, numeric) -> BigDecimal.

• .cos(x, prec) ⇒ Object

  call-seq: cos(decimal, numeric) -> BigDecimal.

• .E(prec) ⇒ Object

  call-seq: E(numeric) -> BigDecimal.

• .exp(x, vprec) ⇒ Object

  BigMath.exp(decimal, numeric) -> BigDecimal.

• .log(x, vprec) ⇒ Object

  BigMath.log(decimal, numeric) -> BigDecimal.

• .PI(prec) ⇒ Object

  call-seq: PI(numeric) -> BigDecimal.

• .sin(x, prec) ⇒ Object

  call-seq: sin(decimal, numeric) -> BigDecimal.

• .sqrt(x, prec) ⇒ Object

  call-seq: sqrt(decimal, numeric) -> BigDecimal.

Class Method Details

.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.new('-1'), 16).to_s
#=> "-0.785398163397448309615660845819878471907514682065E0"

Raises:

• (ArgumentError)

# File 'lib/bigdecimal/math.rb', line 145

def atan(x, prec)
  raise ArgumentError, "Zero or negative precision for atan" if prec <= 0
  return BigDecimal("NaN") if x.nan?
  pi = PI(prec)
  x = -x if neg = x < 0
  return pi.div(neg ? -2 : 2, prec) if x.infinite?
  return pi / (neg ? -4 : 4) if x.round(prec) == 1
  x = BigDecimal("1").div(x, prec) if inv = x > 1
  x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
  n    = prec + BigDecimal.double_fig
  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
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(4), 16).to_s
#=> "-0.999999999999999999999999999999856613163740061349E0"

Raises:

• (ArgumentError)

# File 'lib/bigdecimal/math.rb', line 101

def cos(x, prec)
  raise ArgumentError, "Zero or negative precision for cos" if prec <= 0
  return BigDecimal("NaN") if x.infinite? || x.nan?
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  x = -x if x < 0
  if x > (twopi = two * BigMath.PI(prec))
    if x > 30
      x %= twopi
    else
      x -= twopi while x > twopi
    end
  end
  x1 = one
  x2 = x.mult(x,n)
  sign = 1
  y = one
  d = y
  i = BigDecimal("0")
  z = one
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    sign = -sign
    x1  = x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = sign * x1.div(z,m)
    y  += d
  end
  y
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(10).to_s
#=> "0.271828182845904523536028752390026306410273E1"

Raises:

• (ArgumentError)

# File 'lib/bigdecimal/math.rb', line 227

def E(prec)
  raise ArgumentError, "Zero or negative precision for E" if prec <= 0
  BigMath.exp(1, prec)
end

.exp(x, vprec) ⇒ Object

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.

# File 'bigdecimal.c', line 2720

static VALUE
BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
{
    ssize_t prec, n, i;
    Real* vx = NULL;
    VALUE one, d, y;
    int negative = 0;
    int infinite = 0;
    int nan = 0;
    double flo;

    prec = NUM2SSIZET(vprec);
    if (prec <= 0) {
	rb_raise(rb_eArgError, "Zero or negative precision for exp");
    }

    /* TODO: the following switch statement is almostly the same as one in the
     *       BigDecimalCmp function. */
    switch (TYPE(x)) {
      case T_DATA:
	if (!is_kind_of_BigDecimal(x)) break;
	vx = DATA_PTR(x);
	negative = VpGetSign(vx) < 0;
	infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
	nan = VpIsNaN(vx);
	break;

      case T_FIXNUM:
	/* fall through */
      case T_BIGNUM:
	vx = GetVpValue(x, 0);
	break;

      case T_FLOAT:
	flo = RFLOAT_VALUE(x);
	negative = flo < 0;
	infinite = isinf(flo);
	nan = isnan(flo);
	if (!infinite && !nan) {
	    vx = GetVpValueWithPrec(x, DBL_DIG+1, 0);
	}
	break;

      case T_RATIONAL:
	vx = GetVpValueWithPrec(x, prec, 0);
	break;

      default:
	break;
    }
    if (infinite) {
	if (negative) {
	    return ToValue(GetVpValueWithPrec(INT2FIX(0), prec, 1));
	}
	else {
	    Real* vy;
	    vy = VpCreateRbObject(prec, "#0");
	    VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
	    RB_GC_GUARD(vy->obj);
	    return ToValue(vy);
	}
    }
    else if (nan) {
	Real* vy;
	vy = VpCreateRbObject(prec, "#0");
	VpSetNaN(vy);
	RB_GC_GUARD(vy->obj);
	return ToValue(vy);
    }
    else if (vx == NULL) {
	cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
    }
    x = vx->obj;

    n = prec + rmpd_double_figures();
    negative = VpGetSign(vx) < 0;
    if (negative) {
	VpSetSign(vx, 1);
    }

    one = ToValue(VpCreateRbObject(1, "1"));
    y   = one;
    d   = y;
    i   = 1;

    while (!VpIsZero((Real*)DATA_PTR(d))) {
	SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
	SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
	ssize_t m = n - vabs(ey - ed);

	rb_thread_check_ints();

	if (m <= 0) {
	    break;
	}
	else if ((size_t)m < rmpd_double_figures()) {
	    m = rmpd_double_figures();
	}

	d = BigDecimal_mult(d, x);                             /* d <- d * x */
	d = BigDecimal_div2(d, SSIZET2NUM(i), SSIZET2NUM(m));  /* d <- d / i */
	y = BigDecimal_add(y, d);                              /* y <- y + d  */
	++i;                                                   /* i  <- i + 1 */
    }

    if (negative) {
	return BigDecimal_div2(one, y, vprec);
    }
    else {
	vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
	return BigDecimal_round(1, &vprec, y);
    }

    RB_GC_GUARD(one);
    RB_GC_GUARD(x);
    RB_GC_GUARD(y);
    RB_GC_GUARD(d);
}

.log(x, vprec) ⇒ Object

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.

# File 'bigdecimal.c', line 2851

static VALUE
BigMath_s_log(VALUE klass, VALUE x, VALUE vprec)
{
    ssize_t prec, n, i;
    SIGNED_VALUE expo;
    Real* vx = NULL;
    VALUE vn, one, two, w, x2, y, d;
    int zero = 0;
    int negative = 0;
    int infinite = 0;
    int nan = 0;
    double flo;
    long fix;

    if (!is_integer(vprec)) {
	rb_raise(rb_eArgError, "precision must be an Integer");
    }

    prec = NUM2SSIZET(vprec);
    if (prec <= 0) {
	rb_raise(rb_eArgError, "Zero or negative precision for exp");
    }

    /* TODO: the following switch statement is almostly the same as one in the
     *       BigDecimalCmp function. */
    switch (TYPE(x)) {
      case T_DATA:
	  if (!is_kind_of_BigDecimal(x)) break;
	  vx = DATA_PTR(x);
	  zero = VpIsZero(vx);
	  negative = VpGetSign(vx) < 0;
	  infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
	  nan = VpIsNaN(vx);
	  break;

      case T_FIXNUM:
	fix = FIX2LONG(x);
	zero = fix == 0;
	negative = fix < 0;
	goto get_vp_value;

      case T_BIGNUM:
        i = FIX2INT(rb_big_cmp(x, INT2FIX(0)));
	zero = i == 0;
	negative = i < 0;
get_vp_value:
	if (zero || negative) break;
	vx = GetVpValue(x, 0);
	break;

      case T_FLOAT:
	flo = RFLOAT_VALUE(x);
	zero = flo == 0;
	negative = flo < 0;
	infinite = isinf(flo);
	nan = isnan(flo);
	if (!zero && !negative && !infinite && !nan) {
	    vx = GetVpValueWithPrec(x, DBL_DIG+1, 1);
	}
	break;

      case T_RATIONAL:
	zero = RRATIONAL_ZERO_P(x);
	negative = RRATIONAL_NEGATIVE_P(x);
	if (zero || negative) break;
	vx = GetVpValueWithPrec(x, prec, 1);
	break;

      case T_COMPLEX:
	rb_raise(rb_eMathDomainError,
		 "Complex argument for BigMath.log");

      default:
	break;
    }
    if (infinite && !negative) {
	Real* vy;
	vy = VpCreateRbObject(prec, "#0");
	RB_GC_GUARD(vy->obj);
	VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
	return ToValue(vy);
    }
    else if (nan) {
	Real* vy;
	vy = VpCreateRbObject(prec, "#0");
	RB_GC_GUARD(vy->obj);
	VpSetNaN(vy);
	return ToValue(vy);
    }
    else if (zero || negative) {
	rb_raise(rb_eMathDomainError,
		 "Zero or negative argument for log");
    }
    else if (vx == NULL) {
	cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
    }
    x = ToValue(vx);

    RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
    RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2"));

    n = prec + rmpd_double_figures();
    RB_GC_GUARD(vn) = SSIZET2NUM(n);
    expo = VpExponent10(vx);
    if (expo < 0 || expo >= 3) {
	char buf[DECIMAL_SIZE_OF_BITS(SIZEOF_VALUE * CHAR_BIT) + 4];
	snprintf(buf, sizeof(buf), "1E%"PRIdVALUE, -expo);
	x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
    }
    else {
	expo = 0;
    }
    w = BigDecimal_sub(x, one);
    x = BigDecimal_div2(w, BigDecimal_add(x, one), vn);
    RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn);
    RB_GC_GUARD(y)  = x;
    RB_GC_GUARD(d)  = y;
    i = 1;
    while (!VpIsZero((Real*)DATA_PTR(d))) {
	SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
	SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
	ssize_t m = n - vabs(ey - ed);
	if (m <= 0) {
	    break;
	}
	else if ((size_t)m < rmpd_double_figures()) {
	    m = rmpd_double_figures();
	}

	x = BigDecimal_mult2(x2, x, vn);
	i += 2;
	d = BigDecimal_div2(x, SSIZET2NUM(i), SSIZET2NUM(m));
	y = BigDecimal_add(y, d);
    }

    y = BigDecimal_mult(y, two);
    if (expo != 0) {
	VALUE log10, vexpo, dy;
	log10 = BigMath_s_log(klass, INT2FIX(10), vprec);
	vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1));
	dy = BigDecimal_mult(log10, vexpo);
	y = BigDecimal_add(y, dy);
    }

    return y;
}

.PI(prec) ⇒ Object

call-seq:

PI(numeric) -> BigDecimal

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

BigMath.PI(10).to_s
#=> "0.3141592653589793238462643388813853786957412E1"

Raises:

• (ArgumentError)

# File 'lib/bigdecimal/math.rb', line 182

def PI(prec)
  raise ArgumentError, "Zero or negative precision for PI" if prec <= 0
  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
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, 5).to_s
#=> "0.70710678118654752440082036563292800375E0"

Raises:

• (ArgumentError)

# File 'lib/bigdecimal/math.rb', line 57

def sin(x, prec)
  raise ArgumentError, "Zero or negative precision for sin" if prec <= 0
  return BigDecimal("NaN") if x.infinite? || x.nan?
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  x = -x if neg = x < 0
  if x > (twopi = two * BigMath.PI(prec))
    if x > 30
      x %= twopi
    else
      x -= twopi while x > twopi
    end
  end
  x1   = x
  x2   = x.mult(x,n)
  sign = 1
  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
    sign = -sign
    x1  = x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = sign * x1.div(z,m)
    y  += d
  end
  neg ? -y : y
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.new('2'), 16).to_s
#=> "0.1414213562373095048801688724E1"

# File 'lib/bigdecimal/math.rb', line 42

def sqrt(x, prec)
  x.sqrt(prec)
end