This commit is contained in:
lz_db
2025-11-16 12:31:03 +08:00
commit 0fab423a18
1451 changed files with 743213 additions and 0 deletions

View File

@@ -0,0 +1,35 @@
"""
The routines here were removed from numbers.py, power.py,
digits.py and factor_.py so they could be imported into core
without raising circular import errors.
Although the name 'intfunc' was chosen to represent functions that
work with integers, it can also be thought of as containing
internal/core functions that are needed by the classes of the core.
"""
from ..external.gmpy import (gcdext)
def igcdex(a, b):
"""Returns x, y, g such that g = x*a + y*b = gcd(a, b).
Examples
========
>>> from sympy.core.intfunc import igcdex
>>> igcdex(2, 3)
(-1, 1, 1)
>>> igcdex(10, 12)
(-1, 1, 2)
>>> x, y, g = igcdex(100, 2004)
>>> x, y, g
(-20, 1, 4)
>>> x*100 + y*2004
4
"""
if (not a) and (not b):
return (0, 1, 0)
g, x, y = gcdext(int(a), int(b))
return x, y, g

View File

View File

@@ -0,0 +1,345 @@
import os
from ctypes import c_long, sizeof
from functools import reduce
from typing import Tuple as tTuple, Type
from warnings import warn
from .importtools import import_module
from .pythonmpq import PythonMPQ
from .ntheory import (
# bit_scan1 as python_bit_scan1,
# bit_scan0 as python_bit_scan0,
# remove as python_remove,
# factorial as python_factorial,
# sqrt as python_sqrt,
# sqrtrem as python_sqrtrem,
# gcd as python_gcd,
# lcm as python_lcm,
gcdext as python_gcdext,
# is_square as python_is_square,
# invert as python_invert,
# legendre as python_legendre,
# jacobi as python_jacobi,
# kronecker as python_kronecker,
# iroot as python_iroot,
# is_fermat_prp as python_is_fermat_prp,
# is_euler_prp as python_is_euler_prp,
# is_strong_prp as python_is_strong_prp,
# is_fibonacci_prp as python_is_fibonacci_prp,
# is_lucas_prp as python_is_lucas_prp,
# is_selfridge_prp as python_is_selfridge_prp,
# is_strong_lucas_prp as python_is_strong_lucas_prp,
# is_strong_selfridge_prp as python_is_strong_selfridge_prp,
# is_bpsw_prp as python_is_bpsw_prp,
# is_strong_bpsw_prp as python_is_strong_bpsw_prp,
)
__all__ = [
# GROUND_TYPES is either 'gmpy' or 'python' depending on which is used. If
# gmpy is installed then it will be used unless the environment variable
# SYMPY_GROUND_TYPES is set to something other than 'auto', 'gmpy', or
# 'gmpy2'.
'GROUND_TYPES',
# If HAS_GMPY is 0, no supported version of gmpy is available. Otherwise,
# HAS_GMPY will be 2 for gmpy2 if GROUND_TYPES is 'gmpy'. It used to be
# possible for HAS_GMPY to be 1 for gmpy but gmpy is no longer supported.
'HAS_GMPY',
# SYMPY_INTS is a tuple containing the base types for valid integer types.
# This is either (int,) or (int, type(mpz(0))) depending on GROUND_TYPES.
'SYMPY_INTS',
# MPQ is either gmpy.mpq or the Python equivalent from
# sympy.external.pythonmpq
'MPQ',
# MPZ is either gmpy.mpz or int.
'MPZ',
# 'bit_scan1',
# 'bit_scan0',
# 'remove',
# 'factorial',
# 'sqrt',
# 'is_square',
# 'sqrtrem',
# 'gcd',
# 'lcm',
'gcdext',
# 'invert',
# 'legendre',
# 'jacobi',
# 'kronecker',
# 'iroot',
# 'is_fermat_prp',
# 'is_euler_prp',
# 'is_strong_prp',
# 'is_fibonacci_prp',
# 'is_lucas_prp',
# 'is_selfridge_prp',
# 'is_strong_lucas_prp',
# 'is_strong_selfridge_prp',
# 'is_bpsw_prp',
# 'is_strong_bpsw_prp',
]
#
# Tested python-flint version. Future versions might work but we will only use
# them if explicitly requested by SYMPY_GROUND_TYPES=flint.
#
_PYTHON_FLINT_VERSION_NEEDED = "0.6.*"
def _flint_version_okay(flint_version):
flint_ver = flint_version.split('.')[:2]
needed_ver = _PYTHON_FLINT_VERSION_NEEDED.split('.')[:2]
return flint_ver == needed_ver
#
# We will only use gmpy2 >= 2.0.0
#
_GMPY2_MIN_VERSION = '2.0.0'
def _get_flint(sympy_ground_types):
if sympy_ground_types not in ('auto', 'flint'):
return None
try:
import flint
# Earlier versions of python-flint may not have __version__.
from flint import __version__ as _flint_version
except ImportError:
if sympy_ground_types == 'flint':
warn("SYMPY_GROUND_TYPES was set to flint but python-flint is not "
"installed. Falling back to other ground types.")
return None
if _flint_version_okay(_flint_version):
return flint
elif sympy_ground_types == 'auto':
warn(f"python-flint {_flint_version} is installed but only version "
f"{_PYTHON_FLINT_VERSION_NEEDED} will be used by default. "
f"Falling back to other ground types. Use "
f"SYMPY_GROUND_TYPES=flint to force the use of python-flint.")
return None
else:
warn(f"Using python-flint {_flint_version} because SYMPY_GROUND_TYPES "
f"is set to flint but this version of SymPy has only been tested "
f"with python-flint {_PYTHON_FLINT_VERSION_NEEDED}.")
return flint
def _get_gmpy2(sympy_ground_types):
if sympy_ground_types not in ('auto', 'gmpy', 'gmpy2'):
return None
gmpy = import_module('gmpy2', min_module_version=_GMPY2_MIN_VERSION,
module_version_attr='version', module_version_attr_call_args=())
if sympy_ground_types != 'auto' and gmpy is None:
warn("gmpy2 library is not installed, switching to 'python' ground types")
return gmpy
#
# SYMPY_GROUND_TYPES can be flint, gmpy, gmpy2, python or auto (default)
#
_SYMPY_GROUND_TYPES = os.environ.get('SYMPY_GROUND_TYPES', 'auto').lower()
_flint = None
_gmpy = None
#
# First handle auto-detection of flint/gmpy2. We will prefer flint if available
# or otherwise gmpy2 if available and then lastly the python types.
#
if _SYMPY_GROUND_TYPES in ('auto', 'flint'):
_flint = _get_flint(_SYMPY_GROUND_TYPES)
if _flint is not None:
_SYMPY_GROUND_TYPES = 'flint'
else:
_SYMPY_GROUND_TYPES = 'auto'
if _SYMPY_GROUND_TYPES in ('auto', 'gmpy', 'gmpy2'):
_gmpy = _get_gmpy2(_SYMPY_GROUND_TYPES)
if _gmpy is not None:
_SYMPY_GROUND_TYPES = 'gmpy'
else:
_SYMPY_GROUND_TYPES = 'python'
if _SYMPY_GROUND_TYPES not in ('flint', 'gmpy', 'python'):
warn("SYMPY_GROUND_TYPES environment variable unrecognised. "
"Should be 'auto', 'flint', 'gmpy', 'gmpy2' or 'python'.")
_SYMPY_GROUND_TYPES = 'python'
#
# At this point _SYMPY_GROUND_TYPES is either flint, gmpy or python. The blocks
# below define the values exported by this module in each case.
#
#
# In gmpy2 and flint, there are functions that take a long (or unsigned long)
# argument. That is, it is not possible to input a value larger than that.
#
LONG_MAX = (1 << (8*sizeof(c_long) - 1)) - 1
#
# Type checkers are confused by what SYMPY_INTS is. There may be a better type
# hint for this like Type[Integral] or something.
#
SYMPY_INTS: tTuple[Type, ...]
if _SYMPY_GROUND_TYPES == 'gmpy':
assert _gmpy is not None
flint = None
gmpy = _gmpy
HAS_GMPY = 2
GROUND_TYPES = 'gmpy'
SYMPY_INTS = (int, type(gmpy.mpz(0)))
MPZ = gmpy.mpz
MPQ = gmpy.mpq
# bit_scan1 = gmpy.bit_scan1
# bit_scan0 = gmpy.bit_scan0
# remove = gmpy.remove
# factorial = gmpy.fac
# sqrt = gmpy.isqrt
# is_square = gmpy.is_square
# sqrtrem = gmpy.isqrt_rem
# gcd = gmpy.gcd
# lcm = gmpy.lcm
gcdext = gmpy.gcdext
# invert = gmpy.invert
# legendre = gmpy.legendre
# jacobi = gmpy.jacobi
# kronecker = gmpy.kronecker
# def iroot(x, n):
# # In the latest gmpy2, the threshold for n is ULONG_MAX,
# # but adjust to the older one.
# if n <= LONG_MAX:
# return gmpy.iroot(x, n)
# return python_iroot(x, n)
# is_fermat_prp = gmpy.is_fermat_prp
# is_euler_prp = gmpy.is_euler_prp
# is_strong_prp = gmpy.is_strong_prp
# is_fibonacci_prp = gmpy.is_fibonacci_prp
# is_lucas_prp = gmpy.is_lucas_prp
# is_selfridge_prp = gmpy.is_selfridge_prp
# is_strong_lucas_prp = gmpy.is_strong_lucas_prp
# is_strong_selfridge_prp = gmpy.is_strong_selfridge_prp
# is_bpsw_prp = gmpy.is_bpsw_prp
# is_strong_bpsw_prp = gmpy.is_strong_bpsw_prp
elif _SYMPY_GROUND_TYPES == 'flint':
assert _flint is not None
flint = _flint
gmpy = None
HAS_GMPY = 0
GROUND_TYPES = 'flint'
SYMPY_INTS = (int, flint.fmpz) # type: ignore
MPZ = flint.fmpz # type: ignore
MPQ = flint.fmpq # type: ignore
# bit_scan1 = python_bit_scan1
# bit_scan0 = python_bit_scan0
# remove = python_remove
# factorial = python_factorial
# def sqrt(x):
# return flint.fmpz(x).isqrt()
# def is_square(x):
# if x < 0:
# return False
# return flint.fmpz(x).sqrtrem()[1] == 0
# def sqrtrem(x):
# return flint.fmpz(x).sqrtrem()
# def gcd(*args):
# return reduce(flint.fmpz.gcd, args, flint.fmpz(0))
# def lcm(*args):
# return reduce(flint.fmpz.lcm, args, flint.fmpz(1))
gcdext = python_gcdext
# invert = python_invert
# legendre = python_legendre
# def jacobi(x, y):
# if y <= 0 or not y % 2:
# raise ValueError("y should be an odd positive integer")
# return flint.fmpz(x).jacobi(y)
# kronecker = python_kronecker
# def iroot(x, n):
# if n <= LONG_MAX:
# y = flint.fmpz(x).root(n)
# return y, y**n == x
# return python_iroot(x, n)
# is_fermat_prp = python_is_fermat_prp
# is_euler_prp = python_is_euler_prp
# is_strong_prp = python_is_strong_prp
# is_fibonacci_prp = python_is_fibonacci_prp
# is_lucas_prp = python_is_lucas_prp
# is_selfridge_prp = python_is_selfridge_prp
# is_strong_lucas_prp = python_is_strong_lucas_prp
# is_strong_selfridge_prp = python_is_strong_selfridge_prp
# is_bpsw_prp = python_is_bpsw_prp
# is_strong_bpsw_prp = python_is_strong_bpsw_prp
elif _SYMPY_GROUND_TYPES == 'python':
flint = None
gmpy = None
HAS_GMPY = 0
GROUND_TYPES = 'python'
SYMPY_INTS = (int,)
MPZ = int
MPQ = PythonMPQ
# bit_scan1 = python_bit_scan1
# bit_scan0 = python_bit_scan0
# remove = python_remove
# factorial = python_factorial
# sqrt = python_sqrt
# is_square = python_is_square
# sqrtrem = python_sqrtrem
# gcd = python_gcd
# lcm = python_lcm
gcdext = python_gcdext
# invert = python_invert
# legendre = python_legendre
# jacobi = python_jacobi
# kronecker = python_kronecker
# iroot = python_iroot
# is_fermat_prp = python_is_fermat_prp
# is_euler_prp = python_is_euler_prp
# is_strong_prp = python_is_strong_prp
# is_fibonacci_prp = python_is_fibonacci_prp
# is_lucas_prp = python_is_lucas_prp
# is_selfridge_prp = python_is_selfridge_prp
# is_strong_lucas_prp = python_is_strong_lucas_prp
# is_strong_selfridge_prp = python_is_strong_selfridge_prp
# is_bpsw_prp = python_is_bpsw_prp
# is_strong_bpsw_prp = python_is_strong_bpsw_prp
else:
assert False

View File

@@ -0,0 +1,187 @@
"""Tools to assist importing optional external modules."""
import sys
import re
# Override these in the module to change the default warning behavior.
# For example, you might set both to False before running the tests so that
# warnings are not printed to the console, or set both to True for debugging.
WARN_NOT_INSTALLED = None # Default is False
WARN_OLD_VERSION = None # Default is True
def __sympy_debug():
# helper function from sympy/__init__.py
# We don't just import SYMPY_DEBUG from that file because we don't want to
# import all of SymPy just to use this module.
import os
debug_str = os.getenv('SYMPY_DEBUG', 'False')
if debug_str in ('True', 'False'):
return eval(debug_str)
else:
raise RuntimeError("unrecognized value for SYMPY_DEBUG: %s" %
debug_str)
if __sympy_debug():
WARN_OLD_VERSION = True
WARN_NOT_INSTALLED = True
_component_re = re.compile(r'(\d+ | [a-z]+ | \.)', re.VERBOSE)
def version_tuple(vstring):
# Parse a version string to a tuple e.g. '1.2' -> (1, 2)
# Simplified from distutils.version.LooseVersion which was deprecated in
# Python 3.10.
components = []
for x in _component_re.split(vstring):
if x and x != '.':
try:
x = int(x)
except ValueError:
pass
components.append(x)
return tuple(components)
def import_module(module, min_module_version=None, min_python_version=None,
warn_not_installed=None, warn_old_version=None,
module_version_attr='__version__', module_version_attr_call_args=None,
import_kwargs={}, catch=()):
"""
Import and return a module if it is installed.
If the module is not installed, it returns None.
A minimum version for the module can be given as the keyword argument
min_module_version. This should be comparable against the module version.
By default, module.__version__ is used to get the module version. To
override this, set the module_version_attr keyword argument. If the
attribute of the module to get the version should be called (e.g.,
module.version()), then set module_version_attr_call_args to the args such
that module.module_version_attr(*module_version_attr_call_args) returns the
module's version.
If the module version is less than min_module_version using the Python <
comparison, None will be returned, even if the module is installed. You can
use this to keep from importing an incompatible older version of a module.
You can also specify a minimum Python version by using the
min_python_version keyword argument. This should be comparable against
sys.version_info.
If the keyword argument warn_not_installed is set to True, the function will
emit a UserWarning when the module is not installed.
If the keyword argument warn_old_version is set to True, the function will
emit a UserWarning when the library is installed, but cannot be imported
because of the min_module_version or min_python_version options.
Note that because of the way warnings are handled, a warning will be
emitted for each module only once. You can change the default warning
behavior by overriding the values of WARN_NOT_INSTALLED and WARN_OLD_VERSION
in sympy.external.importtools. By default, WARN_NOT_INSTALLED is False and
WARN_OLD_VERSION is True.
This function uses __import__() to import the module. To pass additional
options to __import__(), use the import_kwargs keyword argument. For
example, to import a submodule A.B, you must pass a nonempty fromlist option
to __import__. See the docstring of __import__().
This catches ImportError to determine if the module is not installed. To
catch additional errors, pass them as a tuple to the catch keyword
argument.
Examples
========
>>> from sympy.external import import_module
>>> numpy = import_module('numpy')
>>> numpy = import_module('numpy', min_python_version=(2, 7),
... warn_old_version=False)
>>> numpy = import_module('numpy', min_module_version='1.5',
... warn_old_version=False) # numpy.__version__ is a string
>>> # gmpy does not have __version__, but it does have gmpy.version()
>>> gmpy = import_module('gmpy', min_module_version='1.14',
... module_version_attr='version', module_version_attr_call_args=(),
... warn_old_version=False)
>>> # To import a submodule, you must pass a nonempty fromlist to
>>> # __import__(). The values do not matter.
>>> p3 = import_module('mpl_toolkits.mplot3d',
... import_kwargs={'fromlist':['something']})
>>> # matplotlib.pyplot can raise RuntimeError when the display cannot be opened
>>> matplotlib = import_module('matplotlib',
... import_kwargs={'fromlist':['pyplot']}, catch=(RuntimeError,))
"""
# keyword argument overrides default, and global variable overrides
# keyword argument.
warn_old_version = (WARN_OLD_VERSION if WARN_OLD_VERSION is not None
else warn_old_version or True)
warn_not_installed = (WARN_NOT_INSTALLED if WARN_NOT_INSTALLED is not None
else warn_not_installed or False)
import warnings
# Check Python first so we don't waste time importing a module we can't use
if min_python_version:
if sys.version_info < min_python_version:
if warn_old_version:
warnings.warn("Python version is too old to use %s "
"(%s or newer required)" % (
module, '.'.join(map(str, min_python_version))),
UserWarning, stacklevel=2)
return
try:
mod = __import__(module, **import_kwargs)
## there's something funny about imports with matplotlib and py3k. doing
## from matplotlib import collections
## gives python's stdlib collections module. explicitly re-importing
## the module fixes this.
from_list = import_kwargs.get('fromlist', ())
for submod in from_list:
if submod == 'collections' and mod.__name__ == 'matplotlib':
__import__(module + '.' + submod)
except ImportError:
if warn_not_installed:
warnings.warn("%s module is not installed" % module, UserWarning,
stacklevel=2)
return
except catch as e:
if warn_not_installed:
warnings.warn(
"%s module could not be used (%s)" % (module, repr(e)),
stacklevel=2)
return
if min_module_version:
modversion = getattr(mod, module_version_attr)
if module_version_attr_call_args is not None:
modversion = modversion(*module_version_attr_call_args)
if version_tuple(modversion) < version_tuple(min_module_version):
if warn_old_version:
# Attempt to create a pretty string version of the version
if isinstance(min_module_version, str):
verstr = min_module_version
elif isinstance(min_module_version, (tuple, list)):
verstr = '.'.join(map(str, min_module_version))
else:
# Either don't know what this is. Hopefully
# it's something that has a nice str version, like an int.
verstr = str(min_module_version)
warnings.warn("%s version is too old to use "
"(%s or newer required)" % (module, verstr),
UserWarning, stacklevel=2)
return
return mod

View File

@@ -0,0 +1,637 @@
# sympy.external.ntheory
#
# This module provides pure Python implementations of some number theory
# functions that are alternately used from gmpy2 if it is installed.
import sys
import math
# import mpmath.libmp as mlib
_small_trailing = [0] * 256
for j in range(1, 8):
_small_trailing[1 << j :: 1 << (j + 1)] = [j] * (1 << (7 - j))
# def bit_scan1(x, n=0):
# if not x:
# return
# x = abs(x >> n)
# low_byte = x & 0xFF
# if low_byte:
# return _small_trailing[low_byte] + n
# t = 8 + n
# x >>= 8
# # 2**m is quick for z up through 2**30
# z = x.bit_length() - 1
# if x == 1 << z:
# return z + t
# if z < 300:
# # fixed 8-byte reduction
# while not x & 0xFF:
# x >>= 8
# t += 8
# else:
# # binary reduction important when there might be a large
# # number of trailing 0s
# p = z >> 1
# while not x & 0xFF:
# while x & ((1 << p) - 1):
# p >>= 1
# x >>= p
# t += p
# return t + _small_trailing[x & 0xFF]
# def bit_scan0(x, n=0):
# return bit_scan1(x + (1 << n), n)
# def remove(x, f):
# if f < 2:
# raise ValueError("factor must be > 1")
# if x == 0:
# return 0, 0
# if f == 2:
# b = bit_scan1(x)
# return x >> b, b
# m = 0
# y, rem = divmod(x, f)
# while not rem:
# x = y
# m += 1
# if m > 5:
# pow_list = [f**2]
# while pow_list:
# _f = pow_list[-1]
# y, rem = divmod(x, _f)
# if not rem:
# m += 1 << len(pow_list)
# x = y
# pow_list.append(_f**2)
# else:
# pow_list.pop()
# y, rem = divmod(x, f)
# return x, m
# def factorial(x):
# """Return x!."""
# return int(mlib.ifac(int(x)))
# def sqrt(x):
# """Integer square root of x."""
# return int(mlib.isqrt(int(x)))
# def sqrtrem(x):
# """Integer square root of x and remainder."""
# s, r = mlib.sqrtrem(int(x))
# return (int(s), int(r))
# if sys.version_info[:2] >= (3, 9):
# # As of Python 3.9 these can take multiple arguments
# gcd = math.gcd
# lcm = math.lcm
# else:
# # Until python 3.8 is no longer supported
# from functools import reduce
# def gcd(*args):
# """gcd of multiple integers."""
# return reduce(math.gcd, args, 0)
# def lcm(*args):
# """lcm of multiple integers."""
# if 0 in args:
# return 0
# return reduce(lambda x, y: x*y//math.gcd(x, y), args, 1)
def _sign(n):
if n < 0:
return -1, -n
return 1, n
def gcdext(a, b):
if not a or not b:
g = abs(a) or abs(b)
if not g:
return (0, 0, 0)
return (g, a // g, b // g)
x_sign, a = _sign(a)
y_sign, b = _sign(b)
x, r = 1, 0
y, s = 0, 1
while b:
q, c = divmod(a, b)
a, b = b, c
x, r = r, x - q*r
y, s = s, y - q*s
return (a, x * x_sign, y * y_sign)
# def is_square(x):
# """Return True if x is a square number."""
# if x < 0:
# return False
# # Note that the possible values of y**2 % n for a given n are limited.
# # For example, when n=4, y**2 % n can only take 0 or 1.
# # In other words, if x % 4 is 2 or 3, then x is not a square number.
# # Mathematically, it determines if it belongs to the set {y**2 % n},
# # but implementationally, it can be realized as a logical conjunction
# # with an n-bit integer.
# # see https://mersenneforum.org/showpost.php?p=110896
# # def magic(n):
# # s = {y**2 % n for y in range(n)}
# # s = set(range(n)) - s
# # return sum(1 << bit for bit in s)
# # >>> print(hex(magic(128)))
# # 0xfdfdfdedfdfdfdecfdfdfdedfdfcfdec
# # >>> print(hex(magic(99)))
# # 0x5f6f9ffb6fb7ddfcb75befdec
# # >>> print(hex(magic(91)))
# # 0x6fd1bfcfed5f3679d3ebdec
# # >>> print(hex(magic(85)))
# # 0xdef9ae771ffe3b9d67dec
# if 0xfdfdfdedfdfdfdecfdfdfdedfdfcfdec & (1 << (x & 127)):
# return False # e.g. 2, 3
# m = x % 765765 # 765765 = 99 * 91 * 85
# if 0x5f6f9ffb6fb7ddfcb75befdec & (1 << (m % 99)):
# return False # e.g. 17, 68
# if 0x6fd1bfcfed5f3679d3ebdec & (1 << (m % 91)):
# return False # e.g. 97, 388
# if 0xdef9ae771ffe3b9d67dec & (1 << (m % 85)):
# return False # e.g. 793, 1408
# return mlib.sqrtrem(int(x))[1] == 0
# def invert(x, m):
# """Modular inverse of x modulo m.
# Returns y such that x*y == 1 mod m.
# Uses ``math.pow`` but reproduces the behaviour of ``gmpy2.invert``
# which raises ZeroDivisionError if no inverse exists.
# """
# try:
# return pow(x, -1, m)
# except ValueError:
# raise ZeroDivisionError("invert() no inverse exists")
# def legendre(x, y):
# """Legendre symbol (x / y).
# Following the implementation of gmpy2,
# the error is raised only when y is an even number.
# """
# if y <= 0 or not y % 2:
# raise ValueError("y should be an odd prime")
# x %= y
# if not x:
# return 0
# if pow(x, (y - 1) // 2, y) == 1:
# return 1
# return -1
# def jacobi(x, y):
# """Jacobi symbol (x / y)."""
# if y <= 0 or not y % 2:
# raise ValueError("y should be an odd positive integer")
# x %= y
# if not x:
# return int(y == 1)
# if y == 1 or x == 1:
# return 1
# if gcd(x, y) != 1:
# return 0
# j = 1
# while x != 0:
# while x % 2 == 0 and x > 0:
# x >>= 1
# if y % 8 in [3, 5]:
# j = -j
# x, y = y, x
# if x % 4 == y % 4 == 3:
# j = -j
# x %= y
# return j
# def kronecker(x, y):
# """Kronecker symbol (x / y)."""
# if gcd(x, y) != 1:
# return 0
# if y == 0:
# return 1
# sign = -1 if y < 0 and x < 0 else 1
# y = abs(y)
# s = bit_scan1(y)
# y >>= s
# if s % 2 and x % 8 in [3, 5]:
# sign = -sign
# return sign * jacobi(x, y)
# def iroot(y, n):
# if y < 0:
# raise ValueError("y must be nonnegative")
# if n < 1:
# raise ValueError("n must be positive")
# if y in (0, 1):
# return y, True
# if n == 1:
# return y, True
# if n == 2:
# x, rem = mlib.sqrtrem(y)
# return int(x), not rem
# if n >= y.bit_length():
# return 1, False
# # Get initial estimate for Newton's method. Care must be taken to
# # avoid overflow
# try:
# guess = int(y**(1./n) + 0.5)
# except OverflowError:
# exp = math.log2(y)/n
# if exp > 53:
# shift = int(exp - 53)
# guess = int(2.0**(exp - shift) + 1) << shift
# else:
# guess = int(2.0**exp)
# if guess > 2**50:
# # Newton iteration
# xprev, x = -1, guess
# while 1:
# t = x**(n - 1)
# xprev, x = x, ((n - 1)*x + y//t)//n
# if abs(x - xprev) < 2:
# break
# else:
# x = guess
# # Compensate
# t = x**n
# while t < y:
# x += 1
# t = x**n
# while t > y:
# x -= 1
# t = x**n
# return x, t == y
# def is_fermat_prp(n, a):
# if a < 2:
# raise ValueError("is_fermat_prp() requires 'a' greater than or equal to 2")
# if n < 1:
# raise ValueError("is_fermat_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# a %= n
# if gcd(n, a) != 1:
# raise ValueError("is_fermat_prp() requires gcd(n,a) == 1")
# return pow(a, n - 1, n) == 1
# def is_euler_prp(n, a):
# if a < 2:
# raise ValueError("is_euler_prp() requires 'a' greater than or equal to 2")
# if n < 1:
# raise ValueError("is_euler_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# a %= n
# if gcd(n, a) != 1:
# raise ValueError("is_euler_prp() requires gcd(n,a) == 1")
# return pow(a, n >> 1, n) == jacobi(a, n) % n
# def _is_strong_prp(n, a):
# s = bit_scan1(n - 1)
# a = pow(a, n >> s, n)
# if a == 1 or a == n - 1:
# return True
# for _ in range(s - 1):
# a = pow(a, 2, n)
# if a == n - 1:
# return True
# if a == 1:
# return False
# return False
# def is_strong_prp(n, a):
# if a < 2:
# raise ValueError("is_strong_prp() requires 'a' greater than or equal to 2")
# if n < 1:
# raise ValueError("is_strong_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# a %= n
# if gcd(n, a) != 1:
# raise ValueError("is_strong_prp() requires gcd(n,a) == 1")
# return _is_strong_prp(n, a)
# def _lucas_sequence(n, P, Q, k):
# r"""Return the modular Lucas sequence (U_k, V_k, Q_k).
# Explanation
# ===========
# Given a Lucas sequence defined by P, Q, returns the kth values for
# U and V, along with Q^k, all modulo n. This is intended for use with
# possibly very large values of n and k, where the combinatorial functions
# would be completely unusable.
# .. math ::
# U_k = \begin{cases}
# 0 & \text{if } k = 0\\
# 1 & \text{if } k = 1\\
# PU_{k-1} - QU_{k-2} & \text{if } k > 1
# \end{cases}\\
# V_k = \begin{cases}
# 2 & \text{if } k = 0\\
# P & \text{if } k = 1\\
# PV_{k-1} - QV_{k-2} & \text{if } k > 1
# \end{cases}
# The modular Lucas sequences are used in numerous places in number theory,
# especially in the Lucas compositeness tests and the various n + 1 proofs.
# Parameters
# ==========
# n : int
# n is an odd number greater than or equal to 3
# P : int
# Q : int
# D determined by D = P**2 - 4*Q is non-zero
# k : int
# k is a nonnegative integer
# Returns
# =======
# U, V, Qk : (int, int, int)
# `(U_k \bmod{n}, V_k \bmod{n}, Q^k \bmod{n})`
# Examples
# ========
# >>> from sympy.external.ntheory import _lucas_sequence
# >>> N = 10**2000 + 4561
# >>> sol = U, V, Qk = _lucas_sequence(N, 3, 1, N//2); sol
# (0, 2, 1)
# References
# ==========
# .. [1] https://en.wikipedia.org/wiki/Lucas_sequence
# """
# if k == 0:
# return (0, 2, 1)
# D = P**2 - 4*Q
# U = 1
# V = P
# Qk = Q % n
# if Q == 1:
# # Optimization for extra strong tests.
# for b in bin(k)[3:]:
# U = (U*V) % n
# V = (V*V - 2) % n
# if b == "1":
# U, V = U*P + V, V*P + U*D
# if U & 1:
# U += n
# if V & 1:
# V += n
# U, V = U >> 1, V >> 1
# elif P == 1 and Q == -1:
# # Small optimization for 50% of Selfridge parameters.
# for b in bin(k)[3:]:
# U = (U*V) % n
# if Qk == 1:
# V = (V*V - 2) % n
# else:
# V = (V*V + 2) % n
# Qk = 1
# if b == "1":
# # new_U = (U + V) // 2
# # new_V = (5*U + V) // 2 = 2*U + new_U
# U, V = U + V, U << 1
# if U & 1:
# U += n
# U >>= 1
# V += U
# Qk = -1
# Qk %= n
# elif P == 1:
# for b in bin(k)[3:]:
# U = (U*V) % n
# V = (V*V - 2*Qk) % n
# Qk *= Qk
# if b == "1":
# # new_U = (U + V) // 2
# # new_V = new_U - 2*Q*U
# U, V = U + V, (Q*U) << 1
# if U & 1:
# U += n
# U >>= 1
# V = U - V
# Qk *= Q
# Qk %= n
# else:
# # The general case with any P and Q.
# for b in bin(k)[3:]:
# U = (U*V) % n
# V = (V*V - 2*Qk) % n
# Qk *= Qk
# if b == "1":
# U, V = U*P + V, V*P + U*D
# if U & 1:
# U += n
# if V & 1:
# V += n
# U, V = U >> 1, V >> 1
# Qk *= Q
# Qk %= n
# return (U % n, V % n, Qk)
# def is_fibonacci_prp(n, p, q):
# d = p**2 - 4*q
# if d == 0 or p <= 0 or q not in [1, -1]:
# raise ValueError("invalid values for p,q in is_fibonacci_prp()")
# if n < 1:
# raise ValueError("is_fibonacci_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# return _lucas_sequence(n, p, q, n)[1] == p % n
# def is_lucas_prp(n, p, q):
# d = p**2 - 4*q
# if d == 0:
# raise ValueError("invalid values for p,q in is_lucas_prp()")
# if n < 1:
# raise ValueError("is_lucas_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# if gcd(n, q*d) not in [1, n]:
# raise ValueError("is_lucas_prp() requires gcd(n,2*q*D) == 1")
# return _lucas_sequence(n, p, q, n - jacobi(d, n))[0] == 0
# def _is_selfridge_prp(n):
# """Lucas compositeness test with the Selfridge parameters for n.
# Explanation
# ===========
# The Lucas compositeness test checks whether n is a prime number.
# The test can be run with arbitrary parameters ``P`` and ``Q``, which also change the performance of the test.
# So, which parameters are most effective for running the Lucas compositeness test?
# As an algorithm for determining ``P`` and ``Q``, Selfridge proposed method A [1]_ page 1401
# (Since two methods were proposed, referred to simply as A and B in the paper,
# we will refer to one of them as "method A").
# method A fixes ``P = 1``. Then, ``D`` defined by ``D = P**2 - 4Q`` is varied from 5, -7, 9, -11, 13, and so on,
# with the first ``D`` being ``jacobi(D, n) == -1``. Once ``D`` is determined,
# ``Q`` is determined to be ``(P**2 - D)//4``.
# References
# ==========
# .. [1] Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes,
# Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417,
# https://doi.org/10.1090%2FS0025-5718-1980-0583518-6
# http://mpqs.free.fr/LucasPseudoprimes.pdf
# """
# for D in range(5, 1_000_000, 2):
# if D & 2: # if D % 4 == 3
# D = -D
# j = jacobi(D, n)
# if j == -1:
# return _lucas_sequence(n, 1, (1-D) // 4, n + 1)[0] == 0
# if j == 0 and D % n:
# return False
# # When j == -1 is hard to find, suspect a square number
# if D == 13 and is_square(n):
# return False
# raise ValueError("appropriate value for D cannot be found in is_selfridge_prp()")
# def is_selfridge_prp(n):
# if n < 1:
# raise ValueError("is_selfridge_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# return _is_selfridge_prp(n)
# def is_strong_lucas_prp(n, p, q):
# D = p**2 - 4*q
# if D == 0:
# raise ValueError("invalid values for p,q in is_strong_lucas_prp()")
# if n < 1:
# raise ValueError("is_selfridge_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# if gcd(n, q*D) not in [1, n]:
# raise ValueError("is_strong_lucas_prp() requires gcd(n,2*q*D) == 1")
# j = jacobi(D, n)
# s = bit_scan1(n - j)
# U, V, Qk = _lucas_sequence(n, p, q, (n - j) >> s)
# if U == 0 or V == 0:
# return True
# for _ in range(s - 1):
# V = (V*V - 2*Qk) % n
# if V == 0:
# return True
# Qk = pow(Qk, 2, n)
# return False
# def _is_strong_selfridge_prp(n):
# for D in range(5, 1_000_000, 2):
# if D & 2: # if D % 4 == 3
# D = -D
# j = jacobi(D, n)
# if j == -1:
# s = bit_scan1(n + 1)
# U, V, Qk = _lucas_sequence(n, 1, (1-D) // 4, (n + 1) >> s)
# if U == 0 or V == 0:
# return True
# for _ in range(s - 1):
# V = (V*V - 2*Qk) % n
# if V == 0:
# return True
# Qk = pow(Qk, 2, n)
# return False
# if j == 0 and D % n:
# return False
# # When j == -1 is hard to find, suspect a square number
# if D == 13 and is_square(n):
# return False
# raise ValueError("appropriate value for D cannot be found in is_strong_selfridge_prp()")
# def is_strong_selfridge_prp(n):
# if n < 1:
# raise ValueError("is_strong_selfridge_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# return _is_strong_selfridge_prp(n)
# def is_bpsw_prp(n):
# if n < 1:
# raise ValueError("is_bpsw_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# return _is_strong_prp(n, 2) and _is_selfridge_prp(n)
# def is_strong_bpsw_prp(n):
# if n < 1:
# raise ValueError("is_strong_bpsw_prp() requires 'n' be greater than 0")
# if n == 1:
# return False
# if n % 2 == 0:
# return n == 2
# return _is_strong_prp(n, 2) and _is_strong_selfridge_prp(n)

View File

@@ -0,0 +1,341 @@
"""
PythonMPQ: Rational number type based on Python integers.
This class is intended as a pure Python fallback for when gmpy2 is not
installed. If gmpy2 is installed then its mpq type will be used instead. The
mpq type is around 20x faster. We could just use the stdlib Fraction class
here but that is slower:
from fractions import Fraction
from sympy.external.pythonmpq import PythonMPQ
nums = range(1000)
dens = range(5, 1005)
rats = [Fraction(n, d) for n, d in zip(nums, dens)]
sum(rats) # <--- 24 milliseconds
rats = [PythonMPQ(n, d) for n, d in zip(nums, dens)]
sum(rats) # <--- 7 milliseconds
Both mpq and Fraction have some awkward features like the behaviour of
division with // and %:
>>> from fractions import Fraction
>>> Fraction(2, 3) % Fraction(1, 4)
1/6
For the QQ domain we do not want this behaviour because there should be no
remainder when dividing rational numbers. SymPy does not make use of this
aspect of mpq when gmpy2 is installed. Since this class is a fallback for that
case we do not bother implementing e.g. __mod__ so that we can be sure we
are not using it when gmpy2 is installed either.
"""
import operator
from math import gcd
from decimal import Decimal
from fractions import Fraction
import sys
from typing import Tuple as tTuple, Type
# Used for __hash__
_PyHASH_MODULUS = sys.hash_info.modulus
_PyHASH_INF = sys.hash_info.inf
class PythonMPQ:
"""Rational number implementation that is intended to be compatible with
gmpy2's mpq.
Also slightly faster than fractions.Fraction.
PythonMPQ should be treated as immutable although no effort is made to
prevent mutation (since that might slow down calculations).
"""
__slots__ = ('numerator', 'denominator')
def __new__(cls, numerator, denominator=None):
"""Construct PythonMPQ with gcd computation and checks"""
if denominator is not None:
#
# PythonMPQ(n, d): require n and d to be int and d != 0
#
if isinstance(numerator, int) and isinstance(denominator, int):
# This is the slow part:
divisor = gcd(numerator, denominator)
numerator //= divisor
denominator //= divisor
return cls._new_check(numerator, denominator)
else:
#
# PythonMPQ(q)
#
# Here q can be PythonMPQ, int, Decimal, float, Fraction or str
#
if isinstance(numerator, int):
return cls._new(numerator, 1)
elif isinstance(numerator, PythonMPQ):
return cls._new(numerator.numerator, numerator.denominator)
# Let Fraction handle Decimal/float conversion and str parsing
if isinstance(numerator, (Decimal, float, str)):
numerator = Fraction(numerator)
if isinstance(numerator, Fraction):
return cls._new(numerator.numerator, numerator.denominator)
#
# Reject everything else. This is more strict than mpq which allows
# things like mpq(Fraction, Fraction) or mpq(Decimal, any). The mpq
# behaviour is somewhat inconsistent so we choose to accept only a
# more strict subset of what mpq allows.
#
raise TypeError("PythonMPQ() requires numeric or string argument")
@classmethod
def _new_check(cls, numerator, denominator):
"""Construct PythonMPQ, check divide by zero and canonicalize signs"""
if not denominator:
raise ZeroDivisionError(f'Zero divisor {numerator}/{denominator}')
elif denominator < 0:
numerator = -numerator
denominator = -denominator
return cls._new(numerator, denominator)
@classmethod
def _new(cls, numerator, denominator):
"""Construct PythonMPQ efficiently (no checks)"""
obj = super().__new__(cls)
obj.numerator = numerator
obj.denominator = denominator
return obj
def __int__(self):
"""Convert to int (truncates towards zero)"""
p, q = self.numerator, self.denominator
if p < 0:
return -(-p//q)
return p//q
def __float__(self):
"""Convert to float (approximately)"""
return self.numerator / self.denominator
def __bool__(self):
"""True/False if nonzero/zero"""
return bool(self.numerator)
def __eq__(self, other):
"""Compare equal with PythonMPQ, int, float, Decimal or Fraction"""
if isinstance(other, PythonMPQ):
return (self.numerator == other.numerator
and self.denominator == other.denominator)
elif isinstance(other, self._compatible_types):
return self.__eq__(PythonMPQ(other))
else:
return NotImplemented
def __hash__(self):
"""hash - same as mpq/Fraction"""
try:
dinv = pow(self.denominator, -1, _PyHASH_MODULUS)
except ValueError:
hash_ = _PyHASH_INF
else:
hash_ = hash(hash(abs(self.numerator)) * dinv)
result = hash_ if self.numerator >= 0 else -hash_
return -2 if result == -1 else result
def __reduce__(self):
"""Deconstruct for pickling"""
return type(self), (self.numerator, self.denominator)
def __str__(self):
"""Convert to string"""
if self.denominator != 1:
return f"{self.numerator}/{self.denominator}"
else:
return f"{self.numerator}"
def __repr__(self):
"""Convert to string"""
return f"MPQ({self.numerator},{self.denominator})"
def _cmp(self, other, op):
"""Helper for lt/le/gt/ge"""
if not isinstance(other, self._compatible_types):
return NotImplemented
lhs = self.numerator * other.denominator
rhs = other.numerator * self.denominator
return op(lhs, rhs)
def __lt__(self, other):
"""self < other"""
return self._cmp(other, operator.lt)
def __le__(self, other):
"""self <= other"""
return self._cmp(other, operator.le)
def __gt__(self, other):
"""self > other"""
return self._cmp(other, operator.gt)
def __ge__(self, other):
"""self >= other"""
return self._cmp(other, operator.ge)
def __abs__(self):
"""abs(q)"""
return self._new(abs(self.numerator), self.denominator)
def __pos__(self):
"""+q"""
return self
def __neg__(self):
"""-q"""
return self._new(-self.numerator, self.denominator)
def __add__(self, other):
"""q1 + q2"""
if isinstance(other, PythonMPQ):
#
# This is much faster than the naive method used in the stdlib
# fractions module. Not sure where this method comes from
# though...
#
# Compare timings for something like:
# nums = range(1000)
# rats = [PythonMPQ(n, d) for n, d in zip(nums[:-5], nums[5:])]
# sum(rats) # <-- time this
#
ap, aq = self.numerator, self.denominator
bp, bq = other.numerator, other.denominator
g = gcd(aq, bq)
if g == 1:
p = ap*bq + aq*bp
q = bq*aq
else:
q1, q2 = aq//g, bq//g
p, q = ap*q2 + bp*q1, q1*q2
g2 = gcd(p, g)
p, q = (p // g2), q * (g // g2)
elif isinstance(other, int):
p = self.numerator + self.denominator * other
q = self.denominator
else:
return NotImplemented
return self._new(p, q)
def __radd__(self, other):
"""z1 + q2"""
if isinstance(other, int):
p = self.numerator + self.denominator * other
q = self.denominator
return self._new(p, q)
else:
return NotImplemented
def __sub__(self ,other):
"""q1 - q2"""
if isinstance(other, PythonMPQ):
ap, aq = self.numerator, self.denominator
bp, bq = other.numerator, other.denominator
g = gcd(aq, bq)
if g == 1:
p = ap*bq - aq*bp
q = bq*aq
else:
q1, q2 = aq//g, bq//g
p, q = ap*q2 - bp*q1, q1*q2
g2 = gcd(p, g)
p, q = (p // g2), q * (g // g2)
elif isinstance(other, int):
p = self.numerator - self.denominator*other
q = self.denominator
else:
return NotImplemented
return self._new(p, q)
def __rsub__(self, other):
"""z1 - q2"""
if isinstance(other, int):
p = self.denominator * other - self.numerator
q = self.denominator
return self._new(p, q)
else:
return NotImplemented
def __mul__(self, other):
"""q1 * q2"""
if isinstance(other, PythonMPQ):
ap, aq = self.numerator, self.denominator
bp, bq = other.numerator, other.denominator
x1 = gcd(ap, bq)
x2 = gcd(bp, aq)
p, q = ((ap//x1)*(bp//x2), (aq//x2)*(bq//x1))
elif isinstance(other, int):
x = gcd(other, self.denominator)
p = self.numerator*(other//x)
q = self.denominator//x
else:
return NotImplemented
return self._new(p, q)
def __rmul__(self, other):
"""z1 * q2"""
if isinstance(other, int):
x = gcd(self.denominator, other)
p = self.numerator*(other//x)
q = self.denominator//x
return self._new(p, q)
else:
return NotImplemented
def __pow__(self, exp):
"""q ** z"""
p, q = self.numerator, self.denominator
if exp < 0:
p, q, exp = q, p, -exp
return self._new_check(p**exp, q**exp)
def __truediv__(self, other):
"""q1 / q2"""
if isinstance(other, PythonMPQ):
ap, aq = self.numerator, self.denominator
bp, bq = other.numerator, other.denominator
x1 = gcd(ap, bp)
x2 = gcd(bq, aq)
p, q = ((ap//x1)*(bq//x2), (aq//x2)*(bp//x1))
elif isinstance(other, int):
x = gcd(other, self.numerator)
p = self.numerator//x
q = self.denominator*(other//x)
else:
return NotImplemented
return self._new_check(p, q)
def __rtruediv__(self, other):
"""z / q"""
if isinstance(other, int):
x = gcd(self.numerator, other)
p = self.denominator*(other//x)
q = self.numerator//x
return self._new_check(p, q)
else:
return NotImplemented
_compatible_types: tTuple[Type, ...] = ()
#
# These are the types that PythonMPQ will interoperate with for operations
# and comparisons such as ==, + etc. We define this down here so that we can
# include PythonMPQ in the list as well.
#
PythonMPQ._compatible_types = (PythonMPQ, int, Decimal, Fraction)