sympycore/arithmetic/numbers.py

   1 """
   2 This module implements efficient "low-level" number types for purely
   3 numerical tasks.
   4
   5 We use the following types:
   6
   7 int/long -- (Python built-ins) exact integers
   8 mpf -- arbitrary-precision floats
   9 mpc -- arbitrary-precision complex floats
  10 mpq -- nonintegral fractions
  11 mpqc -- nonreal complex rationals
  12
  13 In the interest of speed, there are some quirks:
  14
  15   * mpq derives from tuple. mpq.__init__ does *not* validate its input.
  16     mpqs *should only* be created from fully reduced (p, q) pairs,
  17     and p/q should *not* be integer-valued.
  18
  19   * To safely create a fraction from two integers, use the function
  20     normalized_fraction()
  21
  22   * Arithmetic operations on mpqs automatically normalize back to Python
  23     ints when the results are integer-valued. Likewise, mpqc with zero
  24     real part normalize to their real parts.
  25
  26   * Note that ``mpq(2)/mpq(3)`` does *not* work as expected; it does the
  27     same thing as ``2/3`` in Python. To perform safe division, use the
  28     provided ``div`` function.
  29
  30   * mpqs can be compared with Python ints, but cannot be (correctly)
  31     compared to Python floats. (Use the mpf class instead.)
  32
  33 Powers, except when the exponent is a positive integer, should be
  34 computed with the ``try_power()`` function which detects when the
  35 result is not exact.
  36 """
  37 #
  38 # Author: Fredrik Johansson
  39 # Created: January 2008
  40
  41 import math
  42
  43 from . import mpmath
  44
  45 __docformat__ = "restructuredtext"
  46
  47 from ..core import init_module
  48
  49 init_module.import_heads()
  50
  51 @init_module
  52 def _init(module):
  53     from ..basealgebra.verbatim import Verbatim
  54     module.Verbatim = Verbatim
  55
  56 from ..utils import str_SUM, str_PRODUCT, str_POWER, str_APPLY, str_SYMBOL, str_NUMBER, NUMBER, SYMBOL
  57
  58 inttypes = (int, long)
  59
  60 from .mpmath import mpf, mpc, mp
  61 from .mpmath.libmp import from_rational, round_nearest
  62
  63
  64
  65 def mpf_to_str_data(self, sort=True):
  66     if self < 0:
  67         return str_SUM, str(self)
  68     return str_NUMBER, str(self)
  69
  70 def mpc_to_str_data(self, sort=True):
  71     return str_NUMBER, str(self)
  72
  73 mpf.to_str_data = mpf_to_str_data
  74 mpc.to_str_data = mpc_to_str_data
  75
  76 rounding = round_nearest
  77
  78 def getdps():
  79     return mp.dps
  80
  81 def setdps(n):
  82     p = mp.dps
  83     mp.dps = int(n)
  84     return p
  85
  86 def getprec():
  87     return mp.prec
  88
  89 #----------------------------------------------------------------------------
  90 # Fractions
  91 #
  92
  93 def normalized_fraction(p, q=1):
  94     """ Return a normalized fraction.
  95     """
  96     x, y = p, q
  97     while y:
  98         x, y = y, x % y
  99     if x != 1:
 100         p //= x
 101         q //= x
 102     if q == 1:
 103         return p
 104     return mpq((p, q))
 105
 106 class mpq(tuple):
 107     """Represents a fraction."""
 108
 109     # These methods are inherited directly from tuple for speed. This works
 110     # as long as all mpqs are normalized:
 111     # __new__/__init__
 112     # __nonzero__
 113     # __eq__
 114     # __hash__
 115
 116     # These methods are generated and defined in methods.py:
 117     # __add__, __sub__, __rsub__, __mul__, __div__, __rdiv__, __pow__
 118     # __lt__, __le__, __gt__, __ge__
 119
 120     __slots__ = []
 121
 122     @property
 123     def _mpf_(self):
 124         p, q = self
 125         return from_rational(p, q, getprec(), rounding)
 126
 127     def as_verbatim(self):
 128         p, q = self
 129         if p<0:
 130             return -(Verbatim(NUMBER, -p) / Verbatim(NUMBER, q))
 131         return Verbatim(NUMBER, p) / Verbatim(NUMBER, q)
 132
 133     def to_str_data(self,sort=True):
 134         if self[0]<0:
 135             return str_SUM, str(self)
 136         return str_PRODUCT, str(self)
 137
 138     def __str__(self):
 139         return "%i/%i" % self
 140
 141     def __repr__(self):
 142         p, q = self
 143         return "%s((%s, %s))" % (type(self).__name__, p, q)
 144
 145     def __float__(self):
 146         p, q = self
 147         return float(p) / q
 148
 149     def __int__(self):
 150         p, q = self
 151         return p // q
 152
 153     def __neg__(self):
 154         p, q = self
 155         return mpq((-p, q))
 156
 157     def __pos__(self):
 158         return self
 159
 160     def __abs__(self):
 161         p, q = self
 162         if p < 0:
 163             return mpq((-p, q))
 164         return self
 165
 166     def __floordiv__(a, b):
 167         return int(a / b)
 168
 169     def __rfloordiv__(a, b):
 170         return int(b / a)
 171
 172     def __mod__(a, b):
 173         return a - (a//b)*b
 174
 175     def __rmod__(a, b):
 176         return b - (b//a)*a
 177
 178     def __divmod__(a, b):
 179         return (a-a%b)/b, a%b
 180
 181     def __rdivmod__(a, b):
 182         return (b-b%a)/a, b%a
 183
 184     def __rpow__(a, b):
 185         z, sym = try_power(b, a)
 186         if not sym:
 187             return z
 188         return NotImplemented
 189
 190
 191 #----------------------------------------------------------------------------
 192 # Complex numbers
 193 #
 194
 195 def innerstr(x):
 196     if isinstance(x, mpq):
 197         return "%s/%s" % x
 198     return str(x)
 199
 200
 201 class mpqc(object):
 202     """Represents an exact complex number.
 203
 204     The integer power of a complex number is computed using
 205     `exponentiation by squaring`__ method.
 206
 207     __ http://en.wikipedia.org/wiki/Exponentiation_by_squaring
 208     """
 209
 210     __slots__ = ['real', 'imag']
 211
 212     def __new__(cls, real, imag=0):
 213         if not imag:
 214             return real
 215         if isinstance(real, (float, mpf)) or isinstance(imag, (float, mpf)):
 216             return mpc(real, imag)
 217         self = object.__new__(mpqc)
 218         self.real = real
 219         self.imag = imag
 220         return self
 221
 222     def __repr__(self):
 223         return "%s(%r, %r)" % (self.__class__.__name__, self.real, self.imag)
 224
 225     def __hash__(self):
 226         return hash((self.real, self.imag))
 227
 228     def __mpcval__(self):
 229         return mpc(self.real, self.imag)
 230
 231     def __getstate__(self):
 232         return (self.real, self.imag)
 233
 234     def __setstate__(self, (real, imag)):
 235         self.real = real
 236         self.imag = imag
 237
 238     def to_str_data(self,sort=True):
 239         re, im = self.real, self.imag
 240         if re==0:
 241             if im == 1: return str_SYMBOL,"I"
 242             if im == -1: return str_SUM, "-I"
 243             return str_PRODUCT, str(self.imag) + "*I"
 244         restr = innerstr(self.real)
 245         if im == 1: return str_SUM, "%s + I" % restr
 246         if im == -1: return str_SUM, "%s - I" % restr
 247         if im > 0: return str_SUM, "%s + %s*I" % (restr, innerstr(self.imag))
 248         if im < 0: return str_SUM, "%s - %s*I" % (restr, innerstr(-self.imag))
 249         raise NotImplementedError(`self`)
 250
 251     def __str__(self):
 252         return self.to_str_data()[1]
 253
 254     def as_verbatim(self):
 255         re, im = self.real, self.imag
 256         if re < 0:
 257             r = -Verbatim(NUMBER, -self.real)
 258         else:
 259             r = Verbatim(NUMBER, self.real)
 260         if im < 0:
 261             i = -Verbatim(NUMBER, -self.imag)
 262             ni = Verbatim(NUMBER, -self.imag)
 263         else:
 264             i = Verbatim(NUMBER, self.imag)
 265         I = Verbatim(NUMBER, mpqc(0,1))
 266         if re==0:
 267             if im == 1: return I
 268             if im == -1: return -I
 269             return i * I
 270         if im == 1: return r + I
 271         if im == -1: return r - I
 272         if im < 0:
 273             return r - ni * I
 274         else:
 275             return r + i * I
 276
 277     def as_numer_denom(self):
 278         re, im = self.real, self.imag
 279         if type(re) is mpq:
 280             r1, r2 = re
 281         else:
 282             r1, r2 = re, 1
 283         if type(im) is mpq:
 284             i1, i2 = im
 285         else:
 286             i1, i2 = im, 1
 287         r = r1*i2
 288         i = i1*r2
 289         d = r2*i2
 290         c = gcd(r,i,d)
 291         return mpqc(r//c,i//c), d//c
 292
 293     def __eq__(self, other):
 294         if hasattr(other, "imag"):
 295             return self.real == other.real and self.imag == other.imag
 296         return False
 297
 298     def __pos__(self): return self
 299     def __neg__(self): return mpqc(-self.real, -self.imag)
 300
 301     def __abs__(self):
 302         re, im = self.real, self.imag
 303         if not re:
 304             return abs(im)
 305         m2 = re**2 + im**2
 306         if isinstance(m2, inttypes):
 307             m, flag = int_root(m2,2)
 308             if flag:
 309                 return m
 310         if isinstance(m2, mpq):
 311             p,fp = int_root(m2[0],2)
 312             q,fq = int_root(m2[1],2)
 313             if fp and fq:
 314                 return mpq((p,q))
 315         raise NotImplementedError('abs(%r)' % (self))
 316
 317 #----------------------------------------------------------------------------
 318 # Interface functions
 319 #
 320
 321 inttypes = (int, long)
 322 rationaltypes = (mpq,)
 323 numbertypes = (int, long, float, complex, mpq, mpqc, mpf, mpc)
 324 realtypes = (int, long, float, mpq, mpf)
 325 complextypes = (complex, mpqc, mpc)
 326
 327 inttypes_set = frozenset(inttypes)
 328 rationaltypes_set = frozenset(rationaltypes)
 329 realtypes_set = frozenset(realtypes)
 330 complextypes_set = frozenset(complextypes)
 331 numbertypes_set = frozenset(numbertypes)
 332
 333 def number_div(Algebra, a, b):
 334     if type(b) in inttypes_set:
 335         if not b:
 336             if a:
 337                 return Infinity(Infinity(0))
 338             return Infinity(0)
 339         if b == 1:
 340             return a
 341         if type(a) in inttypes_set:
 342             return normalized_fraction(a, b)
 343     return a / b
 344
 345 def div(a, b):
 346     """Safely compute a/b.
 347
 348     If a or b is an integer, this function makes sure to convert it to
 349     a rational.
 350     """
 351     if type(b) in inttypes_set:
 352         if not b:
 353             return Infinity(a)
 354             raise ZeroDivisionError('%r / %r' % (a, b))
 355         if b == 1:
 356             return a
 357         if type(a) in inttypes_set:
 358             return normalized_fraction(a, b)
 359     return a / b
 360
 361 def int_root(y, n):
 362     """ Return a pair ``(floor(y**(1/n)), x**n == y)``.
 363
 364     Given integers y and n, return a tuple containing ``x =
 365     floor(y**(1/n))`` and a boolean indicating whether ``x**n == y``
 366     exactly.
 367     """
 368     if y < 0: raise ValueError, "y must not be negative"
 369     if n < 1: raise ValueError, "n must be positive"
 370     if y in (0, 1): return y, True
 371     if n == 1: return y, True
 372     if n > y: return 1, False
 373     # Get initial estimate for Newton's method. Care must be taken to
 374     # avoid overflow
 375     try:
 376         guess = int(y ** (1./n)+0.5)
 377     except OverflowError:
 378         try:
 379             guess = int(math.exp(math.log(y)/n)+0.5)
 380         except OverflowError:
 381             guess = 1 << int(math.log(y, 2)/n)
 382     # Newton iteration
 383     xprev, x = -1, guess
 384     while 1:
 385         t = x**(n-1)
 386         xprev, x = x, x - (t*x-y)//(n*t)
 387         if abs(x - xprev) < 1:
 388             break
 389     # Compensate
 390     t = x**n
 391     while t > y:
 392         x -= 1
 393         t = x**n
 394     return x, t == y
 395
 396 def try_power(x, y):
 397     """\
 398     Attempt to compute ``x**y`` where ``x`` and ``y`` must be of the
 399     types int, long, mpq, mpf, or mpqc. The function
 400     returns::
 401
 402       z, symbolic
 403
 404     where ``z`` is a number (i.e. a complex rational) and ``symbolic``
 405     is a list of ``(b, e)`` pairs representing symbolic factors
 406     ``b**e``.
 407
 408     Examples::
 409
 410       try_power(3, 2) --> (9, [])
 411       try_power(2, 1/2) --> (1, [(2, 1/2)])
 412       try_power(45, 1/2) --> (3, [(5, 1/2)])
 413     """
 414     if not y or x == 1:
 415         # (anything)**0 -> 1
 416         # 1**(anything) -> 1
 417         return 1, []
 418     if isinstance(x, Infinity) or isinstance(y, Infinity):
 419         return x**y, []
 420     if isinstance(y, inttypes):
 421         if y >= 0:
 422             return x**y, []
 423         elif not x:
 424             return Infinity.get_zoo(), []
 425         elif isinstance(x, inttypes):
 426             return mpq((1, x**(-y))), []
 427         return x**y, []
 428     elif isinstance(x, inttypes) and isinstance(y, mpq):
 429         if x < 0:
 430             if x==-1:
 431                 p, q = y
 432                 if q==2:
 433                     return mpqc(0, 1)**p, []
 434                 return 1, [(x,y)]
 435             else:
 436                 z, sym = try_power(-x, y)
 437                 z1, sym1 = try_power(-1, y)
 438                 return z * z1, sym+sym1
 439         else:
 440             p, q = y
 441             r, exact = int_root(x, q)
 442             if exact:
 443                 if r==1 or not p:
 444                     return 1, []
 445                 g = r**p if p>0 else mpq((1, r**(-p)))
 446                 return g, []
 447     elif isinstance(x, mpq) and isinstance(y, mpq):
 448         a, b = x
 449         r, rsym = try_power(a, y)
 450         s, ssym = try_power(b, y)
 451         ssym = [(b, -e) for b, e in ssym]
 452         return (div(r,s), rsym + ssym)
 453     elif isinstance(x, mpf) or isinstance(y, mpf):
 454         return x ** y, []
 455     return 1, [(x, y)]
 456
 457 from .number_theory import gcd
 458 from .evalf import evalf
 459 from .infinity import Infinity
 460
 461 from .methods import (\
 462     fraction_add, fraction_sub, fraction_rsub, fraction_mul,
 463     fraction_div, fraction_rdiv, fraction_pow,
 464     fraction_lt, fraction_le, fraction_gt, fraction_ge,
 465     complex_add, complex_sub, complex_rsub, complex_mul,
 466     complex_div, complex_rdiv, complex_pow)
 467
 468 mpq.__add__ = mpq.__radd__ = fraction_add
 469 mpq.__sub__ = fraction_sub
 470 mpq.__rsub__ = fraction_rsub
 471 mpq.__mul__ = mpq.__rmul__ = fraction_mul
 472 mpq.__div__ = fraction_div
 473 mpq.__rdiv__ = fraction_rdiv
 474 mpq.__pow__ = fraction_pow
 475 mpq.__lt__ = fraction_lt
 476 mpq.__le__ = fraction_le
 477 mpq.__gt__ = fraction_gt
 478 mpq.__ge__ = fraction_ge
 479
 480 mpqc.__add__ = mpqc.__radd__ = complex_add
 481 mpqc.__sub__ = complex_sub
 482 mpqc.__rsub__ = complex_rsub
 483 mpqc.__mul__ = mpqc.__rmul__ = complex_mul
 484 mpqc.__div__ = complex_div
 485 mpqc.__rdiv__ = complex_rdiv
 486 mpqc.__pow__ = complex_pow