Module: BigMath

Defined in:
lib/bigdecimal/math.rb,
ext/bigdecimal/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.99999999999999999999......e0

Class Method Summary collapse

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

Raises:

  • (ArgumentError)


146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
# File 'lib/bigdecimal/math.rb', line 146

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)


102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# File 'lib/bigdecimal/math.rb', line 102

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)


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

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.



3934
3935
3936
3937
3938
3939
3940
3941
3942
3943
3944
3945
3946
3947
3948
3949
3950
3951
3952
3953
3954
3955
3956
3957
3958
3959
3960
3961
3962
3963
3964
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979
3980
3981
3982
3983
3984
3985
3986
3987
3988
3989
3990
3991
3992
3993
3994
3995
3996
3997
3998
3999
4000
4001
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
4030
4031
4032
4033
4034
4035
4036
4037
4038
4039
4040
4041
4042
4043
4044
4045
4046
4047
4048
4049
4050
4051
4052
4053
# File 'ext/bigdecimal/bigdecimal.c', line 3934

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 almost 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 = BIGDECIMAL_NEGATIVE_P(vx);
	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, 0, 0);
	}
	break;

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

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

    n = prec + BIGDECIMAL_DOUBLE_FIGURES;
    negative = BIGDECIMAL_NEGATIVE_P(vx);
    if (negative) {
        VALUE x_zero = INT2NUM(1);
        VALUE x_copy = f_BigDecimal(1, &x_zero, klass);
        x = BigDecimal_initialize_copy(x_copy, x);
        vx = DATA_PTR(x);
	VpSetSign(vx, 1);
    }

    one = VpCheckGetValue(NewOneWrapLimited(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 < BIGDECIMAL_DOUBLE_FIGURES) {
	    m = BIGDECIMAL_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.



4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080
4081
4082
4083
4084
4085
4086
4087
4088
4089
4090
4091
4092
4093
4094
4095
4096
4097
4098
4099
4100
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112
4113
4114
4115
4116
4117
4118
4119
4120
4121
4122
4123
4124
4125
4126
4127
4128
4129
4130
4131
4132
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
4150
4151
4152
4153
4154
4155
4156
4157
4158
4159
4160
4161
4162
4163
4164
4165
4166
4167
4168
4169
4170
4171
4172
4173
4174
4175
4176
4177
4178
4179
4180
4181
4182
4183
4184
4185
4186
4187
4188
4189
4190
4191
4192
4193
4194
4195
4196
4197
4198
4199
4200
4201
4202
4203
4204
4205
4206
4207
4208
4209
4210
4211
4212
4213
4214
4215
4216
4217
# File 'ext/bigdecimal/bigdecimal.c', line 4067

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 almost 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 = BIGDECIMAL_NEGATIVE_P(vx);
	  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, 0, 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 = NewZeroWrapNolimit(1, prec);
	RB_GC_GUARD(vy->obj);
	VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
        return VpCheckGetValue(vy);
    }
    else if (nan) {
	Real* vy = NewZeroWrapNolimit(1, prec);
	RB_GC_GUARD(vy->obj);
	VpSetNaN(vy);
        return VpCheckGetValue(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 = VpCheckGetValue(vx);

    one = VpCheckGetValue(NewOneWrapLimited(1, 1));
    two = VpCheckGetValue(VpCreateRbObject(1, "2", true));

    n = prec + BIGDECIMAL_DOUBLE_FIGURES;
    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, VpCheckGetValue(VpCreateRbObject(1, buf, true)), vn);
    }
    else {
	expo = 0;
    }
    w = BigDecimal_sub(x, one);
    x = BigDecimal_div2(w, BigDecimal_add(x, one), vn);
    x2 = BigDecimal_mult2(x, x, vn);
    y = x;
    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 < BIGDECIMAL_DOUBLE_FIGURES) {
	    m = BIGDECIMAL_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 = VpCheckGetValue(GetVpValue(SSIZET2NUM(expo), 1));
	dy = BigDecimal_mult(log10, vexpo);
	y = BigDecimal_add(y, dy);
    }

    RB_GC_GUARD(one);
    RB_GC_GUARD(two);
    RB_GC_GUARD(vn);
    RB_GC_GUARD(x2);
    RB_GC_GUARD(y);
    RB_GC_GUARD(d);

    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)


183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# File 'lib/bigdecimal/math.rb', line 183

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)


58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
# File 'lib/bigdecimal/math.rb', line 58

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('2'), 16).to_s
#=> "0.1414213562373095048801688724e1"


43
44
45
# File 'lib/bigdecimal/math.rb', line 43

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