gaussiandomains.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645
  1. """Domains of Gaussian type."""
  2. from sympy.core.numbers import I
  3. from sympy.polys.polyerrors import CoercionFailed
  4. from sympy.polys.domains.integerring import ZZ
  5. from sympy.polys.domains.rationalfield import QQ
  6. from sympy.polys.domains.algebraicfield import AlgebraicField
  7. from sympy.polys.domains.domain import Domain
  8. from sympy.polys.domains.domainelement import DomainElement
  9. from sympy.polys.domains.field import Field
  10. from sympy.polys.domains.ring import Ring
  11. class GaussianElement(DomainElement):
  12. """Base class for elements of Gaussian type domains."""
  13. base = None # type: Domain
  14. _parent = None # type: Domain
  15. __slots__ = ('x', 'y')
  16. def __new__(cls, x, y=0):
  17. conv = cls.base.convert
  18. return cls.new(conv(x), conv(y))
  19. @classmethod
  20. def new(cls, x, y):
  21. """Create a new GaussianElement of the same domain."""
  22. obj = super().__new__(cls)
  23. obj.x = x
  24. obj.y = y
  25. return obj
  26. def parent(self):
  27. """The domain that this is an element of (ZZ_I or QQ_I)"""
  28. return self._parent
  29. def __hash__(self):
  30. return hash((self.x, self.y))
  31. def __eq__(self, other):
  32. if isinstance(other, self.__class__):
  33. return self.x == other.x and self.y == other.y
  34. else:
  35. return NotImplemented
  36. def __lt__(self, other):
  37. if not isinstance(other, GaussianElement):
  38. return NotImplemented
  39. return [self.y, self.x] < [other.y, other.x]
  40. def __pos__(self):
  41. return self
  42. def __neg__(self):
  43. return self.new(-self.x, -self.y)
  44. def __repr__(self):
  45. return "%s(%s, %s)" % (self._parent.rep, self.x, self.y)
  46. def __str__(self):
  47. return str(self._parent.to_sympy(self))
  48. @classmethod
  49. def _get_xy(cls, other):
  50. if not isinstance(other, cls):
  51. try:
  52. other = cls._parent.convert(other)
  53. except CoercionFailed:
  54. return None, None
  55. return other.x, other.y
  56. def __add__(self, other):
  57. x, y = self._get_xy(other)
  58. if x is not None:
  59. return self.new(self.x + x, self.y + y)
  60. else:
  61. return NotImplemented
  62. __radd__ = __add__
  63. def __sub__(self, other):
  64. x, y = self._get_xy(other)
  65. if x is not None:
  66. return self.new(self.x - x, self.y - y)
  67. else:
  68. return NotImplemented
  69. def __rsub__(self, other):
  70. x, y = self._get_xy(other)
  71. if x is not None:
  72. return self.new(x - self.x, y - self.y)
  73. else:
  74. return NotImplemented
  75. def __mul__(self, other):
  76. x, y = self._get_xy(other)
  77. if x is not None:
  78. return self.new(self.x*x - self.y*y, self.x*y + self.y*x)
  79. else:
  80. return NotImplemented
  81. __rmul__ = __mul__
  82. def __pow__(self, exp):
  83. if exp == 0:
  84. return self.new(1, 0)
  85. if exp < 0:
  86. self, exp = 1/self, -exp
  87. if exp == 1:
  88. return self
  89. pow2 = self
  90. prod = self if exp % 2 else self._parent.one
  91. exp //= 2
  92. while exp:
  93. pow2 *= pow2
  94. if exp % 2:
  95. prod *= pow2
  96. exp //= 2
  97. return prod
  98. def __bool__(self):
  99. return bool(self.x) or bool(self.y)
  100. def quadrant(self):
  101. """Return quadrant index 0-3.
  102. 0 is included in quadrant 0.
  103. """
  104. if self.y > 0:
  105. return 0 if self.x > 0 else 1
  106. elif self.y < 0:
  107. return 2 if self.x < 0 else 3
  108. else:
  109. return 0 if self.x >= 0 else 2
  110. def __rdivmod__(self, other):
  111. try:
  112. other = self._parent.convert(other)
  113. except CoercionFailed:
  114. return NotImplemented
  115. else:
  116. return other.__divmod__(self)
  117. def __rtruediv__(self, other):
  118. try:
  119. other = QQ_I.convert(other)
  120. except CoercionFailed:
  121. return NotImplemented
  122. else:
  123. return other.__truediv__(self)
  124. def __floordiv__(self, other):
  125. qr = self.__divmod__(other)
  126. return qr if qr is NotImplemented else qr[0]
  127. def __rfloordiv__(self, other):
  128. qr = self.__rdivmod__(other)
  129. return qr if qr is NotImplemented else qr[0]
  130. def __mod__(self, other):
  131. qr = self.__divmod__(other)
  132. return qr if qr is NotImplemented else qr[1]
  133. def __rmod__(self, other):
  134. qr = self.__rdivmod__(other)
  135. return qr if qr is NotImplemented else qr[1]
  136. class GaussianInteger(GaussianElement):
  137. """Gaussian integer: domain element for :ref:`ZZ_I`
  138. >>> from sympy import ZZ_I
  139. >>> z = ZZ_I(2, 3)
  140. >>> z
  141. (2 + 3*I)
  142. >>> type(z)
  143. <class 'sympy.polys.domains.gaussiandomains.GaussianInteger'>
  144. """
  145. base = ZZ
  146. def __truediv__(self, other):
  147. """Return a Gaussian rational."""
  148. return QQ_I.convert(self)/other
  149. def __divmod__(self, other):
  150. if not other:
  151. raise ZeroDivisionError('divmod({}, 0)'.format(self))
  152. x, y = self._get_xy(other)
  153. if x is None:
  154. return NotImplemented
  155. # multiply self and other by x - I*y
  156. # self/other == (a + I*b)/c
  157. a, b = self.x*x + self.y*y, -self.x*y + self.y*x
  158. c = x*x + y*y
  159. # find integers qx and qy such that
  160. # |a - qx*c| <= c/2 and |b - qy*c| <= c/2
  161. qx = (2*a + c) // (2*c) # -c <= 2*a - qx*2*c < c
  162. qy = (2*b + c) // (2*c)
  163. q = GaussianInteger(qx, qy)
  164. # |self/other - q| < 1 since
  165. # |a/c - qx|**2 + |b/c - qy|**2 <= 1/4 + 1/4 < 1
  166. return q, self - q*other # |r| < |other|
  167. class GaussianRational(GaussianElement):
  168. """Gaussian rational: domain element for :ref:`QQ_I`
  169. >>> from sympy import QQ_I, QQ
  170. >>> z = QQ_I(QQ(2, 3), QQ(4, 5))
  171. >>> z
  172. (2/3 + 4/5*I)
  173. >>> type(z)
  174. <class 'sympy.polys.domains.gaussiandomains.GaussianRational'>
  175. """
  176. base = QQ
  177. def __truediv__(self, other):
  178. """Return a Gaussian rational."""
  179. if not other:
  180. raise ZeroDivisionError('{} / 0'.format(self))
  181. x, y = self._get_xy(other)
  182. if x is None:
  183. return NotImplemented
  184. c = x*x + y*y
  185. return GaussianRational((self.x*x + self.y*y)/c,
  186. (-self.x*y + self.y*x)/c)
  187. def __divmod__(self, other):
  188. try:
  189. other = self._parent.convert(other)
  190. except CoercionFailed:
  191. return NotImplemented
  192. if not other:
  193. raise ZeroDivisionError('{} % 0'.format(self))
  194. else:
  195. return self/other, QQ_I.zero
  196. class GaussianDomain():
  197. """Base class for Gaussian domains."""
  198. dom = None # type: Domain
  199. is_Numerical = True
  200. is_Exact = True
  201. has_assoc_Ring = True
  202. has_assoc_Field = True
  203. def to_sympy(self, a):
  204. """Convert ``a`` to a SymPy object. """
  205. conv = self.dom.to_sympy
  206. return conv(a.x) + I*conv(a.y)
  207. def from_sympy(self, a):
  208. """Convert a SymPy object to ``self.dtype``."""
  209. r, b = a.as_coeff_Add()
  210. x = self.dom.from_sympy(r) # may raise CoercionFailed
  211. if not b:
  212. return self.new(x, 0)
  213. r, b = b.as_coeff_Mul()
  214. y = self.dom.from_sympy(r)
  215. if b is I:
  216. return self.new(x, y)
  217. else:
  218. raise CoercionFailed("{} is not Gaussian".format(a))
  219. def inject(self, *gens):
  220. """Inject generators into this domain. """
  221. return self.poly_ring(*gens)
  222. def canonical_unit(self, d):
  223. unit = self.units[-d.quadrant()] # - for inverse power
  224. return unit
  225. def is_negative(self, element):
  226. """Returns ``False`` for any ``GaussianElement``. """
  227. return False
  228. def is_positive(self, element):
  229. """Returns ``False`` for any ``GaussianElement``. """
  230. return False
  231. def is_nonnegative(self, element):
  232. """Returns ``False`` for any ``GaussianElement``. """
  233. return False
  234. def is_nonpositive(self, element):
  235. """Returns ``False`` for any ``GaussianElement``. """
  236. return False
  237. def from_ZZ_gmpy(K1, a, K0):
  238. """Convert a GMPY mpz to ``self.dtype``."""
  239. return K1(a)
  240. def from_ZZ(K1, a, K0):
  241. """Convert a ZZ_python element to ``self.dtype``."""
  242. return K1(a)
  243. def from_ZZ_python(K1, a, K0):
  244. """Convert a ZZ_python element to ``self.dtype``."""
  245. return K1(a)
  246. def from_QQ(K1, a, K0):
  247. """Convert a GMPY mpq to ``self.dtype``."""
  248. return K1(a)
  249. def from_QQ_gmpy(K1, a, K0):
  250. """Convert a GMPY mpq to ``self.dtype``."""
  251. return K1(a)
  252. def from_QQ_python(K1, a, K0):
  253. """Convert a QQ_python element to ``self.dtype``."""
  254. return K1(a)
  255. def from_AlgebraicField(K1, a, K0):
  256. """Convert an element from ZZ<I> or QQ<I> to ``self.dtype``."""
  257. if K0.ext.args[0] == I:
  258. return K1.from_sympy(K0.to_sympy(a))
  259. class GaussianIntegerRing(GaussianDomain, Ring):
  260. r"""Ring of Gaussian integers ``ZZ_I``
  261. The :ref:`ZZ_I` domain represents the `Gaussian integers`_ `\mathbb{Z}[i]`
  262. as a :py:class:`~.Domain` in the domain system (see
  263. :ref:`polys-domainsintro`).
  264. By default a :py:class:`~.Poly` created from an expression with
  265. coefficients that are combinations of integers and ``I`` (`\sqrt{-1}`)
  266. will have the domain :ref:`ZZ_I`.
  267. >>> from sympy import Poly, Symbol, I
  268. >>> x = Symbol('x')
  269. >>> p = Poly(x**2 + I)
  270. >>> p
  271. Poly(x**2 + I, x, domain='ZZ_I')
  272. >>> p.domain
  273. ZZ_I
  274. The :ref:`ZZ_I` domain can be used to factorise polynomials that are
  275. reducible over the Gaussian integers.
  276. >>> from sympy import factor
  277. >>> factor(x**2 + 1)
  278. x**2 + 1
  279. >>> factor(x**2 + 1, domain='ZZ_I')
  280. (x - I)*(x + I)
  281. The corresponding `field of fractions`_ is the domain of the Gaussian
  282. rationals :ref:`QQ_I`. Conversely :ref:`ZZ_I` is the `ring of integers`_
  283. of :ref:`QQ_I`.
  284. >>> from sympy import ZZ_I, QQ_I
  285. >>> ZZ_I.get_field()
  286. QQ_I
  287. >>> QQ_I.get_ring()
  288. ZZ_I
  289. When using the domain directly :ref:`ZZ_I` can be used as a constructor.
  290. >>> ZZ_I(3, 4)
  291. (3 + 4*I)
  292. >>> ZZ_I(5)
  293. (5 + 0*I)
  294. The domain elements of :ref:`ZZ_I` are instances of
  295. :py:class:`~.GaussianInteger` which support the rings operations
  296. ``+,-,*,**``.
  297. >>> z1 = ZZ_I(5, 1)
  298. >>> z2 = ZZ_I(2, 3)
  299. >>> z1
  300. (5 + 1*I)
  301. >>> z2
  302. (2 + 3*I)
  303. >>> z1 + z2
  304. (7 + 4*I)
  305. >>> z1 * z2
  306. (7 + 17*I)
  307. >>> z1 ** 2
  308. (24 + 10*I)
  309. Both floor (``//``) and modulo (``%``) division work with
  310. :py:class:`~.GaussianInteger` (see the :py:meth:`~.Domain.div` method).
  311. >>> z3, z4 = ZZ_I(5), ZZ_I(1, 3)
  312. >>> z3 // z4 # floor division
  313. (1 + -1*I)
  314. >>> z3 % z4 # modulo division (remainder)
  315. (1 + -2*I)
  316. >>> (z3//z4)*z4 + z3%z4 == z3
  317. True
  318. True division (``/``) in :ref:`ZZ_I` gives an element of :ref:`QQ_I`. The
  319. :py:meth:`~.Domain.exquo` method can be used to divide in :ref:`ZZ_I` when
  320. exact division is possible.
  321. >>> z1 / z2
  322. (1 + -1*I)
  323. >>> ZZ_I.exquo(z1, z2)
  324. (1 + -1*I)
  325. >>> z3 / z4
  326. (1/2 + -3/2*I)
  327. >>> ZZ_I.exquo(z3, z4)
  328. Traceback (most recent call last):
  329. ...
  330. ExactQuotientFailed: (1 + 3*I) does not divide (5 + 0*I) in ZZ_I
  331. The :py:meth:`~.Domain.gcd` method can be used to compute the `gcd`_ of any
  332. two elements.
  333. >>> ZZ_I.gcd(ZZ_I(10), ZZ_I(2))
  334. (2 + 0*I)
  335. >>> ZZ_I.gcd(ZZ_I(5), ZZ_I(2, 1))
  336. (2 + 1*I)
  337. .. _Gaussian integers: https://en.wikipedia.org/wiki/Gaussian_integer
  338. .. _gcd: https://en.wikipedia.org/wiki/Greatest_common_divisor
  339. """
  340. dom = ZZ
  341. dtype = GaussianInteger
  342. zero = dtype(ZZ(0), ZZ(0))
  343. one = dtype(ZZ(1), ZZ(0))
  344. imag_unit = dtype(ZZ(0), ZZ(1))
  345. units = (one, imag_unit, -one, -imag_unit) # powers of i
  346. rep = 'ZZ_I'
  347. is_GaussianRing = True
  348. is_ZZ_I = True
  349. def __init__(self): # override Domain.__init__
  350. """For constructing ZZ_I."""
  351. def get_ring(self):
  352. """Returns a ring associated with ``self``. """
  353. return self
  354. def get_field(self):
  355. """Returns a field associated with ``self``. """
  356. return QQ_I
  357. def normalize(self, d, *args):
  358. """Return first quadrant element associated with ``d``.
  359. Also multiply the other arguments by the same power of i.
  360. """
  361. unit = self.canonical_unit(d)
  362. d *= unit
  363. args = tuple(a*unit for a in args)
  364. return (d,) + args if args else d
  365. def gcd(self, a, b):
  366. """Greatest common divisor of a and b over ZZ_I."""
  367. while b:
  368. a, b = b, a % b
  369. return self.normalize(a)
  370. def lcm(self, a, b):
  371. """Least common multiple of a and b over ZZ_I."""
  372. return (a * b) // self.gcd(a, b)
  373. def from_GaussianIntegerRing(K1, a, K0):
  374. """Convert a ZZ_I element to ZZ_I."""
  375. return a
  376. def from_GaussianRationalField(K1, a, K0):
  377. """Convert a QQ_I element to ZZ_I."""
  378. return K1.new(ZZ.convert(a.x), ZZ.convert(a.y))
  379. ZZ_I = GaussianInteger._parent = GaussianIntegerRing()
  380. class GaussianRationalField(GaussianDomain, Field):
  381. r"""Field of Gaussian rationals ``QQ_I``
  382. The :ref:`QQ_I` domain represents the `Gaussian rationals`_ `\mathbb{Q}(i)`
  383. as a :py:class:`~.Domain` in the domain system (see
  384. :ref:`polys-domainsintro`).
  385. By default a :py:class:`~.Poly` created from an expression with
  386. coefficients that are combinations of rationals and ``I`` (`\sqrt{-1}`)
  387. will have the domain :ref:`QQ_I`.
  388. >>> from sympy import Poly, Symbol, I
  389. >>> x = Symbol('x')
  390. >>> p = Poly(x**2 + I/2)
  391. >>> p
  392. Poly(x**2 + I/2, x, domain='QQ_I')
  393. >>> p.domain
  394. QQ_I
  395. The polys option ``gaussian=True`` can be used to specify that the domain
  396. should be :ref:`QQ_I` even if the coefficients do not contain ``I`` or are
  397. all integers.
  398. >>> Poly(x**2)
  399. Poly(x**2, x, domain='ZZ')
  400. >>> Poly(x**2 + I)
  401. Poly(x**2 + I, x, domain='ZZ_I')
  402. >>> Poly(x**2/2)
  403. Poly(1/2*x**2, x, domain='QQ')
  404. >>> Poly(x**2, gaussian=True)
  405. Poly(x**2, x, domain='QQ_I')
  406. >>> Poly(x**2 + I, gaussian=True)
  407. Poly(x**2 + I, x, domain='QQ_I')
  408. >>> Poly(x**2/2, gaussian=True)
  409. Poly(1/2*x**2, x, domain='QQ_I')
  410. The :ref:`QQ_I` domain can be used to factorise polynomials that are
  411. reducible over the Gaussian rationals.
  412. >>> from sympy import factor, QQ_I
  413. >>> factor(x**2/4 + 1)
  414. (x**2 + 4)/4
  415. >>> factor(x**2/4 + 1, domain='QQ_I')
  416. (x - 2*I)*(x + 2*I)/4
  417. >>> factor(x**2/4 + 1, domain=QQ_I)
  418. (x - 2*I)*(x + 2*I)/4
  419. It is also possible to specify the :ref:`QQ_I` domain explicitly with
  420. polys functions like :py:func:`~.apart`.
  421. >>> from sympy import apart
  422. >>> apart(1/(1 + x**2))
  423. 1/(x**2 + 1)
  424. >>> apart(1/(1 + x**2), domain=QQ_I)
  425. I/(2*(x + I)) - I/(2*(x - I))
  426. The corresponding `ring of integers`_ is the domain of the Gaussian
  427. integers :ref:`ZZ_I`. Conversely :ref:`QQ_I` is the `field of fractions`_
  428. of :ref:`ZZ_I`.
  429. >>> from sympy import ZZ_I, QQ_I, QQ
  430. >>> ZZ_I.get_field()
  431. QQ_I
  432. >>> QQ_I.get_ring()
  433. ZZ_I
  434. When using the domain directly :ref:`QQ_I` can be used as a constructor.
  435. >>> QQ_I(3, 4)
  436. (3 + 4*I)
  437. >>> QQ_I(5)
  438. (5 + 0*I)
  439. >>> QQ_I(QQ(2, 3), QQ(4, 5))
  440. (2/3 + 4/5*I)
  441. The domain elements of :ref:`QQ_I` are instances of
  442. :py:class:`~.GaussianRational` which support the field operations
  443. ``+,-,*,**,/``.
  444. >>> z1 = QQ_I(5, 1)
  445. >>> z2 = QQ_I(2, QQ(1, 2))
  446. >>> z1
  447. (5 + 1*I)
  448. >>> z2
  449. (2 + 1/2*I)
  450. >>> z1 + z2
  451. (7 + 3/2*I)
  452. >>> z1 * z2
  453. (19/2 + 9/2*I)
  454. >>> z2 ** 2
  455. (15/4 + 2*I)
  456. True division (``/``) in :ref:`QQ_I` gives an element of :ref:`QQ_I` and
  457. is always exact.
  458. >>> z1 / z2
  459. (42/17 + -2/17*I)
  460. >>> QQ_I.exquo(z1, z2)
  461. (42/17 + -2/17*I)
  462. >>> z1 == (z1/z2)*z2
  463. True
  464. Both floor (``//``) and modulo (``%``) division can be used with
  465. :py:class:`~.GaussianRational` (see :py:meth:`~.Domain.div`)
  466. but division is always exact so there is no remainder.
  467. >>> z1 // z2
  468. (42/17 + -2/17*I)
  469. >>> z1 % z2
  470. (0 + 0*I)
  471. >>> QQ_I.div(z1, z2)
  472. ((42/17 + -2/17*I), (0 + 0*I))
  473. >>> (z1//z2)*z2 + z1%z2 == z1
  474. True
  475. .. _Gaussian rationals: https://en.wikipedia.org/wiki/Gaussian_rational
  476. """
  477. dom = QQ
  478. dtype = GaussianRational
  479. zero = dtype(QQ(0), QQ(0))
  480. one = dtype(QQ(1), QQ(0))
  481. imag_unit = dtype(QQ(0), QQ(1))
  482. units = (one, imag_unit, -one, -imag_unit) # powers of i
  483. rep = 'QQ_I'
  484. is_GaussianField = True
  485. is_QQ_I = True
  486. def __init__(self): # override Domain.__init__
  487. """For constructing QQ_I."""
  488. def get_ring(self):
  489. """Returns a ring associated with ``self``. """
  490. return ZZ_I
  491. def get_field(self):
  492. """Returns a field associated with ``self``. """
  493. return self
  494. def as_AlgebraicField(self):
  495. """Get equivalent domain as an ``AlgebraicField``. """
  496. return AlgebraicField(self.dom, I)
  497. def numer(self, a):
  498. """Get the numerator of ``a``."""
  499. ZZ_I = self.get_ring()
  500. return ZZ_I.convert(a * self.denom(a))
  501. def denom(self, a):
  502. """Get the denominator of ``a``."""
  503. ZZ = self.dom.get_ring()
  504. QQ = self.dom
  505. ZZ_I = self.get_ring()
  506. denom_ZZ = ZZ.lcm(QQ.denom(a.x), QQ.denom(a.y))
  507. return ZZ_I(denom_ZZ, ZZ.zero)
  508. def from_GaussianIntegerRing(K1, a, K0):
  509. """Convert a ZZ_I element to QQ_I."""
  510. return K1.new(a.x, a.y)
  511. def from_GaussianRationalField(K1, a, K0):
  512. """Convert a QQ_I element to QQ_I."""
  513. return a
  514. QQ_I = GaussianRational._parent = GaussianRationalField()