Coverage for /usr/lib/python3/dist-packages/mpmath/identification.py: 6%
282 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Implements the PSLQ algorithm for integer relation detection,
3and derivative algorithms for constant recognition.
4"""
6from .libmp.backend import xrange
7from .libmp import int_types, sqrt_fixed
9# round to nearest integer (can be done more elegantly...)
10def round_fixed(x, prec):
11 return ((x + (1<<(prec-1))) >> prec) << prec
13class IdentificationMethods(object):
14 pass
17def pslq(ctx, x, tol=None, maxcoeff=1000, maxsteps=100, verbose=False):
18 r"""
19 Given a vector of real numbers `x = [x_0, x_1, ..., x_n]`, ``pslq(x)``
20 uses the PSLQ algorithm to find a list of integers
21 `[c_0, c_1, ..., c_n]` such that
23 .. math ::
25 |c_1 x_1 + c_2 x_2 + ... + c_n x_n| < \mathrm{tol}
27 and such that `\max |c_k| < \mathrm{maxcoeff}`. If no such vector
28 exists, :func:`~mpmath.pslq` returns ``None``. The tolerance defaults to
29 3/4 of the working precision.
31 **Examples**
33 Find rational approximations for `\pi`::
35 >>> from mpmath import *
36 >>> mp.dps = 15; mp.pretty = True
37 >>> pslq([-1, pi], tol=0.01)
38 [22, 7]
39 >>> pslq([-1, pi], tol=0.001)
40 [355, 113]
41 >>> mpf(22)/7; mpf(355)/113; +pi
42 3.14285714285714
43 3.14159292035398
44 3.14159265358979
46 Pi is not a rational number with denominator less than 1000::
48 >>> pslq([-1, pi])
49 >>>
51 To within the standard precision, it can however be approximated
52 by at least one rational number with denominator less than `10^{12}`::
54 >>> p, q = pslq([-1, pi], maxcoeff=10**12)
55 >>> print(p); print(q)
56 238410049439
57 75888275702
58 >>> mpf(p)/q
59 3.14159265358979
61 The PSLQ algorithm can be applied to long vectors. For example,
62 we can investigate the rational (in)dependence of integer square
63 roots::
65 >>> mp.dps = 30
66 >>> pslq([sqrt(n) for n in range(2, 5+1)])
67 >>>
68 >>> pslq([sqrt(n) for n in range(2, 6+1)])
69 >>>
70 >>> pslq([sqrt(n) for n in range(2, 8+1)])
71 [2, 0, 0, 0, 0, 0, -1]
73 **Machin formulas**
75 A famous formula for `\pi` is Machin's,
77 .. math ::
79 \frac{\pi}{4} = 4 \operatorname{acot} 5 - \operatorname{acot} 239
81 There are actually infinitely many formulas of this type. Two
82 others are
84 .. math ::
86 \frac{\pi}{4} = \operatorname{acot} 1
88 \frac{\pi}{4} = 12 \operatorname{acot} 49 + 32 \operatorname{acot} 57
89 + 5 \operatorname{acot} 239 + 12 \operatorname{acot} 110443
91 We can easily verify the formulas using the PSLQ algorithm::
93 >>> mp.dps = 30
94 >>> pslq([pi/4, acot(1)])
95 [1, -1]
96 >>> pslq([pi/4, acot(5), acot(239)])
97 [1, -4, 1]
98 >>> pslq([pi/4, acot(49), acot(57), acot(239), acot(110443)])
99 [1, -12, -32, 5, -12]
101 We could try to generate a custom Machin-like formula by running
102 the PSLQ algorithm with a few inverse cotangent values, for example
103 acot(2), acot(3) ... acot(10). Unfortunately, there is a linear
104 dependence among these values, resulting in only that dependence
105 being detected, with a zero coefficient for `\pi`::
107 >>> pslq([pi] + [acot(n) for n in range(2,11)])
108 [0, 1, -1, 0, 0, 0, -1, 0, 0, 0]
110 We get better luck by removing linearly dependent terms::
112 >>> pslq([pi] + [acot(n) for n in range(2,11) if n not in (3, 5)])
113 [1, -8, 0, 0, 4, 0, 0, 0]
115 In other words, we found the following formula::
117 >>> 8*acot(2) - 4*acot(7)
118 3.14159265358979323846264338328
119 >>> +pi
120 3.14159265358979323846264338328
122 **Algorithm**
124 This is a fairly direct translation to Python of the pseudocode given by
125 David Bailey, "The PSLQ Integer Relation Algorithm":
126 http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/node3.html
128 The present implementation uses fixed-point instead of floating-point
129 arithmetic, since this is significantly (about 7x) faster.
130 """
132 n = len(x)
133 if n < 2:
134 raise ValueError("n cannot be less than 2")
136 # At too low precision, the algorithm becomes meaningless
137 prec = ctx.prec
138 if prec < 53:
139 raise ValueError("prec cannot be less than 53")
141 if verbose and prec // max(2,n) < 5:
142 print("Warning: precision for PSLQ may be too low")
144 target = int(prec * 0.75)
146 if tol is None:
147 tol = ctx.mpf(2)**(-target)
148 else:
149 tol = ctx.convert(tol)
151 extra = 60
152 prec += extra
154 if verbose:
155 print("PSLQ using prec %i and tol %s" % (prec, ctx.nstr(tol)))
157 tol = ctx.to_fixed(tol, prec)
158 assert tol
160 # Convert to fixed-point numbers. The dummy None is added so we can
161 # use 1-based indexing. (This just allows us to be consistent with
162 # Bailey's indexing. The algorithm is 100 lines long, so debugging
163 # a single wrong index can be painful.)
164 x = [None] + [ctx.to_fixed(ctx.mpf(xk), prec) for xk in x]
166 # Sanity check on magnitudes
167 minx = min(abs(xx) for xx in x[1:])
168 if not minx:
169 raise ValueError("PSLQ requires a vector of nonzero numbers")
170 if minx < tol//100:
171 if verbose:
172 print("STOPPING: (one number is too small)")
173 return None
175 g = sqrt_fixed((4<<prec)//3, prec)
176 A = {}
177 B = {}
178 H = {}
179 # Initialization
180 # step 1
181 for i in xrange(1, n+1):
182 for j in xrange(1, n+1):
183 A[i,j] = B[i,j] = (i==j) << prec
184 H[i,j] = 0
185 # step 2
186 s = [None] + [0] * n
187 for k in xrange(1, n+1):
188 t = 0
189 for j in xrange(k, n+1):
190 t += (x[j]**2 >> prec)
191 s[k] = sqrt_fixed(t, prec)
192 t = s[1]
193 y = x[:]
194 for k in xrange(1, n+1):
195 y[k] = (x[k] << prec) // t
196 s[k] = (s[k] << prec) // t
197 # step 3
198 for i in xrange(1, n+1):
199 for j in xrange(i+1, n):
200 H[i,j] = 0
201 if i <= n-1:
202 if s[i]:
203 H[i,i] = (s[i+1] << prec) // s[i]
204 else:
205 H[i,i] = 0
206 for j in range(1, i):
207 sjj1 = s[j]*s[j+1]
208 if sjj1:
209 H[i,j] = ((-y[i]*y[j])<<prec)//sjj1
210 else:
211 H[i,j] = 0
212 # step 4
213 for i in xrange(2, n+1):
214 for j in xrange(i-1, 0, -1):
215 #t = floor(H[i,j]/H[j,j] + 0.5)
216 if H[j,j]:
217 t = round_fixed((H[i,j] << prec)//H[j,j], prec)
218 else:
219 #t = 0
220 continue
221 y[j] = y[j] + (t*y[i] >> prec)
222 for k in xrange(1, j+1):
223 H[i,k] = H[i,k] - (t*H[j,k] >> prec)
224 for k in xrange(1, n+1):
225 A[i,k] = A[i,k] - (t*A[j,k] >> prec)
226 B[k,j] = B[k,j] + (t*B[k,i] >> prec)
227 # Main algorithm
228 for REP in range(maxsteps):
229 # Step 1
230 m = -1
231 szmax = -1
232 for i in range(1, n):
233 h = H[i,i]
234 sz = (g**i * abs(h)) >> (prec*(i-1))
235 if sz > szmax:
236 m = i
237 szmax = sz
238 # Step 2
239 y[m], y[m+1] = y[m+1], y[m]
240 for i in xrange(1,n+1): H[m,i], H[m+1,i] = H[m+1,i], H[m,i]
241 for i in xrange(1,n+1): A[m,i], A[m+1,i] = A[m+1,i], A[m,i]
242 for i in xrange(1,n+1): B[i,m], B[i,m+1] = B[i,m+1], B[i,m]
243 # Step 3
244 if m <= n - 2:
245 t0 = sqrt_fixed((H[m,m]**2 + H[m,m+1]**2)>>prec, prec)
246 # A zero element probably indicates that the precision has
247 # been exhausted. XXX: this could be spurious, due to
248 # using fixed-point arithmetic
249 if not t0:
250 break
251 t1 = (H[m,m] << prec) // t0
252 t2 = (H[m,m+1] << prec) // t0
253 for i in xrange(m, n+1):
254 t3 = H[i,m]
255 t4 = H[i,m+1]
256 H[i,m] = (t1*t3+t2*t4) >> prec
257 H[i,m+1] = (-t2*t3+t1*t4) >> prec
258 # Step 4
259 for i in xrange(m+1, n+1):
260 for j in xrange(min(i-1, m+1), 0, -1):
261 try:
262 t = round_fixed((H[i,j] << prec)//H[j,j], prec)
263 # Precision probably exhausted
264 except ZeroDivisionError:
265 break
266 y[j] = y[j] + ((t*y[i]) >> prec)
267 for k in xrange(1, j+1):
268 H[i,k] = H[i,k] - (t*H[j,k] >> prec)
269 for k in xrange(1, n+1):
270 A[i,k] = A[i,k] - (t*A[j,k] >> prec)
271 B[k,j] = B[k,j] + (t*B[k,i] >> prec)
272 # Until a relation is found, the error typically decreases
273 # slowly (e.g. a factor 1-10) with each step TODO: we could
274 # compare err from two successive iterations. If there is a
275 # large drop (several orders of magnitude), that indicates a
276 # "high quality" relation was detected. Reporting this to
277 # the user somehow might be useful.
278 best_err = maxcoeff<<prec
279 for i in xrange(1, n+1):
280 err = abs(y[i])
281 # Maybe we are done?
282 if err < tol:
283 # We are done if the coefficients are acceptable
284 vec = [int(round_fixed(B[j,i], prec) >> prec) for j in \
285 range(1,n+1)]
286 if max(abs(v) for v in vec) < maxcoeff:
287 if verbose:
288 print("FOUND relation at iter %i/%i, error: %s" % \
289 (REP, maxsteps, ctx.nstr(err / ctx.mpf(2)**prec, 1)))
290 return vec
291 best_err = min(err, best_err)
292 # Calculate a lower bound for the norm. We could do this
293 # more exactly (using the Euclidean norm) but there is probably
294 # no practical benefit.
295 recnorm = max(abs(h) for h in H.values())
296 if recnorm:
297 norm = ((1 << (2*prec)) // recnorm) >> prec
298 norm //= 100
299 else:
300 norm = ctx.inf
301 if verbose:
302 print("%i/%i: Error: %8s Norm: %s" % \
303 (REP, maxsteps, ctx.nstr(best_err / ctx.mpf(2)**prec, 1), norm))
304 if norm >= maxcoeff:
305 break
306 if verbose:
307 print("CANCELLING after step %i/%i." % (REP, maxsteps))
308 print("Could not find an integer relation. Norm bound: %s" % norm)
309 return None
311def findpoly(ctx, x, n=1, **kwargs):
312 r"""
313 ``findpoly(x, n)`` returns the coefficients of an integer
314 polynomial `P` of degree at most `n` such that `P(x) \approx 0`.
315 If no polynomial having `x` as a root can be found,
316 :func:`~mpmath.findpoly` returns ``None``.
318 :func:`~mpmath.findpoly` works by successively calling :func:`~mpmath.pslq` with
319 the vectors `[1, x]`, `[1, x, x^2]`, `[1, x, x^2, x^3]`, ...,
320 `[1, x, x^2, .., x^n]` as input. Keyword arguments given to
321 :func:`~mpmath.findpoly` are forwarded verbatim to :func:`~mpmath.pslq`. In
322 particular, you can specify a tolerance for `P(x)` with ``tol``
323 and a maximum permitted coefficient size with ``maxcoeff``.
325 For large values of `n`, it is recommended to run :func:`~mpmath.findpoly`
326 at high precision; preferably 50 digits or more.
328 **Examples**
330 By default (degree `n = 1`), :func:`~mpmath.findpoly` simply finds a linear
331 polynomial with a rational root::
333 >>> from mpmath import *
334 >>> mp.dps = 15; mp.pretty = True
335 >>> findpoly(0.7)
336 [-10, 7]
338 The generated coefficient list is valid input to ``polyval`` and
339 ``polyroots``::
341 >>> nprint(polyval(findpoly(phi, 2), phi), 1)
342 -2.0e-16
343 >>> for r in polyroots(findpoly(phi, 2)):
344 ... print(r)
345 ...
346 -0.618033988749895
347 1.61803398874989
349 Numbers of the form `m + n \sqrt p` for integers `(m, n, p)` are
350 solutions to quadratic equations. As we find here, `1+\sqrt 2`
351 is a root of the polynomial `x^2 - 2x - 1`::
353 >>> findpoly(1+sqrt(2), 2)
354 [1, -2, -1]
355 >>> findroot(lambda x: x**2 - 2*x - 1, 1)
356 2.4142135623731
358 Despite only containing square roots, the following number results
359 in a polynomial of degree 4::
361 >>> findpoly(sqrt(2)+sqrt(3), 4)
362 [1, 0, -10, 0, 1]
364 In fact, `x^4 - 10x^2 + 1` is the *minimal polynomial* of
365 `r = \sqrt 2 + \sqrt 3`, meaning that a rational polynomial of
366 lower degree having `r` as a root does not exist. Given sufficient
367 precision, :func:`~mpmath.findpoly` will usually find the correct
368 minimal polynomial of a given algebraic number.
370 **Non-algebraic numbers**
372 If :func:`~mpmath.findpoly` fails to find a polynomial with given
373 coefficient size and tolerance constraints, that means no such
374 polynomial exists.
376 We can verify that `\pi` is not an algebraic number of degree 3 with
377 coefficients less than 1000::
379 >>> mp.dps = 15
380 >>> findpoly(pi, 3)
381 >>>
383 It is always possible to find an algebraic approximation of a number
384 using one (or several) of the following methods:
386 1. Increasing the permitted degree
387 2. Allowing larger coefficients
388 3. Reducing the tolerance
390 One example of each method is shown below::
392 >>> mp.dps = 15
393 >>> findpoly(pi, 4)
394 [95, -545, 863, -183, -298]
395 >>> findpoly(pi, 3, maxcoeff=10000)
396 [836, -1734, -2658, -457]
397 >>> findpoly(pi, 3, tol=1e-7)
398 [-4, 22, -29, -2]
400 It is unknown whether Euler's constant is transcendental (or even
401 irrational). We can use :func:`~mpmath.findpoly` to check that if is
402 an algebraic number, its minimal polynomial must have degree
403 at least 7 and a coefficient of magnitude at least 1000000::
405 >>> mp.dps = 200
406 >>> findpoly(euler, 6, maxcoeff=10**6, tol=1e-100, maxsteps=1000)
407 >>>
409 Note that the high precision and strict tolerance is necessary
410 for such high-degree runs, since otherwise unwanted low-accuracy
411 approximations will be detected. It may also be necessary to set
412 maxsteps high to prevent a premature exit (before the coefficient
413 bound has been reached). Running with ``verbose=True`` to get an
414 idea what is happening can be useful.
415 """
416 x = ctx.mpf(x)
417 if n < 1:
418 raise ValueError("n cannot be less than 1")
419 if x == 0:
420 return [1, 0]
421 xs = [ctx.mpf(1)]
422 for i in range(1,n+1):
423 xs.append(x**i)
424 a = ctx.pslq(xs, **kwargs)
425 if a is not None:
426 return a[::-1]
428def fracgcd(p, q):
429 x, y = p, q
430 while y:
431 x, y = y, x % y
432 if x != 1:
433 p //= x
434 q //= x
435 if q == 1:
436 return p
437 return p, q
439def pslqstring(r, constants):
440 q = r[0]
441 r = r[1:]
442 s = []
443 for i in range(len(r)):
444 p = r[i]
445 if p:
446 z = fracgcd(-p,q)
447 cs = constants[i][1]
448 if cs == '1':
449 cs = ''
450 else:
451 cs = '*' + cs
452 if isinstance(z, int_types):
453 if z > 0: term = str(z) + cs
454 else: term = ("(%s)" % z) + cs
455 else:
456 term = ("(%s/%s)" % z) + cs
457 s.append(term)
458 s = ' + '.join(s)
459 if '+' in s or '*' in s:
460 s = '(' + s + ')'
461 return s or '0'
463def prodstring(r, constants):
464 q = r[0]
465 r = r[1:]
466 num = []
467 den = []
468 for i in range(len(r)):
469 p = r[i]
470 if p:
471 z = fracgcd(-p,q)
472 cs = constants[i][1]
473 if isinstance(z, int_types):
474 if abs(z) == 1: t = cs
475 else: t = '%s**%s' % (cs, abs(z))
476 ([num,den][z<0]).append(t)
477 else:
478 t = '%s**(%s/%s)' % (cs, abs(z[0]), z[1])
479 ([num,den][z[0]<0]).append(t)
480 num = '*'.join(num)
481 den = '*'.join(den)
482 if num and den: return "(%s)/(%s)" % (num, den)
483 if num: return num
484 if den: return "1/(%s)" % den
486def quadraticstring(ctx,t,a,b,c):
487 if c < 0:
488 a,b,c = -a,-b,-c
489 u1 = (-b+ctx.sqrt(b**2-4*a*c))/(2*c)
490 u2 = (-b-ctx.sqrt(b**2-4*a*c))/(2*c)
491 if abs(u1-t) < abs(u2-t):
492 if b: s = '((%s+sqrt(%s))/%s)' % (-b,b**2-4*a*c,2*c)
493 else: s = '(sqrt(%s)/%s)' % (-4*a*c,2*c)
494 else:
495 if b: s = '((%s-sqrt(%s))/%s)' % (-b,b**2-4*a*c,2*c)
496 else: s = '(-sqrt(%s)/%s)' % (-4*a*c,2*c)
497 return s
499# Transformation y = f(x,c), with inverse function x = f(y,c)
500# The third entry indicates whether the transformation is
501# redundant when c = 1
502transforms = [
503 (lambda ctx,x,c: x*c, '$y/$c', 0),
504 (lambda ctx,x,c: x/c, '$c*$y', 1),
505 (lambda ctx,x,c: c/x, '$c/$y', 0),
506 (lambda ctx,x,c: (x*c)**2, 'sqrt($y)/$c', 0),
507 (lambda ctx,x,c: (x/c)**2, '$c*sqrt($y)', 1),
508 (lambda ctx,x,c: (c/x)**2, '$c/sqrt($y)', 0),
509 (lambda ctx,x,c: c*x**2, 'sqrt($y)/sqrt($c)', 1),
510 (lambda ctx,x,c: x**2/c, 'sqrt($c)*sqrt($y)', 1),
511 (lambda ctx,x,c: c/x**2, 'sqrt($c)/sqrt($y)', 1),
512 (lambda ctx,x,c: ctx.sqrt(x*c), '$y**2/$c', 0),
513 (lambda ctx,x,c: ctx.sqrt(x/c), '$c*$y**2', 1),
514 (lambda ctx,x,c: ctx.sqrt(c/x), '$c/$y**2', 0),
515 (lambda ctx,x,c: c*ctx.sqrt(x), '$y**2/$c**2', 1),
516 (lambda ctx,x,c: ctx.sqrt(x)/c, '$c**2*$y**2', 1),
517 (lambda ctx,x,c: c/ctx.sqrt(x), '$c**2/$y**2', 1),
518 (lambda ctx,x,c: ctx.exp(x*c), 'log($y)/$c', 0),
519 (lambda ctx,x,c: ctx.exp(x/c), '$c*log($y)', 1),
520 (lambda ctx,x,c: ctx.exp(c/x), '$c/log($y)', 0),
521 (lambda ctx,x,c: c*ctx.exp(x), 'log($y/$c)', 1),
522 (lambda ctx,x,c: ctx.exp(x)/c, 'log($c*$y)', 1),
523 (lambda ctx,x,c: c/ctx.exp(x), 'log($c/$y)', 0),
524 (lambda ctx,x,c: ctx.ln(x*c), 'exp($y)/$c', 0),
525 (lambda ctx,x,c: ctx.ln(x/c), '$c*exp($y)', 1),
526 (lambda ctx,x,c: ctx.ln(c/x), '$c/exp($y)', 0),
527 (lambda ctx,x,c: c*ctx.ln(x), 'exp($y/$c)', 1),
528 (lambda ctx,x,c: ctx.ln(x)/c, 'exp($c*$y)', 1),
529 (lambda ctx,x,c: c/ctx.ln(x), 'exp($c/$y)', 0),
530]
532def identify(ctx, x, constants=[], tol=None, maxcoeff=1000, full=False,
533 verbose=False):
534 r"""
535 Given a real number `x`, ``identify(x)`` attempts to find an exact
536 formula for `x`. This formula is returned as a string. If no match
537 is found, ``None`` is returned. With ``full=True``, a list of
538 matching formulas is returned.
540 As a simple example, :func:`~mpmath.identify` will find an algebraic
541 formula for the golden ratio::
543 >>> from mpmath import *
544 >>> mp.dps = 15; mp.pretty = True
545 >>> identify(phi)
546 '((1+sqrt(5))/2)'
548 :func:`~mpmath.identify` can identify simple algebraic numbers and simple
549 combinations of given base constants, as well as certain basic
550 transformations thereof. More specifically, :func:`~mpmath.identify`
551 looks for the following:
553 1. Fractions
554 2. Quadratic algebraic numbers
555 3. Rational linear combinations of the base constants
556 4. Any of the above after first transforming `x` into `f(x)` where
557 `f(x)` is `1/x`, `\sqrt x`, `x^2`, `\log x` or `\exp x`, either
558 directly or with `x` or `f(x)` multiplied or divided by one of
559 the base constants
560 5. Products of fractional powers of the base constants and
561 small integers
563 Base constants can be given as a list of strings representing mpmath
564 expressions (:func:`~mpmath.identify` will ``eval`` the strings to numerical
565 values and use the original strings for the output), or as a dict of
566 formula:value pairs.
568 In order not to produce spurious results, :func:`~mpmath.identify` should
569 be used with high precision; preferably 50 digits or more.
571 **Examples**
573 Simple identifications can be performed safely at standard
574 precision. Here the default recognition of rational, algebraic,
575 and exp/log of algebraic numbers is demonstrated::
577 >>> mp.dps = 15
578 >>> identify(0.22222222222222222)
579 '(2/9)'
580 >>> identify(1.9662210973805663)
581 'sqrt(((24+sqrt(48))/8))'
582 >>> identify(4.1132503787829275)
583 'exp((sqrt(8)/2))'
584 >>> identify(0.881373587019543)
585 'log(((2+sqrt(8))/2))'
587 By default, :func:`~mpmath.identify` does not recognize `\pi`. At standard
588 precision it finds a not too useful approximation. At slightly
589 increased precision, this approximation is no longer accurate
590 enough and :func:`~mpmath.identify` more correctly returns ``None``::
592 >>> identify(pi)
593 '(2**(176/117)*3**(20/117)*5**(35/39))/(7**(92/117))'
594 >>> mp.dps = 30
595 >>> identify(pi)
596 >>>
598 Numbers such as `\pi`, and simple combinations of user-defined
599 constants, can be identified if they are provided explicitly::
601 >>> identify(3*pi-2*e, ['pi', 'e'])
602 '(3*pi + (-2)*e)'
604 Here is an example using a dict of constants. Note that the
605 constants need not be "atomic"; :func:`~mpmath.identify` can just
606 as well express the given number in terms of expressions
607 given by formulas::
609 >>> identify(pi+e, {'a':pi+2, 'b':2*e})
610 '((-2) + 1*a + (1/2)*b)'
612 Next, we attempt some identifications with a set of base constants.
613 It is necessary to increase the precision a bit.
615 >>> mp.dps = 50
616 >>> base = ['sqrt(2)','pi','log(2)']
617 >>> identify(0.25, base)
618 '(1/4)'
619 >>> identify(3*pi + 2*sqrt(2) + 5*log(2)/7, base)
620 '(2*sqrt(2) + 3*pi + (5/7)*log(2))'
621 >>> identify(exp(pi+2), base)
622 'exp((2 + 1*pi))'
623 >>> identify(1/(3+sqrt(2)), base)
624 '((3/7) + (-1/7)*sqrt(2))'
625 >>> identify(sqrt(2)/(3*pi+4), base)
626 'sqrt(2)/(4 + 3*pi)'
627 >>> identify(5**(mpf(1)/3)*pi*log(2)**2, base)
628 '5**(1/3)*pi*log(2)**2'
630 An example of an erroneous solution being found when too low
631 precision is used::
633 >>> mp.dps = 15
634 >>> identify(1/(3*pi-4*e+sqrt(8)), ['pi', 'e', 'sqrt(2)'])
635 '((11/25) + (-158/75)*pi + (76/75)*e + (44/15)*sqrt(2))'
636 >>> mp.dps = 50
637 >>> identify(1/(3*pi-4*e+sqrt(8)), ['pi', 'e', 'sqrt(2)'])
638 '1/(3*pi + (-4)*e + 2*sqrt(2))'
640 **Finding approximate solutions**
642 The tolerance ``tol`` defaults to 3/4 of the working precision.
643 Lowering the tolerance is useful for finding approximate matches.
644 We can for example try to generate approximations for pi::
646 >>> mp.dps = 15
647 >>> identify(pi, tol=1e-2)
648 '(22/7)'
649 >>> identify(pi, tol=1e-3)
650 '(355/113)'
651 >>> identify(pi, tol=1e-10)
652 '(5**(339/269))/(2**(64/269)*3**(13/269)*7**(92/269))'
654 With ``full=True``, and by supplying a few base constants,
655 ``identify`` can generate almost endless lists of approximations
656 for any number (the output below has been truncated to show only
657 the first few)::
659 >>> for p in identify(pi, ['e', 'catalan'], tol=1e-5, full=True):
660 ... print(p)
661 ... # doctest: +ELLIPSIS
662 e/log((6 + (-4/3)*e))
663 (3**3*5*e*catalan**2)/(2*7**2)
664 sqrt(((-13) + 1*e + 22*catalan))
665 log(((-6) + 24*e + 4*catalan)/e)
666 exp(catalan*((-1/5) + (8/15)*e))
667 catalan*(6 + (-6)*e + 15*catalan)
668 sqrt((5 + 26*e + (-3)*catalan))/e
669 e*sqrt(((-27) + 2*e + 25*catalan))
670 log(((-1) + (-11)*e + 59*catalan))
671 ((3/20) + (21/20)*e + (3/20)*catalan)
672 ...
674 The numerical values are roughly as close to `\pi` as permitted by the
675 specified tolerance:
677 >>> e/log(6-4*e/3)
678 3.14157719846001
679 >>> 135*e*catalan**2/98
680 3.14166950419369
681 >>> sqrt(e-13+22*catalan)
682 3.14158000062992
683 >>> log(24*e-6+4*catalan)-1
684 3.14158791577159
686 **Symbolic processing**
688 The output formula can be evaluated as a Python expression.
689 Note however that if fractions (like '2/3') are present in
690 the formula, Python's :func:`~mpmath.eval()` may erroneously perform
691 integer division. Note also that the output is not necessarily
692 in the algebraically simplest form::
694 >>> identify(sqrt(2))
695 '(sqrt(8)/2)'
697 As a solution to both problems, consider using SymPy's
698 :func:`~mpmath.sympify` to convert the formula into a symbolic expression.
699 SymPy can be used to pretty-print or further simplify the formula
700 symbolically::
702 >>> from sympy import sympify # doctest: +SKIP
703 >>> sympify(identify(sqrt(2))) # doctest: +SKIP
704 2**(1/2)
706 Sometimes :func:`~mpmath.identify` can simplify an expression further than
707 a symbolic algorithm::
709 >>> from sympy import simplify # doctest: +SKIP
710 >>> x = sympify('-1/(-3/2+(1/2)*5**(1/2))*(3/2-1/2*5**(1/2))**(1/2)') # doctest: +SKIP
711 >>> x # doctest: +SKIP
712 (3/2 - 5**(1/2)/2)**(-1/2)
713 >>> x = simplify(x) # doctest: +SKIP
714 >>> x # doctest: +SKIP
715 2/(6 - 2*5**(1/2))**(1/2)
716 >>> mp.dps = 30 # doctest: +SKIP
717 >>> x = sympify(identify(x.evalf(30))) # doctest: +SKIP
718 >>> x # doctest: +SKIP
719 1/2 + 5**(1/2)/2
721 (In fact, this functionality is available directly in SymPy as the
722 function :func:`~mpmath.nsimplify`, which is essentially a wrapper for
723 :func:`~mpmath.identify`.)
725 **Miscellaneous issues and limitations**
727 The input `x` must be a real number. All base constants must be
728 positive real numbers and must not be rationals or rational linear
729 combinations of each other.
731 The worst-case computation time grows quickly with the number of
732 base constants. Already with 3 or 4 base constants,
733 :func:`~mpmath.identify` may require several seconds to finish. To search
734 for relations among a large number of constants, you should
735 consider using :func:`~mpmath.pslq` directly.
737 The extended transformations are applied to x, not the constants
738 separately. As a result, ``identify`` will for example be able to
739 recognize ``exp(2*pi+3)`` with ``pi`` given as a base constant, but
740 not ``2*exp(pi)+3``. It will be able to recognize the latter if
741 ``exp(pi)`` is given explicitly as a base constant.
743 """
745 solutions = []
747 def addsolution(s):
748 if verbose: print("Found: ", s)
749 solutions.append(s)
751 x = ctx.mpf(x)
753 # Further along, x will be assumed positive
754 if x == 0:
755 if full: return ['0']
756 else: return '0'
757 if x < 0:
758 sol = ctx.identify(-x, constants, tol, maxcoeff, full, verbose)
759 if sol is None:
760 return sol
761 if full:
762 return ["-(%s)"%s for s in sol]
763 else:
764 return "-(%s)" % sol
766 if tol:
767 tol = ctx.mpf(tol)
768 else:
769 tol = ctx.eps**0.7
770 M = maxcoeff
772 if constants:
773 if isinstance(constants, dict):
774 constants = [(ctx.mpf(v), name) for (name, v) in sorted(constants.items())]
775 else:
776 namespace = dict((name, getattr(ctx,name)) for name in dir(ctx))
777 constants = [(eval(p, namespace), p) for p in constants]
778 else:
779 constants = []
781 # We always want to find at least rational terms
782 if 1 not in [value for (name, value) in constants]:
783 constants = [(ctx.mpf(1), '1')] + constants
785 # PSLQ with simple algebraic and functional transformations
786 for ft, ftn, red in transforms:
787 for c, cn in constants:
788 if red and cn == '1':
789 continue
790 t = ft(ctx,x,c)
791 # Prevent exponential transforms from wreaking havoc
792 if abs(t) > M**2 or abs(t) < tol:
793 continue
794 # Linear combination of base constants
795 r = ctx.pslq([t] + [a[0] for a in constants], tol, M)
796 s = None
797 if r is not None and max(abs(uw) for uw in r) <= M and r[0]:
798 s = pslqstring(r, constants)
799 # Quadratic algebraic numbers
800 else:
801 q = ctx.pslq([ctx.one, t, t**2], tol, M)
802 if q is not None and len(q) == 3 and q[2]:
803 aa, bb, cc = q
804 if max(abs(aa),abs(bb),abs(cc)) <= M:
805 s = quadraticstring(ctx,t,aa,bb,cc)
806 if s:
807 if cn == '1' and ('/$c' in ftn):
808 s = ftn.replace('$y', s).replace('/$c', '')
809 else:
810 s = ftn.replace('$y', s).replace('$c', cn)
811 addsolution(s)
812 if not full: return solutions[0]
814 if verbose:
815 print(".")
817 # Check for a direct multiplicative formula
818 if x != 1:
819 # Allow fractional powers of fractions
820 ilogs = [2,3,5,7]
821 # Watch out for existing fractional powers of fractions
822 logs = []
823 for a, s in constants:
824 if not sum(bool(ctx.findpoly(ctx.ln(a)/ctx.ln(i),1)) for i in ilogs):
825 logs.append((ctx.ln(a), s))
826 logs = [(ctx.ln(i),str(i)) for i in ilogs] + logs
827 r = ctx.pslq([ctx.ln(x)] + [a[0] for a in logs], tol, M)
828 if r is not None and max(abs(uw) for uw in r) <= M and r[0]:
829 addsolution(prodstring(r, logs))
830 if not full: return solutions[0]
832 if full:
833 return sorted(solutions, key=len)
834 else:
835 return None
837IdentificationMethods.pslq = pslq
838IdentificationMethods.findpoly = findpoly
839IdentificationMethods.identify = identify
842if __name__ == '__main__':
843 import doctest
844 doctest.testmod()