"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/identification.py" between
mpmath-0.18.tar.gz and mpmath-0.19.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

identification.py  (mpmath-0.18):identification.py  (mpmath-0.19)
skipping to change at line 132 skipping to change at line 132
This is a fairly direct translation to Python of the pseudocode given by This is a fairly direct translation to Python of the pseudocode given by
David Bailey, "The PSLQ Integer Relation Algorithm": David Bailey, "The PSLQ Integer Relation Algorithm":
http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/node3.html http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/node3.html
The present implementation uses fixed-point instead of floating-point The present implementation uses fixed-point instead of floating-point
arithmetic, since this is significantly (about 7x) faster. arithmetic, since this is significantly (about 7x) faster.
""" """
n = len(x) n = len(x)
assert n >= 2 if n < 2:
raise ValueError("n cannot be less than 2")
# At too low precision, the algorithm becomes meaningless # At too low precision, the algorithm becomes meaningless
prec = ctx.prec prec = ctx.prec
assert prec >= 53 if prec < 53:
raise ValueError("prec cannot be less than 53")
if verbose and prec // max(2,n) < 5: if verbose and prec // max(2,n) < 5:
print("Warning: precision for PSLQ may be too low") print("Warning: precision for PSLQ may be too low")
target = int(prec * 0.75) target = int(prec * 0.75)
if tol is None: if tol is None:
tol = ctx.mpf(2)**(-target) tol = ctx.mpf(2)**(-target)
else: else:
tol = ctx.convert(tol) tol = ctx.convert(tol)
skipping to change at line 415 skipping to change at line 417
>>> >>>
Note that the high precision and strict tolerance is necessary Note that the high precision and strict tolerance is necessary
for such high-degree runs, since otherwise unwanted low-accuracy for such high-degree runs, since otherwise unwanted low-accuracy
approximations will be detected. It may also be necessary to set approximations will be detected. It may also be necessary to set
maxsteps high to prevent a premature exit (before the coefficient maxsteps high to prevent a premature exit (before the coefficient
bound has been reached). Running with ``verbose=True`` to get an bound has been reached). Running with ``verbose=True`` to get an
idea what is happening can be useful. idea what is happening can be useful.
""" """
x = ctx.mpf(x) x = ctx.mpf(x)
assert n >= 1 if n < 1:
raise ValueError("n cannot be less than 1")
if x == 0: if x == 0:
return [1, 0] return [1, 0]
xs = [ctx.mpf(1)] xs = [ctx.mpf(1)]
for i in range(1,n+1): for i in range(1,n+1):
xs.append(x**i) xs.append(x**i)
a = ctx.pslq(xs, **kwargs) a = ctx.pslq(xs, **kwargs)
if a is not None: if a is not None:
return a[::-1] return a[::-1]
def fracgcd(p, q): def fracgcd(p, q):
skipping to change at line 699 skipping to change at line 702
in the algebraically simplest form:: in the algebraically simplest form::
>>> identify(sqrt(2)) >>> identify(sqrt(2))
'(sqrt(8)/2)' '(sqrt(8)/2)'
As a solution to both problems, consider using SymPy's As a solution to both problems, consider using SymPy's
:func:`~mpmath.sympify` to convert the formula into a symbolic expression. :func:`~mpmath.sympify` to convert the formula into a symbolic expression.
SymPy can be used to pretty-print or further simplify the formula SymPy can be used to pretty-print or further simplify the formula
symbolically:: symbolically::
>>> from sympy import sympify >>> from sympy import sympify # doctest: +SKIP
>>> sympify(identify(sqrt(2))) >>> sympify(identify(sqrt(2))) # doctest: +SKIP
2**(1/2) 2**(1/2)
Sometimes :func:`~mpmath.identify` can simplify an expression further than Sometimes :func:`~mpmath.identify` can simplify an expression further than
a symbolic algorithm:: a symbolic algorithm::
>>> from sympy import simplify >>> from sympy import simplify # doctest: +SKIP
>>> x = sympify('-1/(-3/2+(1/2)*5**(1/2))*(3/2-1/2*5**(1/2))**(1/2)') >>> x = sympify('-1/(-3/2+(1/2)*5**(1/2))*(3/2-1/2*5**(1/2))**(1/2)') #
>>> x doctest: +SKIP
>>> x # doctest: +SKIP
(3/2 - 5**(1/2)/2)**(-1/2) (3/2 - 5**(1/2)/2)**(-1/2)
>>> x = simplify(x) >>> x = simplify(x) # doctest: +SKIP
>>> x >>> x # doctest: +SKIP
2/(6 - 2*5**(1/2))**(1/2) 2/(6 - 2*5**(1/2))**(1/2)
>>> mp.dps = 30 >>> mp.dps = 30 # doctest: +SKIP
>>> x = sympify(identify(x.evalf(30))) >>> x = sympify(identify(x.evalf(30))) # doctest: +SKIP
>>> x >>> x # doctest: +SKIP
1/2 + 5**(1/2)/2 1/2 + 5**(1/2)/2
(In fact, this functionality is available directly in SymPy as the (In fact, this functionality is available directly in SymPy as the
function :func:`~mpmath.nsimplify`, which is essentially a wrapper for function :func:`~mpmath.nsimplify`, which is essentially a wrapper for
:func:`~mpmath.identify`.) :func:`~mpmath.identify`.)
**Miscellaneous issues and limitations** **Miscellaneous issues and limitations**
The input `x` must be a real number. All base constants must be The input `x` must be a real number. All base constants must be
positive real numbers and must not be rationals or rational linear positive real numbers and must not be rationals or rational linear
skipping to change at line 771 skipping to change at line 774
return "-(%s)" % sol return "-(%s)" % sol
if tol: if tol:
tol = ctx.mpf(tol) tol = ctx.mpf(tol)
else: else:
tol = ctx.eps**0.7 tol = ctx.eps**0.7
M = maxcoeff M = maxcoeff
if constants: if constants:
if isinstance(constants, dict): if isinstance(constants, dict):
constants = [(ctx.mpf(v), name) for (name, v) in constants.items()] constants = [(ctx.mpf(v), name) for (name, v) in sorted(constants.it ems())]
else: else:
namespace = dict((name, getattr(ctx,name)) for name in dir(ctx)) namespace = dict((name, getattr(ctx,name)) for name in dir(ctx))
constants = [(eval(p, namespace), p) for p in constants] constants = [(eval(p, namespace), p) for p in constants]
else: else:
constants = [] constants = []
# We always want to find at least rational terms # We always want to find at least rational terms
if 1 not in [value for (name, value) in constants]: if 1 not in [value for (name, value) in constants]:
constants = [(ctx.mpf(1), '1')] + constants constants = [(ctx.mpf(1), '1')] + constants
 End of changes. 8 change blocks. 
14 lines changed or deleted 18 lines changed or added

Home  |  About  |  All  |  Newest  |  Fossies Dox  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTPS