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 |