finitefield.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368
  1. """Implementation of :class:`FiniteField` class. """
  2. import operator
  3. from sympy.external.gmpy import GROUND_TYPES
  4. from sympy.utilities.decorator import doctest_depends_on
  5. from sympy.core.numbers import int_valued
  6. from sympy.polys.domains.field import Field
  7. from sympy.polys.domains.modularinteger import ModularIntegerFactory
  8. from sympy.polys.domains.simpledomain import SimpleDomain
  9. from sympy.polys.galoistools import gf_zassenhaus, gf_irred_p_rabin
  10. from sympy.polys.polyerrors import CoercionFailed
  11. from sympy.utilities import public
  12. from sympy.polys.domains.groundtypes import SymPyInteger
  13. if GROUND_TYPES == 'flint':
  14. __doctest_skip__ = ['FiniteField']
  15. if GROUND_TYPES == 'flint':
  16. import flint
  17. # Don't use python-flint < 0.5.0 because nmod was missing some features in
  18. # previous versions of python-flint and fmpz_mod was not yet added.
  19. _major, _minor, *_ = flint.__version__.split('.')
  20. if (int(_major), int(_minor)) < (0, 5):
  21. flint = None
  22. else:
  23. flint = None
  24. def _modular_int_factory_nmod(mod):
  25. # nmod only recognises int
  26. index = operator.index
  27. mod = index(mod)
  28. nmod = flint.nmod
  29. nmod_poly = flint.nmod_poly
  30. # flint's nmod is only for moduli up to 2^64-1 (on a 64-bit machine)
  31. try:
  32. nmod(0, mod)
  33. except OverflowError:
  34. return None, None
  35. def ctx(x):
  36. try:
  37. return nmod(x, mod)
  38. except TypeError:
  39. return nmod(index(x), mod)
  40. def poly_ctx(cs):
  41. return nmod_poly(cs, mod)
  42. return ctx, poly_ctx
  43. def _modular_int_factory_fmpz_mod(mod):
  44. index = operator.index
  45. fctx = flint.fmpz_mod_ctx(mod)
  46. fctx_poly = flint.fmpz_mod_poly_ctx(mod)
  47. fmpz_mod_poly = flint.fmpz_mod_poly
  48. def ctx(x):
  49. try:
  50. return fctx(x)
  51. except TypeError:
  52. # x might be Integer
  53. return fctx(index(x))
  54. def poly_ctx(cs):
  55. return fmpz_mod_poly(cs, fctx_poly)
  56. return ctx, poly_ctx
  57. def _modular_int_factory(mod, dom, symmetric, self):
  58. # Convert the modulus to ZZ
  59. try:
  60. mod = dom.convert(mod)
  61. except CoercionFailed:
  62. raise ValueError('modulus must be an integer, got %s' % mod)
  63. ctx, poly_ctx, is_flint = None, None, False
  64. # Don't use flint if the modulus is not prime as it often crashes.
  65. if flint is not None and mod.is_prime():
  66. is_flint = True
  67. # Try to use flint's nmod first
  68. ctx, poly_ctx = _modular_int_factory_nmod(mod)
  69. if ctx is None:
  70. # Use fmpz_mod for larger moduli
  71. ctx, poly_ctx = _modular_int_factory_fmpz_mod(mod)
  72. if ctx is None:
  73. # Use the Python implementation if flint is not available or the
  74. # modulus is not prime.
  75. ctx = ModularIntegerFactory(mod, dom, symmetric, self)
  76. poly_ctx = None # not used
  77. return ctx, poly_ctx, is_flint
  78. @public
  79. @doctest_depends_on(modules=['python', 'gmpy'])
  80. class FiniteField(Field, SimpleDomain):
  81. r"""Finite field of prime order :ref:`GF(p)`
  82. A :ref:`GF(p)` domain represents a `finite field`_ `\mathbb{F}_p` of prime
  83. order as :py:class:`~.Domain` in the domain system (see
  84. :ref:`polys-domainsintro`).
  85. A :py:class:`~.Poly` created from an expression with integer
  86. coefficients will have the domain :ref:`ZZ`. However, if the ``modulus=p``
  87. option is given then the domain will be a finite field instead.
  88. >>> from sympy import Poly, Symbol
  89. >>> x = Symbol('x')
  90. >>> p = Poly(x**2 + 1)
  91. >>> p
  92. Poly(x**2 + 1, x, domain='ZZ')
  93. >>> p.domain
  94. ZZ
  95. >>> p2 = Poly(x**2 + 1, modulus=2)
  96. >>> p2
  97. Poly(x**2 + 1, x, modulus=2)
  98. >>> p2.domain
  99. GF(2)
  100. It is possible to factorise a polynomial over :ref:`GF(p)` using the
  101. modulus argument to :py:func:`~.factor` or by specifying the domain
  102. explicitly. The domain can also be given as a string.
  103. >>> from sympy import factor, GF
  104. >>> factor(x**2 + 1)
  105. x**2 + 1
  106. >>> factor(x**2 + 1, modulus=2)
  107. (x + 1)**2
  108. >>> factor(x**2 + 1, domain=GF(2))
  109. (x + 1)**2
  110. >>> factor(x**2 + 1, domain='GF(2)')
  111. (x + 1)**2
  112. It is also possible to use :ref:`GF(p)` with the :py:func:`~.cancel`
  113. and :py:func:`~.gcd` functions.
  114. >>> from sympy import cancel, gcd
  115. >>> cancel((x**2 + 1)/(x + 1))
  116. (x**2 + 1)/(x + 1)
  117. >>> cancel((x**2 + 1)/(x + 1), domain=GF(2))
  118. x + 1
  119. >>> gcd(x**2 + 1, x + 1)
  120. 1
  121. >>> gcd(x**2 + 1, x + 1, domain=GF(2))
  122. x + 1
  123. When using the domain directly :ref:`GF(p)` can be used as a constructor
  124. to create instances which then support the operations ``+,-,*,**,/``
  125. >>> from sympy import GF
  126. >>> K = GF(5)
  127. >>> K
  128. GF(5)
  129. >>> x = K(3)
  130. >>> y = K(2)
  131. >>> x
  132. 3 mod 5
  133. >>> y
  134. 2 mod 5
  135. >>> x * y
  136. 1 mod 5
  137. >>> x / y
  138. 4 mod 5
  139. Notes
  140. =====
  141. It is also possible to create a :ref:`GF(p)` domain of **non-prime**
  142. order but the resulting ring is **not** a field: it is just the ring of
  143. the integers modulo ``n``.
  144. >>> K = GF(9)
  145. >>> z = K(3)
  146. >>> z
  147. 3 mod 9
  148. >>> z**2
  149. 0 mod 9
  150. It would be good to have a proper implementation of prime power fields
  151. (``GF(p**n)``) but these are not yet implemented in SymPY.
  152. .. _finite field: https://en.wikipedia.org/wiki/Finite_field
  153. """
  154. rep = 'FF'
  155. alias = 'FF'
  156. is_FiniteField = is_FF = True
  157. is_Numerical = True
  158. has_assoc_Ring = False
  159. has_assoc_Field = True
  160. dom = None
  161. mod = None
  162. def __init__(self, mod, symmetric=True):
  163. from sympy.polys.domains import ZZ
  164. dom = ZZ
  165. if mod <= 0:
  166. raise ValueError('modulus must be a positive integer, got %s' % mod)
  167. ctx, poly_ctx, is_flint = _modular_int_factory(mod, dom, symmetric, self)
  168. self.dtype = ctx
  169. self._poly_ctx = poly_ctx
  170. self._is_flint = is_flint
  171. self.zero = self.dtype(0)
  172. self.one = self.dtype(1)
  173. self.dom = dom
  174. self.mod = mod
  175. self.sym = symmetric
  176. self._tp = type(self.zero)
  177. @property
  178. def tp(self):
  179. return self._tp
  180. @property
  181. def is_Field(self):
  182. is_field = getattr(self, '_is_field', None)
  183. if is_field is None:
  184. from sympy.ntheory.primetest import isprime
  185. self._is_field = is_field = isprime(self.mod)
  186. return is_field
  187. def __str__(self):
  188. return 'GF(%s)' % self.mod
  189. def __hash__(self):
  190. return hash((self.__class__.__name__, self.dtype, self.mod, self.dom))
  191. def __eq__(self, other):
  192. """Returns ``True`` if two domains are equivalent. """
  193. return isinstance(other, FiniteField) and \
  194. self.mod == other.mod and self.dom == other.dom
  195. def characteristic(self):
  196. """Return the characteristic of this domain. """
  197. return self.mod
  198. def get_field(self):
  199. """Returns a field associated with ``self``. """
  200. return self
  201. def to_sympy(self, a):
  202. """Convert ``a`` to a SymPy object. """
  203. return SymPyInteger(self.to_int(a))
  204. def from_sympy(self, a):
  205. """Convert SymPy's Integer to SymPy's ``Integer``. """
  206. if a.is_Integer:
  207. return self.dtype(self.dom.dtype(int(a)))
  208. elif int_valued(a):
  209. return self.dtype(self.dom.dtype(int(a)))
  210. else:
  211. raise CoercionFailed("expected an integer, got %s" % a)
  212. def to_int(self, a):
  213. """Convert ``val`` to a Python ``int`` object. """
  214. aval = int(a)
  215. if self.sym and aval > self.mod // 2:
  216. aval -= self.mod
  217. return aval
  218. def is_positive(self, a):
  219. """Returns True if ``a`` is positive. """
  220. return bool(a)
  221. def is_nonnegative(self, a):
  222. """Returns True if ``a`` is non-negative. """
  223. return True
  224. def is_negative(self, a):
  225. """Returns True if ``a`` is negative. """
  226. return False
  227. def is_nonpositive(self, a):
  228. """Returns True if ``a`` is non-positive. """
  229. return not a
  230. def from_FF(K1, a, K0=None):
  231. """Convert ``ModularInteger(int)`` to ``dtype``. """
  232. return K1.dtype(K1.dom.from_ZZ(int(a), K0.dom))
  233. def from_FF_python(K1, a, K0=None):
  234. """Convert ``ModularInteger(int)`` to ``dtype``. """
  235. return K1.dtype(K1.dom.from_ZZ_python(int(a), K0.dom))
  236. def from_ZZ(K1, a, K0=None):
  237. """Convert Python's ``int`` to ``dtype``. """
  238. return K1.dtype(K1.dom.from_ZZ_python(a, K0))
  239. def from_ZZ_python(K1, a, K0=None):
  240. """Convert Python's ``int`` to ``dtype``. """
  241. return K1.dtype(K1.dom.from_ZZ_python(a, K0))
  242. def from_QQ(K1, a, K0=None):
  243. """Convert Python's ``Fraction`` to ``dtype``. """
  244. if a.denominator == 1:
  245. return K1.from_ZZ_python(a.numerator)
  246. def from_QQ_python(K1, a, K0=None):
  247. """Convert Python's ``Fraction`` to ``dtype``. """
  248. if a.denominator == 1:
  249. return K1.from_ZZ_python(a.numerator)
  250. def from_FF_gmpy(K1, a, K0=None):
  251. """Convert ``ModularInteger(mpz)`` to ``dtype``. """
  252. return K1.dtype(K1.dom.from_ZZ_gmpy(a.val, K0.dom))
  253. def from_ZZ_gmpy(K1, a, K0=None):
  254. """Convert GMPY's ``mpz`` to ``dtype``. """
  255. return K1.dtype(K1.dom.from_ZZ_gmpy(a, K0))
  256. def from_QQ_gmpy(K1, a, K0=None):
  257. """Convert GMPY's ``mpq`` to ``dtype``. """
  258. if a.denominator == 1:
  259. return K1.from_ZZ_gmpy(a.numerator)
  260. def from_RealField(K1, a, K0):
  261. """Convert mpmath's ``mpf`` to ``dtype``. """
  262. p, q = K0.to_rational(a)
  263. if q == 1:
  264. return K1.dtype(K1.dom.dtype(p))
  265. def is_square(self, a):
  266. """Returns True if ``a`` is a quadratic residue modulo p. """
  267. # a is not a square <=> x**2-a is irreducible
  268. poly = [int(x) for x in [self.one, self.zero, -a]]
  269. return not gf_irred_p_rabin(poly, self.mod, self.dom)
  270. def exsqrt(self, a):
  271. """Square root modulo p of ``a`` if it is a quadratic residue.
  272. Explanation
  273. ===========
  274. Always returns the square root that is no larger than ``p // 2``.
  275. """
  276. # x**2-a is not square-free if a=0 or the field is characteristic 2
  277. if self.mod == 2 or a == 0:
  278. return a
  279. # Otherwise, use square-free factorization routine to factorize x**2-a
  280. poly = [int(x) for x in [self.one, self.zero, -a]]
  281. for factor in gf_zassenhaus(poly, self.mod, self.dom):
  282. if len(factor) == 2 and factor[1] <= self.mod // 2:
  283. return self.dtype(factor[1])
  284. return None
  285. FF = GF = FiniteField