sympy/core/decimal_math.py

   1 """
   2 Decimal math functions.
   3
   4 Heavily copied from http://code.google.com/p/dmath/, generalized pow to non-integer exponents.
   5
   6 Author: Pearu Peterson <pearu.peterson@gmail.com>
   7 Created: February, 2007
   8 """
   9
  10 __all__ = ['acos', 'asin', 'atan', 'atan2', 'ceiling', 'cos', 'cosh', 'degrees',
  11            'e', 'exp', 'floor', 'golden_ratio', 'hypot', 'log', 'log10', 'pi',
  12            'pow', 'radians', 'sign', 'sin', 'sinh', 'sqrt', 'tan', 'tanh',
  13            'cot','coth']
  14
  15 import math, decimal
  16
  17 from decimal import getcontext
  18 from decimal import Decimal as D
  19
  20 def pi():
  21     """Compute Pi to the current precision."""
  22     getcontext().prec += 2
  23     lasts, t, s, n, na, d, da = 0, D(3), 3, 1, 0, 0, 24
  24     while s != lasts:
  25         lasts = s
  26         n, na = n + na, na + 8
  27         d, da = d + da, da + 32
  28         t = (t * n) / d
  29         s += t
  30     getcontext().prec -= 2
  31     return +s
  32
  33 def e():
  34     """Compute the base of the natural logarithm to the current precision."""
  35     return exp(D(1))
  36
  37 def golden_ratio():
  38     """Calculate the golden ratio to the current precision."""
  39     return  +((1 + D(5).sqrt()) / 2)
  40
  41
  42 def exp(x):
  43     """Return e raised to the power of x.
  44     """
  45     getcontext().prec += 2
  46     i, lasts, s, fact, num = 0, 0, 1, 1, 1
  47     while s != lasts:
  48         lasts = s
  49         i += 1
  50         fact *= i
  51         num *= x
  52         s += num / fact
  53     getcontext().prec -= 2
  54     return +s
  55
  56 def cos(x):
  57     """Return the cosine of x as measured in radians.
  58     """
  59     getcontext().prec += 2
  60     i, lasts, s, fact, num, sign = 0, 0, 1, 1, 1, 1
  61     while s != lasts:
  62         lasts = s
  63         i += 2
  64         fact *= i * (i - 1)
  65         num *= x * x
  66         sign *= -1
  67         s += num / fact * sign
  68     getcontext().prec -= 2
  69     return +s
  70
  71 def sin(x):
  72     """Return the sine of x as measured in radians.
  73     """
  74     getcontext().prec += 2
  75     i, lasts, s, fact, num, sign = 1, 0, x, 1, x, 1
  76     while s != lasts:
  77         lasts = s
  78         i += 2
  79         fact *= i * (i - 1)
  80         num *= x * x
  81         sign *= -1
  82         s += num / fact * sign
  83     getcontext().prec -= 2
  84     return +s
  85
  86 def cosh(x):
  87     """Return the hyperbolic cosine of Decimal x."""
  88     if x == 0:
  89         return D(1)
  90
  91     getcontext().prec += 2
  92     i, lasts, s, fact, num = 0, 0, 1, 1, 1
  93     while s != lasts:
  94         lasts = s
  95         i += 2
  96         num *= x * x
  97         fact *= i * (i - 1)
  98         s += num / fact
  99     getcontext().prec -= 2
 100     return +s
 101
 102 def sinh(x):
 103     """Return the hyperbolic sine of Decimal x."""
 104     if x == 0:
 105         return D(0)
 106
 107     getcontext().prec += 2
 108     i, lasts, s, fact, num = 1, 0, x, 1, x
 109     while s != lasts:
 110         lasts = s
 111         i += 2
 112         num *= x * x
 113         fact *= i * (i - 1)
 114         s += num / fact
 115     getcontext().prec -= 2
 116     return +s
 117
 118 def asin(x):
 119     """Return the arc sine (measured in radians) of Decimal x."""
 120     if abs(x) > 1:
 121         raise ValueError("Domain error: asin accepts -1 <= x <= 1")
 122
 123     if x == -1:
 124         return pi() / -2
 125     elif x == 0:
 126         return D(0)
 127     elif x == 1:
 128         return pi() / 2
 129
 130     return atan2(x, D.sqrt(1 - x ** 2))
 131
 132 def acos(x):
 133     """Return the arc cosine (measured in radians) of Decimal x."""
 134     if abs(x) > 1:
 135         raise ValueError("Domain error: acos accepts -1 <= x <= 1")
 136
 137     if x == -1:
 138         return pi()
 139     elif x == 0:
 140         return pi() / 2
 141     elif x == 1:
 142         return D(0)
 143
 144     return pi() / 2 - atan2(x, D.sqrt(1 - x ** 2))
 145
 146 def tan(x):
 147     """Return the tangent of Decimal x (measured in radians)."""
 148     return +(sin(x) / cos(x))
 149
 150 def tanh(x):
 151     """Return the hyperbolic tangent of Decimal x."""
 152     return +(sinh(x) / cosh(x))
 153
 154 def cot(x):
 155     """Return the cotangent of Decimal x (measured in radians)."""
 156     return +(cos(x) / sin(x))
 157
 158 def coth(x):
 159     """Return the hyperbolic cotangent of Decimal x."""
 160     return +(cosh(x) / sinh(x))
 161
 162 def atan(x):
 163     """Return the arc tangent (measured in radians) of Decimal x."""
 164     if x == D('-Inf'):
 165         return pi() / -2
 166     elif x == 0:
 167         return D(0)
 168     elif x == D('Inf'):
 169         return pi() / 2
 170
 171     if x < -1:
 172         c = pi() / -2
 173         x = 1 / x
 174     elif x > 1:
 175         c = pi() / 2
 176         x = 1 / x
 177     else:
 178         c = 0
 179
 180     getcontext().prec += 2
 181     x_squared = x ** 2
 182     y = x_squared / (1 + x_squared)
 183     y_over_x = y / x
 184     i, lasts, s, coeff, num = D(0), 0, y_over_x, 1, y_over_x
 185     while s != lasts:
 186         lasts = s
 187         i += 2
 188         coeff *= i / (i + 1)
 189         num *= y
 190         s += coeff * num
 191     if c:
 192         s = c - s
 193     getcontext().prec -= 2
 194     return +s
 195
 196 def sign(x):
 197     """Return -1 for negative numbers and 1 for positive numbers."""
 198     return 2 * D(x >= 0) - 1
 199
 200 def atan2(y, x):
 201     """Return the arc tangent (measured in radians) of y/x.
 202     Unlike atan(y/x), the signs of both x and y are considered.
 203     """
 204     abs_y = abs(y)
 205     abs_x = abs(x)
 206     y_is_real = abs_y != D('Inf')
 207
 208     if x:
 209         if y_is_real:
 210             a = y and atan(y / x) or D(0)
 211             if x < 0:
 212                 a += sign(y) * pi()
 213             return a
 214         elif abs_y == abs_x:
 215             x = sign(x)
 216             y = sign(y)
 217             return pi() * (D(2) * abs(x) - x) / (D(4) * y)
 218     if y:
 219         return atan(sign(y) * D('Inf'))
 220     elif x < 0:
 221         return sign(y) * pi()
 222     else:
 223         return D(0)
 224
 225 def log(x, base=None):
 226     """log(x[, base]) -> the logarithm of Decimal x to the given Decimal base.
 227     If the base not specified, returns the natural logarithm (base e) of x.
 228     """
 229     if x < 0:
 230         return D('NaN')
 231     elif base == 1:
 232         raise ValueError("Base was 1!")
 233     elif x == base:
 234         return D(1)
 235     elif x == 0:
 236         return D('-Inf')
 237
 238     getcontext().prec += 2
 239
 240     if base is None:
 241         log_base = 1
 242         approx = math.log(x)
 243     else:
 244         log_base = log(base)
 245         approx = math.log(x, base)
 246
 247     lasts, s = 0, D(repr(approx))
 248     while lasts != s:
 249         lasts = s
 250         s = s - 1 + x / exp(s)
 251     s /= log_base
 252     getcontext().prec -= 2
 253     return +s
 254
 255 def log10(x):
 256     """log10(x) -> the base 10 logarithm of Decimal x."""
 257     return log(x, D(10))
 258
 259 def pow(x, y):
 260     """pow(x,y) -> x ** y."""
 261     if isinstance(y,(int,long)):
 262         return D.__pow__(x, y)
 263     if isinstance(y,D) and y==D('0.5'):
 264         return sqrt(x)
 265     return exp(y * log(x)) # this is slow
 266
 267 sqrt = D.sqrt
 268
 269 def degrees(x):
 270     """degrees(x) -> converts Decimal angle x from radians to degrees"""
 271     return +(x * 180 / pi())
 272
 273 def radians(x):
 274     """radians(x) -> converts Decimal angle x from degrees to radians"""
 275     return +(x * pi() / 180)
 276
 277 def ceiling(x):
 278     """Return the smallest integral value >= x."""
 279     return x.to_integral(rounding=decimal.ROUND_CEILING)
 280
 281 def floor(x):
 282     """Return the largest integral value <= x."""
 283     return x.to_integral(rounding=decimal.ROUND_FLOOR)
 284
 285 def hypot(x, y):
 286     """Return the Euclidean distance, sqrt(x*x + y*y)."""
 287     return sqrt(x * x + y * y)
 288
 289 # EOF