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