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

Raises:

  • (ArgumentError)

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

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

227
228
229
230
# 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

.expObject

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.


2696
2697
2698
2699
2700
2701
2702
2703
2704
2705
2706
2707
2708
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778
2779
2780
2781
2782
2783
2784
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803
2804
2805
2806
2807
2808
2809
2810
2811
2812
2813
# File 'bigdecimal.c', line 2696

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);
}

.logObject

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.


2827
2828
2829
2830
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849
2850
2851
2852
2853
2854
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884
2885
2886
2887
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936
2937
2938
2939
2940
2941
2942
2943
2944
2945
2946
2947
2948
2949
2950
2951
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
# File 'bigdecimal.c', line 2827

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:
	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%"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)

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

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

57
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
# 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.14142135623730950488016887242096975E1"

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

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