1require 'bigdecimal'
2
3#
4#--
5# Contents:
6#   sqrt(x, prec)
7#   sin (x, prec)
8#   cos (x, prec)
9#   atan(x, prec)  Note: |x|<1, x=0.9999 may not converge.
10#   PI  (prec)
11#   E   (prec) == exp(1.0,prec)
12#
13# where:
14#   x    ... BigDecimal number to be computed.
15#            |x| must be small enough to get convergence.
16#   prec ... Number of digits to be obtained.
17#++
18#
19# Provides mathematical functions.
20#
21# Example:
22#
23#   require "bigdecimal"
24#   require "bigdecimal/math"
25#
26#   include BigMath
27#
28#   a = BigDecimal((PI(100)/2).to_s)
29#   puts sin(a,100) # -> 0.10000000000000000000......E1
30#
31module BigMath
32  module_function
33
34  # Computes the square root of x to the specified number of digits of
35  # precision.
36  #
37  # BigDecimal.new('2').sqrt(16).to_s
38  #  -> "0.14142135623730950488016887242096975E1"
39  #
40  def sqrt(x,prec)
41    x.sqrt(prec)
42  end
43
44  # Computes the sine of x to the specified number of digits of precision.
45  #
46  # If x is infinite or NaN, returns NaN.
47  def sin(x, prec)
48    raise ArgumentError, "Zero or negative precision for sin" if prec <= 0
49    return BigDecimal("NaN") if x.infinite? || x.nan?
50    n    = prec + BigDecimal.double_fig
51    one  = BigDecimal("1")
52    two  = BigDecimal("2")
53    x = -x if neg = x < 0
54    if x > (twopi = two * BigMath.PI(prec))
55      if x > 30
56        x %= twopi
57      else
58        x -= twopi while x > twopi
59      end
60    end
61    x1   = x
62    x2   = x.mult(x,n)
63    sign = 1
64    y    = x
65    d    = y
66    i    = one
67    z    = one
68    while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
69      m = BigDecimal.double_fig if m < BigDecimal.double_fig
70      sign = -sign
71      x1  = x2.mult(x1,n)
72      i  += two
73      z  *= (i-one) * i
74      d   = sign * x1.div(z,m)
75      y  += d
76    end
77    neg ? -y : y
78  end
79
80  # Computes the cosine of x to the specified number of digits of precision.
81  #
82  # If x is infinite or NaN, returns NaN.
83  def cos(x, prec)
84    raise ArgumentError, "Zero or negative precision for cos" if prec <= 0
85    return BigDecimal("NaN") if x.infinite? || x.nan?
86    n    = prec + BigDecimal.double_fig
87    one  = BigDecimal("1")
88    two  = BigDecimal("2")
89    x = -x if x < 0
90    if x > (twopi = two * BigMath.PI(prec))
91      if x > 30
92        x %= twopi
93      else
94        x -= twopi while x > twopi
95      end
96    end
97    x1 = one
98    x2 = x.mult(x,n)
99    sign = 1
100    y = one
101    d = y
102    i = BigDecimal("0")
103    z = one
104    while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
105      m = BigDecimal.double_fig if m < BigDecimal.double_fig
106      sign = -sign
107      x1  = x2.mult(x1,n)
108      i  += two
109      z  *= (i-one) * i
110      d   = sign * x1.div(z,m)
111      y  += d
112    end
113    y
114  end
115
116  # Computes the arctangent of x to the specified number of digits of precision.
117  #
118  # If x is NaN, returns NaN.
119  def atan(x, prec)
120    raise ArgumentError, "Zero or negative precision for atan" if prec <= 0
121    return BigDecimal("NaN") if x.nan?
122    pi = PI(prec)
123    x = -x if neg = x < 0
124    return pi.div(neg ? -2 : 2, prec) if x.infinite?
125    return pi / (neg ? -4 : 4) if x.round(prec) == 1
126    x = BigDecimal("1").div(x, prec) if inv = x > 1
127    x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
128    n    = prec + BigDecimal.double_fig
129    y = x
130    d = y
131    t = x
132    r = BigDecimal("3")
133    x2 = x.mult(x,n)
134    while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
135      m = BigDecimal.double_fig if m < BigDecimal.double_fig
136      t = -t.mult(x2,n)
137      d = t.div(r,m)
138      y += d
139      r += 2
140    end
141    y *= 2 if dbl
142    y = pi / 2 - y if inv
143    y = -y if neg
144    y
145  end
146
147  # Computes the value of pi to the specified number of digits of precision.
148  def PI(prec)
149    raise ArgumentError, "Zero or negative argument for PI" if prec <= 0
150    n      = prec + BigDecimal.double_fig
151    zero   = BigDecimal("0")
152    one    = BigDecimal("1")
153    two    = BigDecimal("2")
154
155    m25    = BigDecimal("-0.04")
156    m57121 = BigDecimal("-57121")
157
158    pi     = zero
159
160    d = one
161    k = one
162    w = one
163    t = BigDecimal("-80")
164    while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
165      m = BigDecimal.double_fig if m < BigDecimal.double_fig
166      t   = t*m25
167      d   = t.div(k,m)
168      k   = k+two
169      pi  = pi + d
170    end
171
172    d = one
173    k = one
174    w = one
175    t = BigDecimal("956")
176    while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
177      m = BigDecimal.double_fig if m < BigDecimal.double_fig
178      t   = t.div(m57121,n)
179      d   = t.div(k,m)
180      pi  = pi + d
181      k   = k+two
182    end
183    pi
184  end
185
186  # Computes e (the base of natural logarithms) to the specified number of
187  # digits of precision.
188  def E(prec)
189    raise ArgumentError, "Zero or negative precision for E" if prec <= 0
190    n    = prec + BigDecimal.double_fig
191    one  = BigDecimal("1")
192    y  = one
193    d  = y
194    z  = one
195    i  = 0
196    while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
197      m = BigDecimal.double_fig if m < BigDecimal.double_fig
198      i += 1
199      z *= i
200      d  = one.div(z,m)
201      y += d
202    end
203    y
204  end
205end
206