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.
log (x, prec)
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"
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

  Computes the arctangent of x to the specified number of digits of precision.

• .cos(x, prec) ⇒ Object

  Computes the cosine of x to the specified number of digits of precision.

• .E(prec) ⇒ Object

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

• .exp ⇒ Object

  BigMath.exp(x, prec).

• .log ⇒ Object

  BigMath.log(x, prec).

• .PI(prec) ⇒ Object

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

• .sin(x, prec) ⇒ Object

  Computes the sine of x to the specified number of digits of precision.

• .sqrt(x, prec) ⇒ Object

  Computes the square root of x to the specified number of digits of precision.

Class Method Details

.atan(x, prec) ⇒ Object

Computes the arctangent of x to the specified number of digits of precision.

If x is NaN, returns NaN.

Raises:

• (ArgumentError)

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

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

Computes the cosine of x to the specified number of digits of precision.

If x is infinite or NaN, returns NaN.

Raises:

• (ArgumentError)

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

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

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

Raises:

• (ArgumentError)

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

def E(prec)
  raise ArgumentError, "Zero or negative precision for E" if prec <= 0
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  y  = one
  d  = y
  z  = one
  i  = 0
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    i += 1
    z *= i
    d  = one.div(z,m)
    y += d
  end
  y
end

.exp ⇒ Object

BigMath.exp(x, prec)

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

If x is infinite, returns Infinity.

If x is NaN, returns NaN.

# File 'bigdecimal.c', line 2414

static VALUE
BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
{
    ssize_t prec, n, i;
    Real* vx = NULL;
    VALUE one, d, x1, y, z;
    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(INT2NUM(0), prec, 1));
	}
	else {
	    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 (vx == NULL) {
	cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
    }
    RB_GC_GUARD(vx->obj);

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

    RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
    RB_GC_GUARD(x1) = one;
    RB_GC_GUARD(y)  = one;
    RB_GC_GUARD(d)  = y;
    RB_GC_GUARD(z)  = one;
    i  = 0;

    while (!VpIsZero((Real*)DATA_PTR(d))) {
	VALUE argv[2];
	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();
	}

	x1 = BigDecimal_mult2(x1, x, SSIZET2NUM(n));
	++i;
	z = BigDecimal_mult(z, SSIZET2NUM(i));
	argv[0] = z;
	argv[1] = SSIZET2NUM(m);
	d = BigDecimal_div2(2, argv, x1);
	y = BigDecimal_add(y, d);
    }

    if (negative) {
	VALUE argv[2];
	argv[0] = y;
	argv[1] = vprec;
	return BigDecimal_div2(2, argv, one);
    }
    else {
	vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
	return BigDecimal_round(1, &vprec, y);
    }
}

.log ⇒ Object

BigMath.log(x, prec)

Computes the natural logarithm of x to the specified number of digits of precision.

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

If x is positive infinite, returns Infinity.

If x is NaN, returns NaN.

# File 'bigdecimal.c', line 2546

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

    if (TYPE(vprec) != T_FIXNUM && TYPE(vprec) != T_BIGNUM) {
	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:
	zero = RBIGNUM_ZERO_P(x);
	negative = RBIGNUM_NEGATIVE_P(x);
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[16];
	snprintf(buf, 16, "1E%ld", -expo);
	x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
    }
    else {
	expo = 0;
    }
    w = BigDecimal_sub(x, one);
    argv[0] = BigDecimal_add(x, one);
    argv[1] = vn;
    x = BigDecimal_div2(2, argv, w);
    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;
	argv[0] = SSIZET2NUM(i);
	argv[1] = SSIZET2NUM(m);
	d = BigDecimal_div2(2, argv, x);
	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

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

Raises:

• (ArgumentError)

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

def PI(prec)
  raise ArgumentError, "Zero or negative argument 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
  w = 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
  w = 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

Computes the sine of x to the specified number of digits of precision.

If x is infinite or NaN, returns NaN.

Raises:

• (ArgumentError)

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

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

Computes the square root of x to the specified number of digits of precision.

BigDecimal.new('2').sqrt(16).to_s

-> "0.14142135623730950488016887242096975E1"

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

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