Coverage for /usr/lib/python3/dist-packages/sympy/solvers/diophantine/diophantine.py: 10%
1527 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core.add import Add
2from sympy.core.assumptions import check_assumptions
3from sympy.core.containers import Tuple
4from sympy.core.exprtools import factor_terms
5from sympy.core.function import _mexpand
6from sympy.core.mul import Mul
7from sympy.core.numbers import Rational
8from sympy.core.numbers import igcdex, ilcm, igcd
9from sympy.core.power import integer_nthroot, isqrt
10from sympy.core.relational import Eq
11from sympy.core.singleton import S
12from sympy.core.sorting import default_sort_key, ordered
13from sympy.core.symbol import Symbol, symbols
14from sympy.core.sympify import _sympify
15from sympy.functions.elementary.complexes import sign
16from sympy.functions.elementary.integers import floor
17from sympy.functions.elementary.miscellaneous import sqrt
18from sympy.matrices.dense import MutableDenseMatrix as Matrix
19from sympy.ntheory.factor_ import (
20 divisors, factorint, multiplicity, perfect_power)
21from sympy.ntheory.generate import nextprime
22from sympy.ntheory.primetest import is_square, isprime
23from sympy.ntheory.residue_ntheory import sqrt_mod
24from sympy.polys.polyerrors import GeneratorsNeeded
25from sympy.polys.polytools import Poly, factor_list
26from sympy.simplify.simplify import signsimp
27from sympy.solvers.solveset import solveset_real
28from sympy.utilities import numbered_symbols
29from sympy.utilities.misc import as_int, filldedent
30from sympy.utilities.iterables import (is_sequence, subsets, permute_signs,
31 signed_permutations, ordered_partitions)
34# these are imported with 'from sympy.solvers.diophantine import *
35__all__ = ['diophantine', 'classify_diop']
38class DiophantineSolutionSet(set):
39 """
40 Container for a set of solutions to a particular diophantine equation.
42 The base representation is a set of tuples representing each of the solutions.
44 Parameters
45 ==========
47 symbols : list
48 List of free symbols in the original equation.
49 parameters: list
50 List of parameters to be used in the solution.
52 Examples
53 ========
55 Adding solutions:
57 >>> from sympy.solvers.diophantine.diophantine import DiophantineSolutionSet
58 >>> from sympy.abc import x, y, t, u
59 >>> s1 = DiophantineSolutionSet([x, y], [t, u])
60 >>> s1
61 set()
62 >>> s1.add((2, 3))
63 >>> s1.add((-1, u))
64 >>> s1
65 {(-1, u), (2, 3)}
66 >>> s2 = DiophantineSolutionSet([x, y], [t, u])
67 >>> s2.add((3, 4))
68 >>> s1.update(*s2)
69 >>> s1
70 {(-1, u), (2, 3), (3, 4)}
72 Conversion of solutions into dicts:
74 >>> list(s1.dict_iterator())
75 [{x: -1, y: u}, {x: 2, y: 3}, {x: 3, y: 4}]
77 Substituting values:
79 >>> s3 = DiophantineSolutionSet([x, y], [t, u])
80 >>> s3.add((t**2, t + u))
81 >>> s3
82 {(t**2, t + u)}
83 >>> s3.subs({t: 2, u: 3})
84 {(4, 5)}
85 >>> s3.subs(t, -1)
86 {(1, u - 1)}
87 >>> s3.subs(t, 3)
88 {(9, u + 3)}
90 Evaluation at specific values. Positional arguments are given in the same order as the parameters:
92 >>> s3(-2, 3)
93 {(4, 1)}
94 >>> s3(5)
95 {(25, u + 5)}
96 >>> s3(None, 2)
97 {(t**2, t + 2)}
98 """
100 def __init__(self, symbols_seq, parameters):
101 super().__init__()
103 if not is_sequence(symbols_seq):
104 raise ValueError("Symbols must be given as a sequence.")
106 if not is_sequence(parameters):
107 raise ValueError("Parameters must be given as a sequence.")
109 self.symbols = tuple(symbols_seq)
110 self.parameters = tuple(parameters)
112 def add(self, solution):
113 if len(solution) != len(self.symbols):
114 raise ValueError("Solution should have a length of %s, not %s" % (len(self.symbols), len(solution)))
115 super().add(Tuple(*solution))
117 def update(self, *solutions):
118 for solution in solutions:
119 self.add(solution)
121 def dict_iterator(self):
122 for solution in ordered(self):
123 yield dict(zip(self.symbols, solution))
125 def subs(self, *args, **kwargs):
126 result = DiophantineSolutionSet(self.symbols, self.parameters)
127 for solution in self:
128 result.add(solution.subs(*args, **kwargs))
129 return result
131 def __call__(self, *args):
132 if len(args) > len(self.parameters):
133 raise ValueError("Evaluation should have at most %s values, not %s" % (len(self.parameters), len(args)))
134 rep = {p: v for p, v in zip(self.parameters, args) if v is not None}
135 return self.subs(rep)
138class DiophantineEquationType:
139 """
140 Internal representation of a particular diophantine equation type.
142 Parameters
143 ==========
145 equation :
146 The diophantine equation that is being solved.
147 free_symbols : list (optional)
148 The symbols being solved for.
150 Attributes
151 ==========
153 total_degree :
154 The maximum of the degrees of all terms in the equation
155 homogeneous :
156 Does the equation contain a term of degree 0
157 homogeneous_order :
158 Does the equation contain any coefficient that is in the symbols being solved for
159 dimension :
160 The number of symbols being solved for
161 """
162 name = None # type: str
164 def __init__(self, equation, free_symbols=None):
165 self.equation = _sympify(equation).expand(force=True)
167 if free_symbols is not None:
168 self.free_symbols = free_symbols
169 else:
170 self.free_symbols = list(self.equation.free_symbols)
171 self.free_symbols.sort(key=default_sort_key)
173 if not self.free_symbols:
174 raise ValueError('equation should have 1 or more free symbols')
176 self.coeff = self.equation.as_coefficients_dict()
177 if not all(_is_int(c) for c in self.coeff.values()):
178 raise TypeError("Coefficients should be Integers")
180 self.total_degree = Poly(self.equation).total_degree()
181 self.homogeneous = 1 not in self.coeff
182 self.homogeneous_order = not (set(self.coeff) & set(self.free_symbols))
183 self.dimension = len(self.free_symbols)
184 self._parameters = None
186 def matches(self):
187 """
188 Determine whether the given equation can be matched to the particular equation type.
189 """
190 return False
192 @property
193 def n_parameters(self):
194 return self.dimension
196 @property
197 def parameters(self):
198 if self._parameters is None:
199 self._parameters = symbols('t_:%i' % (self.n_parameters,), integer=True)
200 return self._parameters
202 def solve(self, parameters=None, limit=None) -> DiophantineSolutionSet:
203 raise NotImplementedError('No solver has been written for %s.' % self.name)
205 def pre_solve(self, parameters=None):
206 if not self.matches():
207 raise ValueError("This equation does not match the %s equation type." % self.name)
209 if parameters is not None:
210 if len(parameters) != self.n_parameters:
211 raise ValueError("Expected %s parameter(s) but got %s" % (self.n_parameters, len(parameters)))
213 self._parameters = parameters
216class Univariate(DiophantineEquationType):
217 """
218 Representation of a univariate diophantine equation.
220 A univariate diophantine equation is an equation of the form
221 `a_{0} + a_{1}x + a_{2}x^2 + .. + a_{n}x^n = 0` where `a_{1}, a_{2}, ..a_{n}` are
222 integer constants and `x` is an integer variable.
224 Examples
225 ========
227 >>> from sympy.solvers.diophantine.diophantine import Univariate
228 >>> from sympy.abc import x
229 >>> Univariate((x - 2)*(x - 3)**2).solve() # solves equation (x - 2)*(x - 3)**2 == 0
230 {(2,), (3,)}
232 """
234 name = 'univariate'
236 def matches(self):
237 return self.dimension == 1
239 def solve(self, parameters=None, limit=None):
240 self.pre_solve(parameters)
242 result = DiophantineSolutionSet(self.free_symbols, parameters=self.parameters)
243 for i in solveset_real(self.equation, self.free_symbols[0]).intersect(S.Integers):
244 result.add((i,))
245 return result
248class Linear(DiophantineEquationType):
249 """
250 Representation of a linear diophantine equation.
252 A linear diophantine equation is an equation of the form `a_{1}x_{1} +
253 a_{2}x_{2} + .. + a_{n}x_{n} = 0` where `a_{1}, a_{2}, ..a_{n}` are
254 integer constants and `x_{1}, x_{2}, ..x_{n}` are integer variables.
256 Examples
257 ========
259 >>> from sympy.solvers.diophantine.diophantine import Linear
260 >>> from sympy.abc import x, y, z
261 >>> l1 = Linear(2*x - 3*y - 5)
262 >>> l1.matches() # is this equation linear
263 True
264 >>> l1.solve() # solves equation 2*x - 3*y - 5 == 0
265 {(3*t_0 - 5, 2*t_0 - 5)}
267 Here x = -3*t_0 - 5 and y = -2*t_0 - 5
269 >>> Linear(2*x - 3*y - 4*z -3).solve()
270 {(t_0, 2*t_0 + 4*t_1 + 3, -t_0 - 3*t_1 - 3)}
272 """
274 name = 'linear'
276 def matches(self):
277 return self.total_degree == 1
279 def solve(self, parameters=None, limit=None):
280 self.pre_solve(parameters)
282 coeff = self.coeff
283 var = self.free_symbols
285 if 1 in coeff:
286 # negate coeff[] because input is of the form: ax + by + c == 0
287 # but is used as: ax + by == -c
288 c = -coeff[1]
289 else:
290 c = 0
292 result = DiophantineSolutionSet(var, parameters=self.parameters)
293 params = result.parameters
295 if len(var) == 1:
296 q, r = divmod(c, coeff[var[0]])
297 if not r:
298 result.add((q,))
299 return result
300 else:
301 return result
303 '''
304 base_solution_linear() can solve diophantine equations of the form:
306 a*x + b*y == c
308 We break down multivariate linear diophantine equations into a
309 series of bivariate linear diophantine equations which can then
310 be solved individually by base_solution_linear().
312 Consider the following:
314 a_0*x_0 + a_1*x_1 + a_2*x_2 == c
316 which can be re-written as:
318 a_0*x_0 + g_0*y_0 == c
320 where
322 g_0 == gcd(a_1, a_2)
324 and
326 y == (a_1*x_1)/g_0 + (a_2*x_2)/g_0
328 This leaves us with two binary linear diophantine equations.
329 For the first equation:
331 a == a_0
332 b == g_0
333 c == c
335 For the second:
337 a == a_1/g_0
338 b == a_2/g_0
339 c == the solution we find for y_0 in the first equation.
341 The arrays A and B are the arrays of integers used for
342 'a' and 'b' in each of the n-1 bivariate equations we solve.
343 '''
345 A = [coeff[v] for v in var]
346 B = []
347 if len(var) > 2:
348 B.append(igcd(A[-2], A[-1]))
349 A[-2] = A[-2] // B[0]
350 A[-1] = A[-1] // B[0]
351 for i in range(len(A) - 3, 0, -1):
352 gcd = igcd(B[0], A[i])
353 B[0] = B[0] // gcd
354 A[i] = A[i] // gcd
355 B.insert(0, gcd)
356 B.append(A[-1])
358 '''
359 Consider the trivariate linear equation:
361 4*x_0 + 6*x_1 + 3*x_2 == 2
363 This can be re-written as:
365 4*x_0 + 3*y_0 == 2
367 where
369 y_0 == 2*x_1 + x_2
370 (Note that gcd(3, 6) == 3)
372 The complete integral solution to this equation is:
374 x_0 == 2 + 3*t_0
375 y_0 == -2 - 4*t_0
377 where 't_0' is any integer.
379 Now that we have a solution for 'x_0', find 'x_1' and 'x_2':
381 2*x_1 + x_2 == -2 - 4*t_0
383 We can then solve for '-2' and '-4' independently,
384 and combine the results:
386 2*x_1a + x_2a == -2
387 x_1a == 0 + t_0
388 x_2a == -2 - 2*t_0
390 2*x_1b + x_2b == -4*t_0
391 x_1b == 0*t_0 + t_1
392 x_2b == -4*t_0 - 2*t_1
394 ==>
396 x_1 == t_0 + t_1
397 x_2 == -2 - 6*t_0 - 2*t_1
399 where 't_0' and 't_1' are any integers.
401 Note that:
403 4*(2 + 3*t_0) + 6*(t_0 + t_1) + 3*(-2 - 6*t_0 - 2*t_1) == 2
405 for any integral values of 't_0', 't_1'; as required.
407 This method is generalised for many variables, below.
409 '''
410 solutions = []
411 for Ai, Bi in zip(A, B):
412 tot_x, tot_y = [], []
414 for j, arg in enumerate(Add.make_args(c)):
415 if arg.is_Integer:
416 # example: 5 -> k = 5
417 k, p = arg, S.One
418 pnew = params[0]
419 else: # arg is a Mul or Symbol
420 # example: 3*t_1 -> k = 3
421 # example: t_0 -> k = 1
422 k, p = arg.as_coeff_Mul()
423 pnew = params[params.index(p) + 1]
425 sol = sol_x, sol_y = base_solution_linear(k, Ai, Bi, pnew)
427 if p is S.One:
428 if None in sol:
429 return result
430 else:
431 # convert a + b*pnew -> a*p + b*pnew
432 if isinstance(sol_x, Add):
433 sol_x = sol_x.args[0]*p + sol_x.args[1]
434 if isinstance(sol_y, Add):
435 sol_y = sol_y.args[0]*p + sol_y.args[1]
437 tot_x.append(sol_x)
438 tot_y.append(sol_y)
440 solutions.append(Add(*tot_x))
441 c = Add(*tot_y)
443 solutions.append(c)
444 result.add(solutions)
445 return result
448class BinaryQuadratic(DiophantineEquationType):
449 """
450 Representation of a binary quadratic diophantine equation.
452 A binary quadratic diophantine equation is an equation of the
453 form `Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0`, where `A, B, C, D, E,
454 F` are integer constants and `x` and `y` are integer variables.
456 Examples
457 ========
459 >>> from sympy.abc import x, y
460 >>> from sympy.solvers.diophantine.diophantine import BinaryQuadratic
461 >>> b1 = BinaryQuadratic(x**3 + y**2 + 1)
462 >>> b1.matches()
463 False
464 >>> b2 = BinaryQuadratic(x**2 + y**2 + 2*x + 2*y + 2)
465 >>> b2.matches()
466 True
467 >>> b2.solve()
468 {(-1, -1)}
470 References
471 ==========
473 .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
474 Available: https://www.alpertron.com.ar/METHODS.HTM
475 .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
476 Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
478 """
480 name = 'binary_quadratic'
482 def matches(self):
483 return self.total_degree == 2 and self.dimension == 2
485 def solve(self, parameters=None, limit=None) -> DiophantineSolutionSet:
486 self.pre_solve(parameters)
488 var = self.free_symbols
489 coeff = self.coeff
491 x, y = var
493 A = coeff[x**2]
494 B = coeff[x*y]
495 C = coeff[y**2]
496 D = coeff[x]
497 E = coeff[y]
498 F = coeff[S.One]
500 A, B, C, D, E, F = [as_int(i) for i in _remove_gcd(A, B, C, D, E, F)]
502 # (1) Simple-Hyperbolic case: A = C = 0, B != 0
503 # In this case equation can be converted to (Bx + E)(By + D) = DE - BF
504 # We consider two cases; DE - BF = 0 and DE - BF != 0
505 # More details, https://www.alpertron.com.ar/METHODS.HTM#SHyperb
507 result = DiophantineSolutionSet(var, self.parameters)
508 t, u = result.parameters
510 discr = B**2 - 4*A*C
511 if A == 0 and C == 0 and B != 0:
513 if D*E - B*F == 0:
514 q, r = divmod(E, B)
515 if not r:
516 result.add((-q, t))
517 q, r = divmod(D, B)
518 if not r:
519 result.add((t, -q))
520 else:
521 div = divisors(D*E - B*F)
522 div = div + [-term for term in div]
523 for d in div:
524 x0, r = divmod(d - E, B)
525 if not r:
526 q, r = divmod(D*E - B*F, d)
527 if not r:
528 y0, r = divmod(q - D, B)
529 if not r:
530 result.add((x0, y0))
532 # (2) Parabolic case: B**2 - 4*A*C = 0
533 # There are two subcases to be considered in this case.
534 # sqrt(c)D - sqrt(a)E = 0 and sqrt(c)D - sqrt(a)E != 0
535 # More Details, https://www.alpertron.com.ar/METHODS.HTM#Parabol
537 elif discr == 0:
539 if A == 0:
540 s = BinaryQuadratic(self.equation, free_symbols=[y, x]).solve(parameters=[t, u])
541 for soln in s:
542 result.add((soln[1], soln[0]))
544 else:
545 g = sign(A)*igcd(A, C)
546 a = A // g
547 c = C // g
548 e = sign(B / A)
550 sqa = isqrt(a)
551 sqc = isqrt(c)
552 _c = e*sqc*D - sqa*E
553 if not _c:
554 z = Symbol("z", real=True)
555 eq = sqa*g*z**2 + D*z + sqa*F
556 roots = solveset_real(eq, z).intersect(S.Integers)
557 for root in roots:
558 ans = diop_solve(sqa*x + e*sqc*y - root)
559 result.add((ans[0], ans[1]))
561 elif _is_int(c):
562 solve_x = lambda u: -e*sqc*g*_c*t**2 - (E + 2*e*sqc*g*u)*t \
563 - (e*sqc*g*u**2 + E*u + e*sqc*F) // _c
565 solve_y = lambda u: sqa*g*_c*t**2 + (D + 2*sqa*g*u)*t \
566 + (sqa*g*u**2 + D*u + sqa*F) // _c
568 for z0 in range(0, abs(_c)):
569 # Check if the coefficients of y and x obtained are integers or not
570 if (divisible(sqa*g*z0**2 + D*z0 + sqa*F, _c) and
571 divisible(e*sqc*g*z0**2 + E*z0 + e*sqc*F, _c)):
572 result.add((solve_x(z0), solve_y(z0)))
574 # (3) Method used when B**2 - 4*A*C is a square, is described in p. 6 of the below paper
575 # by John P. Robertson.
576 # https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
578 elif is_square(discr):
579 if A != 0:
580 r = sqrt(discr)
581 u, v = symbols("u, v", integer=True)
582 eq = _mexpand(
583 4*A*r*u*v + 4*A*D*(B*v + r*u + r*v - B*u) +
584 2*A*4*A*E*(u - v) + 4*A*r*4*A*F)
586 solution = diop_solve(eq, t)
588 for s0, t0 in solution:
590 num = B*t0 + r*s0 + r*t0 - B*s0
591 x_0 = S(num) / (4*A*r)
592 y_0 = S(s0 - t0) / (2*r)
593 if isinstance(s0, Symbol) or isinstance(t0, Symbol):
594 if len(check_param(x_0, y_0, 4*A*r, parameters)) > 0:
595 ans = check_param(x_0, y_0, 4*A*r, parameters)
596 result.update(*ans)
597 elif x_0.is_Integer and y_0.is_Integer:
598 if is_solution_quad(var, coeff, x_0, y_0):
599 result.add((x_0, y_0))
601 else:
602 s = BinaryQuadratic(self.equation, free_symbols=var[::-1]).solve(parameters=[t, u]) # Interchange x and y
603 while s:
604 result.add(s.pop()[::-1]) # and solution <--------+
606 # (4) B**2 - 4*A*C > 0 and B**2 - 4*A*C not a square or B**2 - 4*A*C < 0
608 else:
610 P, Q = _transformation_to_DN(var, coeff)
611 D, N = _find_DN(var, coeff)
612 solns_pell = diop_DN(D, N)
614 if D < 0:
615 for x0, y0 in solns_pell:
616 for x in [-x0, x0]:
617 for y in [-y0, y0]:
618 s = P*Matrix([x, y]) + Q
619 try:
620 result.add([as_int(_) for _ in s])
621 except ValueError:
622 pass
623 else:
624 # In this case equation can be transformed into a Pell equation
626 solns_pell = set(solns_pell)
627 for X, Y in list(solns_pell):
628 solns_pell.add((-X, -Y))
630 a = diop_DN(D, 1)
631 T = a[0][0]
632 U = a[0][1]
634 if all(_is_int(_) for _ in P[:4] + Q[:2]):
635 for r, s in solns_pell:
636 _a = (r + s*sqrt(D))*(T + U*sqrt(D))**t
637 _b = (r - s*sqrt(D))*(T - U*sqrt(D))**t
638 x_n = _mexpand(S(_a + _b) / 2)
639 y_n = _mexpand(S(_a - _b) / (2*sqrt(D)))
640 s = P*Matrix([x_n, y_n]) + Q
641 result.add(s)
643 else:
644 L = ilcm(*[_.q for _ in P[:4] + Q[:2]])
646 k = 1
648 T_k = T
649 U_k = U
651 while (T_k - 1) % L != 0 or U_k % L != 0:
652 T_k, U_k = T_k*T + D*U_k*U, T_k*U + U_k*T
653 k += 1
655 for X, Y in solns_pell:
657 for i in range(k):
658 if all(_is_int(_) for _ in P*Matrix([X, Y]) + Q):
659 _a = (X + sqrt(D)*Y)*(T_k + sqrt(D)*U_k)**t
660 _b = (X - sqrt(D)*Y)*(T_k - sqrt(D)*U_k)**t
661 Xt = S(_a + _b) / 2
662 Yt = S(_a - _b) / (2*sqrt(D))
663 s = P*Matrix([Xt, Yt]) + Q
664 result.add(s)
666 X, Y = X*T + D*U*Y, X*U + Y*T
668 return result
671class InhomogeneousTernaryQuadratic(DiophantineEquationType):
672 """
674 Representation of an inhomogeneous ternary quadratic.
676 No solver is currently implemented for this equation type.
678 """
680 name = 'inhomogeneous_ternary_quadratic'
682 def matches(self):
683 if not (self.total_degree == 2 and self.dimension == 3):
684 return False
685 if not self.homogeneous:
686 return False
687 return not self.homogeneous_order
690class HomogeneousTernaryQuadraticNormal(DiophantineEquationType):
691 """
692 Representation of a homogeneous ternary quadratic normal diophantine equation.
694 Examples
695 ========
697 >>> from sympy.abc import x, y, z
698 >>> from sympy.solvers.diophantine.diophantine import HomogeneousTernaryQuadraticNormal
699 >>> HomogeneousTernaryQuadraticNormal(4*x**2 - 5*y**2 + z**2).solve()
700 {(1, 2, 4)}
702 """
704 name = 'homogeneous_ternary_quadratic_normal'
706 def matches(self):
707 if not (self.total_degree == 2 and self.dimension == 3):
708 return False
709 if not self.homogeneous:
710 return False
711 if not self.homogeneous_order:
712 return False
714 nonzero = [k for k in self.coeff if self.coeff[k]]
715 return len(nonzero) == 3 and all(i**2 in nonzero for i in self.free_symbols)
717 def solve(self, parameters=None, limit=None) -> DiophantineSolutionSet:
718 self.pre_solve(parameters)
720 var = self.free_symbols
721 coeff = self.coeff
723 x, y, z = var
725 a = coeff[x**2]
726 b = coeff[y**2]
727 c = coeff[z**2]
729 (sqf_of_a, sqf_of_b, sqf_of_c), (a_1, b_1, c_1), (a_2, b_2, c_2) = \
730 sqf_normal(a, b, c, steps=True)
732 A = -a_2*c_2
733 B = -b_2*c_2
735 result = DiophantineSolutionSet(var, parameters=self.parameters)
737 # If following two conditions are satisfied then there are no solutions
738 if A < 0 and B < 0:
739 return result
741 if (
742 sqrt_mod(-b_2*c_2, a_2) is None or
743 sqrt_mod(-c_2*a_2, b_2) is None or
744 sqrt_mod(-a_2*b_2, c_2) is None):
745 return result
747 z_0, x_0, y_0 = descent(A, B)
749 z_0, q = _rational_pq(z_0, abs(c_2))
750 x_0 *= q
751 y_0 *= q
753 x_0, y_0, z_0 = _remove_gcd(x_0, y_0, z_0)
755 # Holzer reduction
756 if sign(a) == sign(b):
757 x_0, y_0, z_0 = holzer(x_0, y_0, z_0, abs(a_2), abs(b_2), abs(c_2))
758 elif sign(a) == sign(c):
759 x_0, z_0, y_0 = holzer(x_0, z_0, y_0, abs(a_2), abs(c_2), abs(b_2))
760 else:
761 y_0, z_0, x_0 = holzer(y_0, z_0, x_0, abs(b_2), abs(c_2), abs(a_2))
763 x_0 = reconstruct(b_1, c_1, x_0)
764 y_0 = reconstruct(a_1, c_1, y_0)
765 z_0 = reconstruct(a_1, b_1, z_0)
767 sq_lcm = ilcm(sqf_of_a, sqf_of_b, sqf_of_c)
769 x_0 = abs(x_0*sq_lcm // sqf_of_a)
770 y_0 = abs(y_0*sq_lcm // sqf_of_b)
771 z_0 = abs(z_0*sq_lcm // sqf_of_c)
773 result.add(_remove_gcd(x_0, y_0, z_0))
774 return result
777class HomogeneousTernaryQuadratic(DiophantineEquationType):
778 """
779 Representation of a homogeneous ternary quadratic diophantine equation.
781 Examples
782 ========
784 >>> from sympy.abc import x, y, z
785 >>> from sympy.solvers.diophantine.diophantine import HomogeneousTernaryQuadratic
786 >>> HomogeneousTernaryQuadratic(x**2 + y**2 - 3*z**2 + x*y).solve()
787 {(-1, 2, 1)}
788 >>> HomogeneousTernaryQuadratic(3*x**2 + y**2 - 3*z**2 + 5*x*y + y*z).solve()
789 {(3, 12, 13)}
791 """
793 name = 'homogeneous_ternary_quadratic'
795 def matches(self):
796 if not (self.total_degree == 2 and self.dimension == 3):
797 return False
798 if not self.homogeneous:
799 return False
800 if not self.homogeneous_order:
801 return False
803 nonzero = [k for k in self.coeff if self.coeff[k]]
804 return not (len(nonzero) == 3 and all(i**2 in nonzero for i in self.free_symbols))
806 def solve(self, parameters=None, limit=None):
807 self.pre_solve(parameters)
809 _var = self.free_symbols
810 coeff = self.coeff
812 x, y, z = _var
813 var = [x, y, z]
815 # Equations of the form B*x*y + C*z*x + E*y*z = 0 and At least two of the
816 # coefficients A, B, C are non-zero.
817 # There are infinitely many solutions for the equation.
818 # Ex: (0, 0, t), (0, t, 0), (t, 0, 0)
819 # Equation can be re-written as y*(B*x + E*z) = -C*x*z and we can find rather
820 # unobvious solutions. Set y = -C and B*x + E*z = x*z. The latter can be solved by
821 # using methods for binary quadratic diophantine equations. Let's select the
822 # solution which minimizes |x| + |z|
824 result = DiophantineSolutionSet(var, parameters=self.parameters)
826 def unpack_sol(sol):
827 if len(sol) > 0:
828 return list(sol)[0]
829 return None, None, None
831 if not any(coeff[i**2] for i in var):
832 if coeff[x*z]:
833 sols = diophantine(coeff[x*y]*x + coeff[y*z]*z - x*z)
834 s = sols.pop()
835 min_sum = abs(s[0]) + abs(s[1])
837 for r in sols:
838 m = abs(r[0]) + abs(r[1])
839 if m < min_sum:
840 s = r
841 min_sum = m
843 result.add(_remove_gcd(s[0], -coeff[x*z], s[1]))
844 return result
846 else:
847 var[0], var[1] = _var[1], _var[0]
848 y_0, x_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
849 if x_0 is not None:
850 result.add((x_0, y_0, z_0))
851 return result
853 if coeff[x**2] == 0:
854 # If the coefficient of x is zero change the variables
855 if coeff[y**2] == 0:
856 var[0], var[2] = _var[2], _var[0]
857 z_0, y_0, x_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
859 else:
860 var[0], var[1] = _var[1], _var[0]
861 y_0, x_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
863 else:
864 if coeff[x*y] or coeff[x*z]:
865 # Apply the transformation x --> X - (B*y + C*z)/(2*A)
866 A = coeff[x**2]
867 B = coeff[x*y]
868 C = coeff[x*z]
869 D = coeff[y**2]
870 E = coeff[y*z]
871 F = coeff[z**2]
873 _coeff = {}
875 _coeff[x**2] = 4*A**2
876 _coeff[y**2] = 4*A*D - B**2
877 _coeff[z**2] = 4*A*F - C**2
878 _coeff[y*z] = 4*A*E - 2*B*C
879 _coeff[x*y] = 0
880 _coeff[x*z] = 0
882 x_0, y_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, _coeff))
884 if x_0 is None:
885 return result
887 p, q = _rational_pq(B*y_0 + C*z_0, 2*A)
888 x_0, y_0, z_0 = x_0*q - p, y_0*q, z_0*q
890 elif coeff[z*y] != 0:
891 if coeff[y**2] == 0:
892 if coeff[z**2] == 0:
893 # Equations of the form A*x**2 + E*yz = 0.
894 A = coeff[x**2]
895 E = coeff[y*z]
897 b, a = _rational_pq(-E, A)
899 x_0, y_0, z_0 = b, a, b
901 else:
902 # Ax**2 + E*y*z + F*z**2 = 0
903 var[0], var[2] = _var[2], _var[0]
904 z_0, y_0, x_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
906 else:
907 # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, C may be zero
908 var[0], var[1] = _var[1], _var[0]
909 y_0, x_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
911 else:
912 # Ax**2 + D*y**2 + F*z**2 = 0, C may be zero
913 x_0, y_0, z_0 = unpack_sol(_diop_ternary_quadratic_normal(var, coeff))
915 if x_0 is None:
916 return result
918 result.add(_remove_gcd(x_0, y_0, z_0))
919 return result
922class InhomogeneousGeneralQuadratic(DiophantineEquationType):
923 """
925 Representation of an inhomogeneous general quadratic.
927 No solver is currently implemented for this equation type.
929 """
931 name = 'inhomogeneous_general_quadratic'
933 def matches(self):
934 if not (self.total_degree == 2 and self.dimension >= 3):
935 return False
936 if not self.homogeneous_order:
937 return True
938 else:
939 # there may be Pow keys like x**2 or Mul keys like x*y
940 if any(k.is_Mul for k in self.coeff): # cross terms
941 return not self.homogeneous
942 return False
945class HomogeneousGeneralQuadratic(DiophantineEquationType):
946 """
948 Representation of a homogeneous general quadratic.
950 No solver is currently implemented for this equation type.
952 """
954 name = 'homogeneous_general_quadratic'
956 def matches(self):
957 if not (self.total_degree == 2 and self.dimension >= 3):
958 return False
959 if not self.homogeneous_order:
960 return False
961 else:
962 # there may be Pow keys like x**2 or Mul keys like x*y
963 if any(k.is_Mul for k in self.coeff): # cross terms
964 return self.homogeneous
965 return False
968class GeneralSumOfSquares(DiophantineEquationType):
969 r"""
970 Representation of the diophantine equation
972 `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`.
974 Details
975 =======
977 When `n = 3` if `k = 4^a(8m + 7)` for some `a, m \in Z` then there will be
978 no solutions. Refer [1]_ for more details.
980 Examples
981 ========
983 >>> from sympy.solvers.diophantine.diophantine import GeneralSumOfSquares
984 >>> from sympy.abc import a, b, c, d, e
985 >>> GeneralSumOfSquares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345).solve()
986 {(15, 22, 22, 24, 24)}
988 By default only 1 solution is returned. Use the `limit` keyword for more:
990 >>> sorted(GeneralSumOfSquares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345).solve(limit=3))
991 [(15, 22, 22, 24, 24), (16, 19, 24, 24, 24), (16, 20, 22, 23, 26)]
993 References
994 ==========
996 .. [1] Representing an integer as a sum of three squares, [online],
997 Available:
998 https://www.proofwiki.org/wiki/Integer_as_Sum_of_Three_Squares
999 """
1001 name = 'general_sum_of_squares'
1003 def matches(self):
1004 if not (self.total_degree == 2 and self.dimension >= 3):
1005 return False
1006 if not self.homogeneous_order:
1007 return False
1008 if any(k.is_Mul for k in self.coeff):
1009 return False
1010 return all(self.coeff[k] == 1 for k in self.coeff if k != 1)
1012 def solve(self, parameters=None, limit=1):
1013 self.pre_solve(parameters)
1015 var = self.free_symbols
1016 k = -int(self.coeff[1])
1017 n = self.dimension
1019 result = DiophantineSolutionSet(var, parameters=self.parameters)
1021 if k < 0 or limit < 1:
1022 return result
1024 signs = [-1 if x.is_nonpositive else 1 for x in var]
1025 negs = signs.count(-1) != 0
1027 took = 0
1028 for t in sum_of_squares(k, n, zeros=True):
1029 if negs:
1030 result.add([signs[i]*j for i, j in enumerate(t)])
1031 else:
1032 result.add(t)
1033 took += 1
1034 if took == limit:
1035 break
1036 return result
1039class GeneralPythagorean(DiophantineEquationType):
1040 """
1041 Representation of the general pythagorean equation,
1042 `a_{1}^2x_{1}^2 + a_{2}^2x_{2}^2 + . . . + a_{n}^2x_{n}^2 - a_{n + 1}^2x_{n + 1}^2 = 0`.
1044 Examples
1045 ========
1047 >>> from sympy.solvers.diophantine.diophantine import GeneralPythagorean
1048 >>> from sympy.abc import a, b, c, d, e, x, y, z, t
1049 >>> GeneralPythagorean(a**2 + b**2 + c**2 - d**2).solve()
1050 {(t_0**2 + t_1**2 - t_2**2, 2*t_0*t_2, 2*t_1*t_2, t_0**2 + t_1**2 + t_2**2)}
1051 >>> GeneralPythagorean(9*a**2 - 4*b**2 + 16*c**2 + 25*d**2 + e**2).solve(parameters=[x, y, z, t])
1052 {(-10*t**2 + 10*x**2 + 10*y**2 + 10*z**2, 15*t**2 + 15*x**2 + 15*y**2 + 15*z**2, 15*t*x, 12*t*y, 60*t*z)}
1053 """
1055 name = 'general_pythagorean'
1057 def matches(self):
1058 if not (self.total_degree == 2 and self.dimension >= 3):
1059 return False
1060 if not self.homogeneous_order:
1061 return False
1062 if any(k.is_Mul for k in self.coeff):
1063 return False
1064 if all(self.coeff[k] == 1 for k in self.coeff if k != 1):
1065 return False
1066 if not all(is_square(abs(self.coeff[k])) for k in self.coeff):
1067 return False
1068 # all but one has the same sign
1069 # e.g. 4*x**2 + y**2 - 4*z**2
1070 return abs(sum(sign(self.coeff[k]) for k in self.coeff)) == self.dimension - 2
1072 @property
1073 def n_parameters(self):
1074 return self.dimension - 1
1076 def solve(self, parameters=None, limit=1):
1077 self.pre_solve(parameters)
1079 coeff = self.coeff
1080 var = self.free_symbols
1081 n = self.dimension
1083 if sign(coeff[var[0] ** 2]) + sign(coeff[var[1] ** 2]) + sign(coeff[var[2] ** 2]) < 0:
1084 for key in coeff.keys():
1085 coeff[key] = -coeff[key]
1087 result = DiophantineSolutionSet(var, parameters=self.parameters)
1089 index = 0
1091 for i, v in enumerate(var):
1092 if sign(coeff[v ** 2]) == -1:
1093 index = i
1095 m = result.parameters
1097 ith = sum(m_i ** 2 for m_i in m)
1098 L = [ith - 2 * m[n - 2] ** 2]
1099 L.extend([2 * m[i] * m[n - 2] for i in range(n - 2)])
1100 sol = L[:index] + [ith] + L[index:]
1102 lcm = 1
1103 for i, v in enumerate(var):
1104 if i == index or (index > 0 and i == 0) or (index == 0 and i == 1):
1105 lcm = ilcm(lcm, sqrt(abs(coeff[v ** 2])))
1106 else:
1107 s = sqrt(coeff[v ** 2])
1108 lcm = ilcm(lcm, s if _odd(s) else s // 2)
1110 for i, v in enumerate(var):
1111 sol[i] = (lcm * sol[i]) / sqrt(abs(coeff[v ** 2]))
1113 result.add(sol)
1114 return result
1117class CubicThue(DiophantineEquationType):
1118 """
1119 Representation of a cubic Thue diophantine equation.
1121 A cubic Thue diophantine equation is a polynomial of the form
1122 `f(x, y) = r` of degree 3, where `x` and `y` are integers
1123 and `r` is a rational number.
1125 No solver is currently implemented for this equation type.
1127 Examples
1128 ========
1130 >>> from sympy.abc import x, y
1131 >>> from sympy.solvers.diophantine.diophantine import CubicThue
1132 >>> c1 = CubicThue(x**3 + y**2 + 1)
1133 >>> c1.matches()
1134 True
1136 """
1138 name = 'cubic_thue'
1140 def matches(self):
1141 return self.total_degree == 3 and self.dimension == 2
1144class GeneralSumOfEvenPowers(DiophantineEquationType):
1145 """
1146 Representation of the diophantine equation
1148 `x_{1}^e + x_{2}^e + . . . + x_{n}^e - k = 0`
1150 where `e` is an even, integer power.
1152 Examples
1153 ========
1155 >>> from sympy.solvers.diophantine.diophantine import GeneralSumOfEvenPowers
1156 >>> from sympy.abc import a, b
1157 >>> GeneralSumOfEvenPowers(a**4 + b**4 - (2**4 + 3**4)).solve()
1158 {(2, 3)}
1160 """
1162 name = 'general_sum_of_even_powers'
1164 def matches(self):
1165 if not self.total_degree > 3:
1166 return False
1167 if self.total_degree % 2 != 0:
1168 return False
1169 if not all(k.is_Pow and k.exp == self.total_degree for k in self.coeff if k != 1):
1170 return False
1171 return all(self.coeff[k] == 1 for k in self.coeff if k != 1)
1173 def solve(self, parameters=None, limit=1):
1174 self.pre_solve(parameters)
1176 var = self.free_symbols
1177 coeff = self.coeff
1179 p = None
1180 for q in coeff.keys():
1181 if q.is_Pow and coeff[q]:
1182 p = q.exp
1184 k = len(var)
1185 n = -coeff[1]
1187 result = DiophantineSolutionSet(var, parameters=self.parameters)
1189 if n < 0 or limit < 1:
1190 return result
1192 sign = [-1 if x.is_nonpositive else 1 for x in var]
1193 negs = sign.count(-1) != 0
1195 took = 0
1196 for t in power_representation(n, p, k):
1197 if negs:
1198 result.add([sign[i]*j for i, j in enumerate(t)])
1199 else:
1200 result.add(t)
1201 took += 1
1202 if took == limit:
1203 break
1204 return result
1206# these types are known (but not necessarily handled)
1207# note that order is important here (in the current solver state)
1208all_diop_classes = [
1209 Linear,
1210 Univariate,
1211 BinaryQuadratic,
1212 InhomogeneousTernaryQuadratic,
1213 HomogeneousTernaryQuadraticNormal,
1214 HomogeneousTernaryQuadratic,
1215 InhomogeneousGeneralQuadratic,
1216 HomogeneousGeneralQuadratic,
1217 GeneralSumOfSquares,
1218 GeneralPythagorean,
1219 CubicThue,
1220 GeneralSumOfEvenPowers,
1221]
1223diop_known = {diop_class.name for diop_class in all_diop_classes}
1226def _is_int(i):
1227 try:
1228 as_int(i)
1229 return True
1230 except ValueError:
1231 pass
1234def _sorted_tuple(*i):
1235 return tuple(sorted(i))
1238def _remove_gcd(*x):
1239 try:
1240 g = igcd(*x)
1241 except ValueError:
1242 fx = list(filter(None, x))
1243 if len(fx) < 2:
1244 return x
1245 g = igcd(*[i.as_content_primitive()[0] for i in fx])
1246 except TypeError:
1247 raise TypeError('_remove_gcd(a,b,c) or _remove_gcd(*container)')
1248 if g == 1:
1249 return x
1250 return tuple([i//g for i in x])
1253def _rational_pq(a, b):
1254 # return `(numer, denom)` for a/b; sign in numer and gcd removed
1255 return _remove_gcd(sign(b)*a, abs(b))
1258def _nint_or_floor(p, q):
1259 # return nearest int to p/q; in case of tie return floor(p/q)
1260 w, r = divmod(p, q)
1261 if abs(r) <= abs(q)//2:
1262 return w
1263 return w + 1
1266def _odd(i):
1267 return i % 2 != 0
1270def _even(i):
1271 return i % 2 == 0
1274def diophantine(eq, param=symbols("t", integer=True), syms=None,
1275 permute=False):
1276 """
1277 Simplify the solution procedure of diophantine equation ``eq`` by
1278 converting it into a product of terms which should equal zero.
1280 Explanation
1281 ===========
1283 For example, when solving, `x^2 - y^2 = 0` this is treated as
1284 `(x + y)(x - y) = 0` and `x + y = 0` and `x - y = 0` are solved
1285 independently and combined. Each term is solved by calling
1286 ``diop_solve()``. (Although it is possible to call ``diop_solve()``
1287 directly, one must be careful to pass an equation in the correct
1288 form and to interpret the output correctly; ``diophantine()`` is
1289 the public-facing function to use in general.)
1291 Output of ``diophantine()`` is a set of tuples. The elements of the
1292 tuple are the solutions for each variable in the equation and
1293 are arranged according to the alphabetic ordering of the variables.
1294 e.g. For an equation with two variables, `a` and `b`, the first
1295 element of the tuple is the solution for `a` and the second for `b`.
1297 Usage
1298 =====
1300 ``diophantine(eq, t, syms)``: Solve the diophantine
1301 equation ``eq``.
1302 ``t`` is the optional parameter to be used by ``diop_solve()``.
1303 ``syms`` is an optional list of symbols which determines the
1304 order of the elements in the returned tuple.
1306 By default, only the base solution is returned. If ``permute`` is set to
1307 True then permutations of the base solution and/or permutations of the
1308 signs of the values will be returned when applicable.
1310 Details
1311 =======
1313 ``eq`` should be an expression which is assumed to be zero.
1314 ``t`` is the parameter to be used in the solution.
1316 Examples
1317 ========
1319 >>> from sympy import diophantine
1320 >>> from sympy.abc import a, b
1321 >>> eq = a**4 + b**4 - (2**4 + 3**4)
1322 >>> diophantine(eq)
1323 {(2, 3)}
1324 >>> diophantine(eq, permute=True)
1325 {(-3, -2), (-3, 2), (-2, -3), (-2, 3), (2, -3), (2, 3), (3, -2), (3, 2)}
1327 >>> from sympy.abc import x, y, z
1328 >>> diophantine(x**2 - y**2)
1329 {(t_0, -t_0), (t_0, t_0)}
1331 >>> diophantine(x*(2*x + 3*y - z))
1332 {(0, n1, n2), (t_0, t_1, 2*t_0 + 3*t_1)}
1333 >>> diophantine(x**2 + 3*x*y + 4*x)
1334 {(0, n1), (3*t_0 - 4, -t_0)}
1336 See Also
1337 ========
1339 diop_solve
1340 sympy.utilities.iterables.permute_signs
1341 sympy.utilities.iterables.signed_permutations
1342 """
1344 eq = _sympify(eq)
1346 if isinstance(eq, Eq):
1347 eq = eq.lhs - eq.rhs
1349 try:
1350 var = list(eq.expand(force=True).free_symbols)
1351 var.sort(key=default_sort_key)
1352 if syms:
1353 if not is_sequence(syms):
1354 raise TypeError(
1355 'syms should be given as a sequence, e.g. a list')
1356 syms = [i for i in syms if i in var]
1357 if syms != var:
1358 dict_sym_index = dict(zip(syms, range(len(syms))))
1359 return {tuple([t[dict_sym_index[i]] for i in var])
1360 for t in diophantine(eq, param, permute=permute)}
1361 n, d = eq.as_numer_denom()
1362 if n.is_number:
1363 return set()
1364 if not d.is_number:
1365 dsol = diophantine(d)
1366 good = diophantine(n) - dsol
1367 return {s for s in good if _mexpand(d.subs(zip(var, s)))}
1368 else:
1369 eq = n
1370 eq = factor_terms(eq)
1371 assert not eq.is_number
1372 eq = eq.as_independent(*var, as_Add=False)[1]
1373 p = Poly(eq)
1374 assert not any(g.is_number for g in p.gens)
1375 eq = p.as_expr()
1376 assert eq.is_polynomial()
1377 except (GeneratorsNeeded, AssertionError):
1378 raise TypeError(filldedent('''
1379 Equation should be a polynomial with Rational coefficients.'''))
1381 # permute only sign
1382 do_permute_signs = False
1383 # permute sign and values
1384 do_permute_signs_var = False
1385 # permute few signs
1386 permute_few_signs = False
1387 try:
1388 # if we know that factoring should not be attempted, skip
1389 # the factoring step
1390 v, c, t = classify_diop(eq)
1392 # check for permute sign
1393 if permute:
1394 len_var = len(v)
1395 permute_signs_for = [
1396 GeneralSumOfSquares.name,
1397 GeneralSumOfEvenPowers.name]
1398 permute_signs_check = [
1399 HomogeneousTernaryQuadratic.name,
1400 HomogeneousTernaryQuadraticNormal.name,
1401 BinaryQuadratic.name]
1402 if t in permute_signs_for:
1403 do_permute_signs_var = True
1404 elif t in permute_signs_check:
1405 # if all the variables in eq have even powers
1406 # then do_permute_sign = True
1407 if len_var == 3:
1408 var_mul = list(subsets(v, 2))
1409 # here var_mul is like [(x, y), (x, z), (y, z)]
1410 xy_coeff = True
1411 x_coeff = True
1412 var1_mul_var2 = (a[0]*a[1] for a in var_mul)
1413 # if coeff(y*z), coeff(y*x), coeff(x*z) is not 0 then
1414 # `xy_coeff` => True and do_permute_sign => False.
1415 # Means no permuted solution.
1416 for v1_mul_v2 in var1_mul_var2:
1417 try:
1418 coeff = c[v1_mul_v2]
1419 except KeyError:
1420 coeff = 0
1421 xy_coeff = bool(xy_coeff) and bool(coeff)
1422 var_mul = list(subsets(v, 1))
1423 # here var_mul is like [(x,), (y, )]
1424 for v1 in var_mul:
1425 try:
1426 coeff = c[v1[0]]
1427 except KeyError:
1428 coeff = 0
1429 x_coeff = bool(x_coeff) and bool(coeff)
1430 if not any((xy_coeff, x_coeff)):
1431 # means only x**2, y**2, z**2, const is present
1432 do_permute_signs = True
1433 elif not x_coeff:
1434 permute_few_signs = True
1435 elif len_var == 2:
1436 var_mul = list(subsets(v, 2))
1437 # here var_mul is like [(x, y)]
1438 xy_coeff = True
1439 x_coeff = True
1440 var1_mul_var2 = (x[0]*x[1] for x in var_mul)
1441 for v1_mul_v2 in var1_mul_var2:
1442 try:
1443 coeff = c[v1_mul_v2]
1444 except KeyError:
1445 coeff = 0
1446 xy_coeff = bool(xy_coeff) and bool(coeff)
1447 var_mul = list(subsets(v, 1))
1448 # here var_mul is like [(x,), (y, )]
1449 for v1 in var_mul:
1450 try:
1451 coeff = c[v1[0]]
1452 except KeyError:
1453 coeff = 0
1454 x_coeff = bool(x_coeff) and bool(coeff)
1455 if not any((xy_coeff, x_coeff)):
1456 # means only x**2, y**2 and const is present
1457 # so we can get more soln by permuting this soln.
1458 do_permute_signs = True
1459 elif not x_coeff:
1460 # when coeff(x), coeff(y) is not present then signs of
1461 # x, y can be permuted such that their sign are same
1462 # as sign of x*y.
1463 # e.g 1. (x_val,y_val)=> (x_val,y_val), (-x_val,-y_val)
1464 # 2. (-x_vall, y_val)=> (-x_val,y_val), (x_val,-y_val)
1465 permute_few_signs = True
1466 if t == 'general_sum_of_squares':
1467 # trying to factor such expressions will sometimes hang
1468 terms = [(eq, 1)]
1469 else:
1470 raise TypeError
1471 except (TypeError, NotImplementedError):
1472 fl = factor_list(eq)
1473 if fl[0].is_Rational and fl[0] != 1:
1474 return diophantine(eq/fl[0], param=param, syms=syms, permute=permute)
1475 terms = fl[1]
1477 sols = set()
1479 for term in terms:
1481 base, _ = term
1482 var_t, _, eq_type = classify_diop(base, _dict=False)
1483 _, base = signsimp(base, evaluate=False).as_coeff_Mul()
1484 solution = diop_solve(base, param)
1486 if eq_type in [
1487 Linear.name,
1488 HomogeneousTernaryQuadratic.name,
1489 HomogeneousTernaryQuadraticNormal.name,
1490 GeneralPythagorean.name]:
1491 sols.add(merge_solution(var, var_t, solution))
1493 elif eq_type in [
1494 BinaryQuadratic.name,
1495 GeneralSumOfSquares.name,
1496 GeneralSumOfEvenPowers.name,
1497 Univariate.name]:
1498 for sol in solution:
1499 sols.add(merge_solution(var, var_t, sol))
1501 else:
1502 raise NotImplementedError('unhandled type: %s' % eq_type)
1504 # remove null merge results
1505 if () in sols:
1506 sols.remove(())
1507 null = tuple([0]*len(var))
1508 # if there is no solution, return trivial solution
1509 if not sols and eq.subs(zip(var, null)).is_zero:
1510 sols.add(null)
1511 final_soln = set()
1512 for sol in sols:
1513 if all(_is_int(s) for s in sol):
1514 if do_permute_signs:
1515 permuted_sign = set(permute_signs(sol))
1516 final_soln.update(permuted_sign)
1517 elif permute_few_signs:
1518 lst = list(permute_signs(sol))
1519 lst = list(filter(lambda x: x[0]*x[1] == sol[1]*sol[0], lst))
1520 permuted_sign = set(lst)
1521 final_soln.update(permuted_sign)
1522 elif do_permute_signs_var:
1523 permuted_sign_var = set(signed_permutations(sol))
1524 final_soln.update(permuted_sign_var)
1525 else:
1526 final_soln.add(sol)
1527 else:
1528 final_soln.add(sol)
1529 return final_soln
1532def merge_solution(var, var_t, solution):
1533 """
1534 This is used to construct the full solution from the solutions of sub
1535 equations.
1537 Explanation
1538 ===========
1540 For example when solving the equation `(x - y)(x^2 + y^2 - z^2) = 0`,
1541 solutions for each of the equations `x - y = 0` and `x^2 + y^2 - z^2` are
1542 found independently. Solutions for `x - y = 0` are `(x, y) = (t, t)`. But
1543 we should introduce a value for z when we output the solution for the
1544 original equation. This function converts `(t, t)` into `(t, t, n_{1})`
1545 where `n_{1}` is an integer parameter.
1546 """
1547 sol = []
1549 if None in solution:
1550 return ()
1552 solution = iter(solution)
1553 params = numbered_symbols("n", integer=True, start=1)
1554 for v in var:
1555 if v in var_t:
1556 sol.append(next(solution))
1557 else:
1558 sol.append(next(params))
1560 for val, symb in zip(sol, var):
1561 if check_assumptions(val, **symb.assumptions0) is False:
1562 return ()
1564 return tuple(sol)
1567def _diop_solve(eq, params=None):
1568 for diop_type in all_diop_classes:
1569 if diop_type(eq).matches():
1570 return diop_type(eq).solve(parameters=params)
1573def diop_solve(eq, param=symbols("t", integer=True)):
1574 """
1575 Solves the diophantine equation ``eq``.
1577 Explanation
1578 ===========
1580 Unlike ``diophantine()``, factoring of ``eq`` is not attempted. Uses
1581 ``classify_diop()`` to determine the type of the equation and calls
1582 the appropriate solver function.
1584 Use of ``diophantine()`` is recommended over other helper functions.
1585 ``diop_solve()`` can return either a set or a tuple depending on the
1586 nature of the equation.
1588 Usage
1589 =====
1591 ``diop_solve(eq, t)``: Solve diophantine equation, ``eq`` using ``t``
1592 as a parameter if needed.
1594 Details
1595 =======
1597 ``eq`` should be an expression which is assumed to be zero.
1598 ``t`` is a parameter to be used in the solution.
1600 Examples
1601 ========
1603 >>> from sympy.solvers.diophantine import diop_solve
1604 >>> from sympy.abc import x, y, z, w
1605 >>> diop_solve(2*x + 3*y - 5)
1606 (3*t_0 - 5, 5 - 2*t_0)
1607 >>> diop_solve(4*x + 3*y - 4*z + 5)
1608 (t_0, 8*t_0 + 4*t_1 + 5, 7*t_0 + 3*t_1 + 5)
1609 >>> diop_solve(x + 3*y - 4*z + w - 6)
1610 (t_0, t_0 + t_1, 6*t_0 + 5*t_1 + 4*t_2 - 6, 5*t_0 + 4*t_1 + 3*t_2 - 6)
1611 >>> diop_solve(x**2 + y**2 - 5)
1612 {(-2, -1), (-2, 1), (-1, -2), (-1, 2), (1, -2), (1, 2), (2, -1), (2, 1)}
1615 See Also
1616 ========
1618 diophantine()
1619 """
1620 var, coeff, eq_type = classify_diop(eq, _dict=False)
1622 if eq_type == Linear.name:
1623 return diop_linear(eq, param)
1625 elif eq_type == BinaryQuadratic.name:
1626 return diop_quadratic(eq, param)
1628 elif eq_type == HomogeneousTernaryQuadratic.name:
1629 return diop_ternary_quadratic(eq, parameterize=True)
1631 elif eq_type == HomogeneousTernaryQuadraticNormal.name:
1632 return diop_ternary_quadratic_normal(eq, parameterize=True)
1634 elif eq_type == GeneralPythagorean.name:
1635 return diop_general_pythagorean(eq, param)
1637 elif eq_type == Univariate.name:
1638 return diop_univariate(eq)
1640 elif eq_type == GeneralSumOfSquares.name:
1641 return diop_general_sum_of_squares(eq, limit=S.Infinity)
1643 elif eq_type == GeneralSumOfEvenPowers.name:
1644 return diop_general_sum_of_even_powers(eq, limit=S.Infinity)
1646 if eq_type is not None and eq_type not in diop_known:
1647 raise ValueError(filldedent('''
1648 Although this type of equation was identified, it is not yet
1649 handled. It should, however, be listed in `diop_known` at the
1650 top of this file. Developers should see comments at the end of
1651 `classify_diop`.
1652 ''')) # pragma: no cover
1653 else:
1654 raise NotImplementedError(
1655 'No solver has been written for %s.' % eq_type)
1658def classify_diop(eq, _dict=True):
1659 # docstring supplied externally
1661 matched = False
1662 diop_type = None
1663 for diop_class in all_diop_classes:
1664 diop_type = diop_class(eq)
1665 if diop_type.matches():
1666 matched = True
1667 break
1669 if matched:
1670 return diop_type.free_symbols, dict(diop_type.coeff) if _dict else diop_type.coeff, diop_type.name
1672 # new diop type instructions
1673 # --------------------------
1674 # if this error raises and the equation *can* be classified,
1675 # * it should be identified in the if-block above
1676 # * the type should be added to the diop_known
1677 # if a solver can be written for it,
1678 # * a dedicated handler should be written (e.g. diop_linear)
1679 # * it should be passed to that handler in diop_solve
1680 raise NotImplementedError(filldedent('''
1681 This equation is not yet recognized or else has not been
1682 simplified sufficiently to put it in a form recognized by
1683 diop_classify().'''))
1686classify_diop.func_doc = ( # type: ignore
1687 '''
1688 Helper routine used by diop_solve() to find information about ``eq``.
1690 Explanation
1691 ===========
1693 Returns a tuple containing the type of the diophantine equation
1694 along with the variables (free symbols) and their coefficients.
1695 Variables are returned as a list and coefficients are returned
1696 as a dict with the key being the respective term and the constant
1697 term is keyed to 1. The type is one of the following:
1699 * %s
1701 Usage
1702 =====
1704 ``classify_diop(eq)``: Return variables, coefficients and type of the
1705 ``eq``.
1707 Details
1708 =======
1710 ``eq`` should be an expression which is assumed to be zero.
1711 ``_dict`` is for internal use: when True (default) a dict is returned,
1712 otherwise a defaultdict which supplies 0 for missing keys is returned.
1714 Examples
1715 ========
1717 >>> from sympy.solvers.diophantine import classify_diop
1718 >>> from sympy.abc import x, y, z, w, t
1719 >>> classify_diop(4*x + 6*y - 4)
1720 ([x, y], {1: -4, x: 4, y: 6}, 'linear')
1721 >>> classify_diop(x + 3*y -4*z + 5)
1722 ([x, y, z], {1: 5, x: 1, y: 3, z: -4}, 'linear')
1723 >>> classify_diop(x**2 + y**2 - x*y + x + 5)
1724 ([x, y], {1: 5, x: 1, x**2: 1, y**2: 1, x*y: -1}, 'binary_quadratic')
1725 ''' % ('\n * '.join(sorted(diop_known))))
1728def diop_linear(eq, param=symbols("t", integer=True)):
1729 """
1730 Solves linear diophantine equations.
1732 A linear diophantine equation is an equation of the form `a_{1}x_{1} +
1733 a_{2}x_{2} + .. + a_{n}x_{n} = 0` where `a_{1}, a_{2}, ..a_{n}` are
1734 integer constants and `x_{1}, x_{2}, ..x_{n}` are integer variables.
1736 Usage
1737 =====
1739 ``diop_linear(eq)``: Returns a tuple containing solutions to the
1740 diophantine equation ``eq``. Values in the tuple is arranged in the same
1741 order as the sorted variables.
1743 Details
1744 =======
1746 ``eq`` is a linear diophantine equation which is assumed to be zero.
1747 ``param`` is the parameter to be used in the solution.
1749 Examples
1750 ========
1752 >>> from sympy.solvers.diophantine.diophantine import diop_linear
1753 >>> from sympy.abc import x, y, z
1754 >>> diop_linear(2*x - 3*y - 5) # solves equation 2*x - 3*y - 5 == 0
1755 (3*t_0 - 5, 2*t_0 - 5)
1757 Here x = -3*t_0 - 5 and y = -2*t_0 - 5
1759 >>> diop_linear(2*x - 3*y - 4*z -3)
1760 (t_0, 2*t_0 + 4*t_1 + 3, -t_0 - 3*t_1 - 3)
1762 See Also
1763 ========
1765 diop_quadratic(), diop_ternary_quadratic(), diop_general_pythagorean(),
1766 diop_general_sum_of_squares()
1767 """
1768 var, coeff, diop_type = classify_diop(eq, _dict=False)
1770 if diop_type == Linear.name:
1771 parameters = None
1772 if param is not None:
1773 parameters = symbols('%s_0:%i' % (param, len(var)), integer=True)
1775 result = Linear(eq).solve(parameters=parameters)
1777 if param is None:
1778 result = result(*[0]*len(result.parameters))
1780 if len(result) > 0:
1781 return list(result)[0]
1782 else:
1783 return tuple([None]*len(result.parameters))
1786def base_solution_linear(c, a, b, t=None):
1787 """
1788 Return the base solution for the linear equation, `ax + by = c`.
1790 Explanation
1791 ===========
1793 Used by ``diop_linear()`` to find the base solution of a linear
1794 Diophantine equation. If ``t`` is given then the parametrized solution is
1795 returned.
1797 Usage
1798 =====
1800 ``base_solution_linear(c, a, b, t)``: ``a``, ``b``, ``c`` are coefficients
1801 in `ax + by = c` and ``t`` is the parameter to be used in the solution.
1803 Examples
1804 ========
1806 >>> from sympy.solvers.diophantine.diophantine import base_solution_linear
1807 >>> from sympy.abc import t
1808 >>> base_solution_linear(5, 2, 3) # equation 2*x + 3*y = 5
1809 (-5, 5)
1810 >>> base_solution_linear(0, 5, 7) # equation 5*x + 7*y = 0
1811 (0, 0)
1812 >>> base_solution_linear(5, 2, 3, t) # equation 2*x + 3*y = 5
1813 (3*t - 5, 5 - 2*t)
1814 >>> base_solution_linear(0, 5, 7, t) # equation 5*x + 7*y = 0
1815 (7*t, -5*t)
1816 """
1817 a, b, c = _remove_gcd(a, b, c)
1819 if c == 0:
1820 if t is not None:
1821 if b < 0:
1822 t = -t
1823 return (b*t, -a*t)
1824 else:
1825 return (0, 0)
1826 else:
1827 x0, y0, d = igcdex(abs(a), abs(b))
1829 x0 *= sign(a)
1830 y0 *= sign(b)
1832 if divisible(c, d):
1833 if t is not None:
1834 if b < 0:
1835 t = -t
1836 return (c*x0 + b*t, c*y0 - a*t)
1837 else:
1838 return (c*x0, c*y0)
1839 else:
1840 return (None, None)
1843def diop_univariate(eq):
1844 """
1845 Solves a univariate diophantine equations.
1847 Explanation
1848 ===========
1850 A univariate diophantine equation is an equation of the form
1851 `a_{0} + a_{1}x + a_{2}x^2 + .. + a_{n}x^n = 0` where `a_{1}, a_{2}, ..a_{n}` are
1852 integer constants and `x` is an integer variable.
1854 Usage
1855 =====
1857 ``diop_univariate(eq)``: Returns a set containing solutions to the
1858 diophantine equation ``eq``.
1860 Details
1861 =======
1863 ``eq`` is a univariate diophantine equation which is assumed to be zero.
1865 Examples
1866 ========
1868 >>> from sympy.solvers.diophantine.diophantine import diop_univariate
1869 >>> from sympy.abc import x
1870 >>> diop_univariate((x - 2)*(x - 3)**2) # solves equation (x - 2)*(x - 3)**2 == 0
1871 {(2,), (3,)}
1873 """
1874 var, coeff, diop_type = classify_diop(eq, _dict=False)
1876 if diop_type == Univariate.name:
1877 return {(int(i),) for i in solveset_real(
1878 eq, var[0]).intersect(S.Integers)}
1881def divisible(a, b):
1882 """
1883 Returns `True` if ``a`` is divisible by ``b`` and `False` otherwise.
1884 """
1885 return not a % b
1888def diop_quadratic(eq, param=symbols("t", integer=True)):
1889 """
1890 Solves quadratic diophantine equations.
1892 i.e. equations of the form `Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0`. Returns a
1893 set containing the tuples `(x, y)` which contains the solutions. If there
1894 are no solutions then `(None, None)` is returned.
1896 Usage
1897 =====
1899 ``diop_quadratic(eq, param)``: ``eq`` is a quadratic binary diophantine
1900 equation. ``param`` is used to indicate the parameter to be used in the
1901 solution.
1903 Details
1904 =======
1906 ``eq`` should be an expression which is assumed to be zero.
1907 ``param`` is a parameter to be used in the solution.
1909 Examples
1910 ========
1912 >>> from sympy.abc import x, y, t
1913 >>> from sympy.solvers.diophantine.diophantine import diop_quadratic
1914 >>> diop_quadratic(x**2 + y**2 + 2*x + 2*y + 2, t)
1915 {(-1, -1)}
1917 References
1918 ==========
1920 .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
1921 Available: https://www.alpertron.com.ar/METHODS.HTM
1922 .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
1923 Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
1925 See Also
1926 ========
1928 diop_linear(), diop_ternary_quadratic(), diop_general_sum_of_squares(),
1929 diop_general_pythagorean()
1930 """
1931 var, coeff, diop_type = classify_diop(eq, _dict=False)
1933 if diop_type == BinaryQuadratic.name:
1934 if param is not None:
1935 parameters = [param, Symbol("u", integer=True)]
1936 else:
1937 parameters = None
1938 return set(BinaryQuadratic(eq).solve(parameters=parameters))
1941def is_solution_quad(var, coeff, u, v):
1942 """
1943 Check whether `(u, v)` is solution to the quadratic binary diophantine
1944 equation with the variable list ``var`` and coefficient dictionary
1945 ``coeff``.
1947 Not intended for use by normal users.
1948 """
1949 reps = dict(zip(var, (u, v)))
1950 eq = Add(*[j*i.xreplace(reps) for i, j in coeff.items()])
1951 return _mexpand(eq) == 0
1954def diop_DN(D, N, t=symbols("t", integer=True)):
1955 """
1956 Solves the equation `x^2 - Dy^2 = N`.
1958 Explanation
1959 ===========
1961 Mainly concerned with the case `D > 0, D` is not a perfect square,
1962 which is the same as the generalized Pell equation. The LMM
1963 algorithm [1]_ is used to solve this equation.
1965 Returns one solution tuple, (`x, y)` for each class of the solutions.
1966 Other solutions of the class can be constructed according to the
1967 values of ``D`` and ``N``.
1969 Usage
1970 =====
1972 ``diop_DN(D, N, t)``: D and N are integers as in `x^2 - Dy^2 = N` and
1973 ``t`` is the parameter to be used in the solutions.
1975 Details
1976 =======
1978 ``D`` and ``N`` correspond to D and N in the equation.
1979 ``t`` is the parameter to be used in the solutions.
1981 Examples
1982 ========
1984 >>> from sympy.solvers.diophantine.diophantine import diop_DN
1985 >>> diop_DN(13, -4) # Solves equation x**2 - 13*y**2 = -4
1986 [(3, 1), (393, 109), (36, 10)]
1988 The output can be interpreted as follows: There are three fundamental
1989 solutions to the equation `x^2 - 13y^2 = -4` given by (3, 1), (393, 109)
1990 and (36, 10). Each tuple is in the form (x, y), i.e. solution (3, 1) means
1991 that `x = 3` and `y = 1`.
1993 >>> diop_DN(986, 1) # Solves equation x**2 - 986*y**2 = 1
1994 [(49299, 1570)]
1996 See Also
1997 ========
1999 find_DN(), diop_bf_DN()
2001 References
2002 ==========
2004 .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
2005 Robertson, July 31, 2004, Pages 16 - 17. [online], Available:
2006 https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
2007 """
2008 if D < 0:
2009 if N == 0:
2010 return [(0, 0)]
2011 elif N < 0:
2012 return []
2013 elif N > 0:
2014 sol = []
2015 for d in divisors(square_factor(N)):
2016 sols = cornacchia(1, -D, N // d**2)
2017 if sols:
2018 for x, y in sols:
2019 sol.append((d*x, d*y))
2020 if D == -1:
2021 sol.append((d*y, d*x))
2022 return sol
2024 elif D == 0:
2025 if N < 0:
2026 return []
2027 if N == 0:
2028 return [(0, t)]
2029 sN, _exact = integer_nthroot(N, 2)
2030 if _exact:
2031 return [(sN, t)]
2032 else:
2033 return []
2035 else: # D > 0
2036 sD, _exact = integer_nthroot(D, 2)
2037 if _exact:
2038 if N == 0:
2039 return [(sD*t, t)]
2040 else:
2041 sol = []
2043 for y in range(floor(sign(N)*(N - 1)/(2*sD)) + 1):
2044 try:
2045 sq, _exact = integer_nthroot(D*y**2 + N, 2)
2046 except ValueError:
2047 _exact = False
2048 if _exact:
2049 sol.append((sq, y))
2051 return sol
2053 elif 1 < N**2 < D:
2054 # It is much faster to call `_special_diop_DN`.
2055 return _special_diop_DN(D, N)
2057 else:
2058 if N == 0:
2059 return [(0, 0)]
2061 elif abs(N) == 1:
2063 pqa = PQa(0, 1, D)
2064 j = 0
2065 G = []
2066 B = []
2068 for i in pqa:
2070 a = i[2]
2071 G.append(i[5])
2072 B.append(i[4])
2074 if j != 0 and a == 2*sD:
2075 break
2076 j = j + 1
2078 if _odd(j):
2080 if N == -1:
2081 x = G[j - 1]
2082 y = B[j - 1]
2083 else:
2084 count = j
2085 while count < 2*j - 1:
2086 i = next(pqa)
2087 G.append(i[5])
2088 B.append(i[4])
2089 count += 1
2091 x = G[count]
2092 y = B[count]
2093 else:
2094 if N == 1:
2095 x = G[j - 1]
2096 y = B[j - 1]
2097 else:
2098 return []
2100 return [(x, y)]
2102 else:
2104 fs = []
2105 sol = []
2106 div = divisors(N)
2108 for d in div:
2109 if divisible(N, d**2):
2110 fs.append(d)
2112 for f in fs:
2113 m = N // f**2
2115 zs = sqrt_mod(D, abs(m), all_roots=True)
2116 zs = [i for i in zs if i <= abs(m) // 2 ]
2118 if abs(m) != 2:
2119 zs = zs + [-i for i in zs if i] # omit dupl 0
2121 for z in zs:
2123 pqa = PQa(z, abs(m), D)
2124 j = 0
2125 G = []
2126 B = []
2128 for i in pqa:
2130 G.append(i[5])
2131 B.append(i[4])
2133 if j != 0 and abs(i[1]) == 1:
2134 r = G[j-1]
2135 s = B[j-1]
2137 if r**2 - D*s**2 == m:
2138 sol.append((f*r, f*s))
2140 elif diop_DN(D, -1) != []:
2141 a = diop_DN(D, -1)
2142 sol.append((f*(r*a[0][0] + a[0][1]*s*D), f*(r*a[0][1] + s*a[0][0])))
2144 break
2146 j = j + 1
2147 if j == length(z, abs(m), D):
2148 break
2150 return sol
2153def _special_diop_DN(D, N):
2154 """
2155 Solves the equation `x^2 - Dy^2 = N` for the special case where
2156 `1 < N**2 < D` and `D` is not a perfect square.
2157 It is better to call `diop_DN` rather than this function, as
2158 the former checks the condition `1 < N**2 < D`, and calls the latter only
2159 if appropriate.
2161 Usage
2162 =====
2164 WARNING: Internal method. Do not call directly!
2166 ``_special_diop_DN(D, N)``: D and N are integers as in `x^2 - Dy^2 = N`.
2168 Details
2169 =======
2171 ``D`` and ``N`` correspond to D and N in the equation.
2173 Examples
2174 ========
2176 >>> from sympy.solvers.diophantine.diophantine import _special_diop_DN
2177 >>> _special_diop_DN(13, -3) # Solves equation x**2 - 13*y**2 = -3
2178 [(7, 2), (137, 38)]
2180 The output can be interpreted as follows: There are two fundamental
2181 solutions to the equation `x^2 - 13y^2 = -3` given by (7, 2) and
2182 (137, 38). Each tuple is in the form (x, y), i.e. solution (7, 2) means
2183 that `x = 7` and `y = 2`.
2185 >>> _special_diop_DN(2445, -20) # Solves equation x**2 - 2445*y**2 = -20
2186 [(445, 9), (17625560, 356454), (698095554475, 14118073569)]
2188 See Also
2189 ========
2191 diop_DN()
2193 References
2194 ==========
2196 .. [1] Section 4.4.4 of the following book:
2197 Quadratic Diophantine Equations, T. Andreescu and D. Andrica,
2198 Springer, 2015.
2199 """
2201 # The following assertion was removed for efficiency, with the understanding
2202 # that this method is not called directly. The parent method, `diop_DN`
2203 # is responsible for performing the appropriate checks.
2204 #
2205 # assert (1 < N**2 < D) and (not integer_nthroot(D, 2)[1])
2207 sqrt_D = sqrt(D)
2208 F = [(N, 1)]
2209 f = 2
2210 while True:
2211 f2 = f**2
2212 if f2 > abs(N):
2213 break
2214 n, r = divmod(N, f2)
2215 if r == 0:
2216 F.append((n, f))
2217 f += 1
2219 P = 0
2220 Q = 1
2221 G0, G1 = 0, 1
2222 B0, B1 = 1, 0
2224 solutions = []
2226 i = 0
2227 while True:
2228 a = floor((P + sqrt_D) / Q)
2229 P = a*Q - P
2230 Q = (D - P**2) // Q
2231 G2 = a*G1 + G0
2232 B2 = a*B1 + B0
2234 for n, f in F:
2235 if G2**2 - D*B2**2 == n:
2236 solutions.append((f*G2, f*B2))
2238 i += 1
2239 if Q == 1 and i % 2 == 0:
2240 break
2242 G0, G1 = G1, G2
2243 B0, B1 = B1, B2
2245 return solutions
2248def cornacchia(a, b, m):
2249 r"""
2250 Solves `ax^2 + by^2 = m` where `\gcd(a, b) = 1 = gcd(a, m)` and `a, b > 0`.
2252 Explanation
2253 ===========
2255 Uses the algorithm due to Cornacchia. The method only finds primitive
2256 solutions, i.e. ones with `\gcd(x, y) = 1`. So this method cannot be used to
2257 find the solutions of `x^2 + y^2 = 20` since the only solution to former is
2258 `(x, y) = (4, 2)` and it is not primitive. When `a = b`, only the
2259 solutions with `x \leq y` are found. For more details, see the References.
2261 Examples
2262 ========
2264 >>> from sympy.solvers.diophantine.diophantine import cornacchia
2265 >>> cornacchia(2, 3, 35) # equation 2x**2 + 3y**2 = 35
2266 {(2, 3), (4, 1)}
2267 >>> cornacchia(1, 1, 25) # equation x**2 + y**2 = 25
2268 {(4, 3)}
2270 References
2271 ===========
2273 .. [1] A. Nitaj, "L'algorithme de Cornacchia"
2274 .. [2] Solving the diophantine equation ax**2 + by**2 = m by Cornacchia's
2275 method, [online], Available:
2276 http://www.numbertheory.org/php/cornacchia.html
2278 See Also
2279 ========
2281 sympy.utilities.iterables.signed_permutations
2282 """
2283 sols = set()
2285 a1 = igcdex(a, m)[0]
2286 v = sqrt_mod(-b*a1, m, all_roots=True)
2287 if not v:
2288 return None
2290 for t in v:
2291 if t < m // 2:
2292 continue
2294 u, r = t, m
2296 while True:
2297 u, r = r, u % r
2298 if a*r**2 < m:
2299 break
2301 m1 = m - a*r**2
2303 if m1 % b == 0:
2304 m1 = m1 // b
2305 s, _exact = integer_nthroot(m1, 2)
2306 if _exact:
2307 if a == b and r < s:
2308 r, s = s, r
2309 sols.add((int(r), int(s)))
2311 return sols
2314def PQa(P_0, Q_0, D):
2315 r"""
2316 Returns useful information needed to solve the Pell equation.
2318 Explanation
2319 ===========
2321 There are six sequences of integers defined related to the continued
2322 fraction representation of `\\frac{P + \sqrt{D}}{Q}`, namely {`P_{i}`},
2323 {`Q_{i}`}, {`a_{i}`},{`A_{i}`}, {`B_{i}`}, {`G_{i}`}. ``PQa()`` Returns
2324 these values as a 6-tuple in the same order as mentioned above. Refer [1]_
2325 for more detailed information.
2327 Usage
2328 =====
2330 ``PQa(P_0, Q_0, D)``: ``P_0``, ``Q_0`` and ``D`` are integers corresponding
2331 to `P_{0}`, `Q_{0}` and `D` in the continued fraction
2332 `\\frac{P_{0} + \sqrt{D}}{Q_{0}}`.
2333 Also it's assumed that `P_{0}^2 == D mod(|Q_{0}|)` and `D` is square free.
2335 Examples
2336 ========
2338 >>> from sympy.solvers.diophantine.diophantine import PQa
2339 >>> pqa = PQa(13, 4, 5) # (13 + sqrt(5))/4
2340 >>> next(pqa) # (P_0, Q_0, a_0, A_0, B_0, G_0)
2341 (13, 4, 3, 3, 1, -1)
2342 >>> next(pqa) # (P_1, Q_1, a_1, A_1, B_1, G_1)
2343 (-1, 1, 1, 4, 1, 3)
2345 References
2346 ==========
2348 .. [1] Solving the generalized Pell equation x^2 - Dy^2 = N, John P.
2349 Robertson, July 31, 2004, Pages 4 - 8. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
2350 """
2351 A_i_2 = B_i_1 = 0
2352 A_i_1 = B_i_2 = 1
2354 G_i_2 = -P_0
2355 G_i_1 = Q_0
2357 P_i = P_0
2358 Q_i = Q_0
2360 while True:
2362 a_i = floor((P_i + sqrt(D))/Q_i)
2363 A_i = a_i*A_i_1 + A_i_2
2364 B_i = a_i*B_i_1 + B_i_2
2365 G_i = a_i*G_i_1 + G_i_2
2367 yield P_i, Q_i, a_i, A_i, B_i, G_i
2369 A_i_1, A_i_2 = A_i, A_i_1
2370 B_i_1, B_i_2 = B_i, B_i_1
2371 G_i_1, G_i_2 = G_i, G_i_1
2373 P_i = a_i*Q_i - P_i
2374 Q_i = (D - P_i**2)/Q_i
2377def diop_bf_DN(D, N, t=symbols("t", integer=True)):
2378 r"""
2379 Uses brute force to solve the equation, `x^2 - Dy^2 = N`.
2381 Explanation
2382 ===========
2384 Mainly concerned with the generalized Pell equation which is the case when
2385 `D > 0, D` is not a perfect square. For more information on the case refer
2386 [1]_. Let `(t, u)` be the minimal positive solution of the equation
2387 `x^2 - Dy^2 = 1`. Then this method requires
2388 `\sqrt{\\frac{\mid N \mid (t \pm 1)}{2D}}` to be small.
2390 Usage
2391 =====
2393 ``diop_bf_DN(D, N, t)``: ``D`` and ``N`` are coefficients in
2394 `x^2 - Dy^2 = N` and ``t`` is the parameter to be used in the solutions.
2396 Details
2397 =======
2399 ``D`` and ``N`` correspond to D and N in the equation.
2400 ``t`` is the parameter to be used in the solutions.
2402 Examples
2403 ========
2405 >>> from sympy.solvers.diophantine.diophantine import diop_bf_DN
2406 >>> diop_bf_DN(13, -4)
2407 [(3, 1), (-3, 1), (36, 10)]
2408 >>> diop_bf_DN(986, 1)
2409 [(49299, 1570)]
2411 See Also
2412 ========
2414 diop_DN()
2416 References
2417 ==========
2419 .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
2420 Robertson, July 31, 2004, Page 15. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
2421 """
2422 D = as_int(D)
2423 N = as_int(N)
2425 sol = []
2426 a = diop_DN(D, 1)
2427 u = a[0][0]
2429 if abs(N) == 1:
2430 return diop_DN(D, N)
2432 elif N > 1:
2433 L1 = 0
2434 L2 = integer_nthroot(int(N*(u - 1)/(2*D)), 2)[0] + 1
2436 elif N < -1:
2437 L1, _exact = integer_nthroot(-int(N/D), 2)
2438 if not _exact:
2439 L1 += 1
2440 L2 = integer_nthroot(-int(N*(u + 1)/(2*D)), 2)[0] + 1
2442 else: # N = 0
2443 if D < 0:
2444 return [(0, 0)]
2445 elif D == 0:
2446 return [(0, t)]
2447 else:
2448 sD, _exact = integer_nthroot(D, 2)
2449 if _exact:
2450 return [(sD*t, t), (-sD*t, t)]
2451 else:
2452 return [(0, 0)]
2455 for y in range(L1, L2):
2456 try:
2457 x, _exact = integer_nthroot(N + D*y**2, 2)
2458 except ValueError:
2459 _exact = False
2460 if _exact:
2461 sol.append((x, y))
2462 if not equivalent(x, y, -x, y, D, N):
2463 sol.append((-x, y))
2465 return sol
2468def equivalent(u, v, r, s, D, N):
2469 """
2470 Returns True if two solutions `(u, v)` and `(r, s)` of `x^2 - Dy^2 = N`
2471 belongs to the same equivalence class and False otherwise.
2473 Explanation
2474 ===========
2476 Two solutions `(u, v)` and `(r, s)` to the above equation fall to the same
2477 equivalence class iff both `(ur - Dvs)` and `(us - vr)` are divisible by
2478 `N`. See reference [1]_. No test is performed to test whether `(u, v)` and
2479 `(r, s)` are actually solutions to the equation. User should take care of
2480 this.
2482 Usage
2483 =====
2485 ``equivalent(u, v, r, s, D, N)``: `(u, v)` and `(r, s)` are two solutions
2486 of the equation `x^2 - Dy^2 = N` and all parameters involved are integers.
2488 Examples
2489 ========
2491 >>> from sympy.solvers.diophantine.diophantine import equivalent
2492 >>> equivalent(18, 5, -18, -5, 13, -1)
2493 True
2494 >>> equivalent(3, 1, -18, 393, 109, -4)
2495 False
2497 References
2498 ==========
2500 .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
2501 Robertson, July 31, 2004, Page 12. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
2503 """
2504 return divisible(u*r - D*v*s, N) and divisible(u*s - v*r, N)
2507def length(P, Q, D):
2508 r"""
2509 Returns the (length of aperiodic part + length of periodic part) of
2510 continued fraction representation of `\\frac{P + \sqrt{D}}{Q}`.
2512 It is important to remember that this does NOT return the length of the
2513 periodic part but the sum of the lengths of the two parts as mentioned
2514 above.
2516 Usage
2517 =====
2519 ``length(P, Q, D)``: ``P``, ``Q`` and ``D`` are integers corresponding to
2520 the continued fraction `\\frac{P + \sqrt{D}}{Q}`.
2522 Details
2523 =======
2525 ``P``, ``D`` and ``Q`` corresponds to P, D and Q in the continued fraction,
2526 `\\frac{P + \sqrt{D}}{Q}`.
2528 Examples
2529 ========
2531 >>> from sympy.solvers.diophantine.diophantine import length
2532 >>> length(-2, 4, 5) # (-2 + sqrt(5))/4
2533 3
2534 >>> length(-5, 4, 17) # (-5 + sqrt(17))/4
2535 4
2537 See Also
2538 ========
2539 sympy.ntheory.continued_fraction.continued_fraction_periodic
2540 """
2541 from sympy.ntheory.continued_fraction import continued_fraction_periodic
2542 v = continued_fraction_periodic(P, Q, D)
2543 if isinstance(v[-1], list):
2544 rpt = len(v[-1])
2545 nonrpt = len(v) - 1
2546 else:
2547 rpt = 0
2548 nonrpt = len(v)
2549 return rpt + nonrpt
2552def transformation_to_DN(eq):
2553 """
2554 This function transforms general quadratic,
2555 `ax^2 + bxy + cy^2 + dx + ey + f = 0`
2556 to more easy to deal with `X^2 - DY^2 = N` form.
2558 Explanation
2559 ===========
2561 This is used to solve the general quadratic equation by transforming it to
2562 the latter form. Refer to [1]_ for more detailed information on the
2563 transformation. This function returns a tuple (A, B) where A is a 2 X 2
2564 matrix and B is a 2 X 1 matrix such that,
2566 Transpose([x y]) = A * Transpose([X Y]) + B
2568 Usage
2569 =====
2571 ``transformation_to_DN(eq)``: where ``eq`` is the quadratic to be
2572 transformed.
2574 Examples
2575 ========
2577 >>> from sympy.abc import x, y
2578 >>> from sympy.solvers.diophantine.diophantine import transformation_to_DN
2579 >>> A, B = transformation_to_DN(x**2 - 3*x*y - y**2 - 2*y + 1)
2580 >>> A
2581 Matrix([
2582 [1/26, 3/26],
2583 [ 0, 1/13]])
2584 >>> B
2585 Matrix([
2586 [-6/13],
2587 [-4/13]])
2589 A, B returned are such that Transpose((x y)) = A * Transpose((X Y)) + B.
2590 Substituting these values for `x` and `y` and a bit of simplifying work
2591 will give an equation of the form `x^2 - Dy^2 = N`.
2593 >>> from sympy.abc import X, Y
2594 >>> from sympy import Matrix, simplify
2595 >>> u = (A*Matrix([X, Y]) + B)[0] # Transformation for x
2596 >>> u
2597 X/26 + 3*Y/26 - 6/13
2598 >>> v = (A*Matrix([X, Y]) + B)[1] # Transformation for y
2599 >>> v
2600 Y/13 - 4/13
2602 Next we will substitute these formulas for `x` and `y` and do
2603 ``simplify()``.
2605 >>> eq = simplify((x**2 - 3*x*y - y**2 - 2*y + 1).subs(zip((x, y), (u, v))))
2606 >>> eq
2607 X**2/676 - Y**2/52 + 17/13
2609 By multiplying the denominator appropriately, we can get a Pell equation
2610 in the standard form.
2612 >>> eq * 676
2613 X**2 - 13*Y**2 + 884
2615 If only the final equation is needed, ``find_DN()`` can be used.
2617 See Also
2618 ========
2620 find_DN()
2622 References
2623 ==========
2625 .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0,
2626 John P.Robertson, May 8, 2003, Page 7 - 11.
2627 https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
2628 """
2630 var, coeff, diop_type = classify_diop(eq, _dict=False)
2631 if diop_type == BinaryQuadratic.name:
2632 return _transformation_to_DN(var, coeff)
2635def _transformation_to_DN(var, coeff):
2637 x, y = var
2639 a = coeff[x**2]
2640 b = coeff[x*y]
2641 c = coeff[y**2]
2642 d = coeff[x]
2643 e = coeff[y]
2644 f = coeff[1]
2646 a, b, c, d, e, f = [as_int(i) for i in _remove_gcd(a, b, c, d, e, f)]
2648 X, Y = symbols("X, Y", integer=True)
2650 if b:
2651 B, C = _rational_pq(2*a, b)
2652 A, T = _rational_pq(a, B**2)
2654 # eq_1 = A*B*X**2 + B*(c*T - A*C**2)*Y**2 + d*T*X + (B*e*T - d*T*C)*Y + f*T*B
2655 coeff = {X**2: A*B, X*Y: 0, Y**2: B*(c*T - A*C**2), X: d*T, Y: B*e*T - d*T*C, 1: f*T*B}
2656 A_0, B_0 = _transformation_to_DN([X, Y], coeff)
2657 return Matrix(2, 2, [S.One/B, -S(C)/B, 0, 1])*A_0, Matrix(2, 2, [S.One/B, -S(C)/B, 0, 1])*B_0
2659 else:
2660 if d:
2661 B, C = _rational_pq(2*a, d)
2662 A, T = _rational_pq(a, B**2)
2664 # eq_2 = A*X**2 + c*T*Y**2 + e*T*Y + f*T - A*C**2
2665 coeff = {X**2: A, X*Y: 0, Y**2: c*T, X: 0, Y: e*T, 1: f*T - A*C**2}
2666 A_0, B_0 = _transformation_to_DN([X, Y], coeff)
2667 return Matrix(2, 2, [S.One/B, 0, 0, 1])*A_0, Matrix(2, 2, [S.One/B, 0, 0, 1])*B_0 + Matrix([-S(C)/B, 0])
2669 else:
2670 if e:
2671 B, C = _rational_pq(2*c, e)
2672 A, T = _rational_pq(c, B**2)
2674 # eq_3 = a*T*X**2 + A*Y**2 + f*T - A*C**2
2675 coeff = {X**2: a*T, X*Y: 0, Y**2: A, X: 0, Y: 0, 1: f*T - A*C**2}
2676 A_0, B_0 = _transformation_to_DN([X, Y], coeff)
2677 return Matrix(2, 2, [1, 0, 0, S.One/B])*A_0, Matrix(2, 2, [1, 0, 0, S.One/B])*B_0 + Matrix([0, -S(C)/B])
2679 else:
2680 # TODO: pre-simplification: Not necessary but may simplify
2681 # the equation.
2683 return Matrix(2, 2, [S.One/a, 0, 0, 1]), Matrix([0, 0])
2686def find_DN(eq):
2687 """
2688 This function returns a tuple, `(D, N)` of the simplified form,
2689 `x^2 - Dy^2 = N`, corresponding to the general quadratic,
2690 `ax^2 + bxy + cy^2 + dx + ey + f = 0`.
2692 Solving the general quadratic is then equivalent to solving the equation
2693 `X^2 - DY^2 = N` and transforming the solutions by using the transformation
2694 matrices returned by ``transformation_to_DN()``.
2696 Usage
2697 =====
2699 ``find_DN(eq)``: where ``eq`` is the quadratic to be transformed.
2701 Examples
2702 ========
2704 >>> from sympy.abc import x, y
2705 >>> from sympy.solvers.diophantine.diophantine import find_DN
2706 >>> find_DN(x**2 - 3*x*y - y**2 - 2*y + 1)
2707 (13, -884)
2709 Interpretation of the output is that we get `X^2 -13Y^2 = -884` after
2710 transforming `x^2 - 3xy - y^2 - 2y + 1` using the transformation returned
2711 by ``transformation_to_DN()``.
2713 See Also
2714 ========
2716 transformation_to_DN()
2718 References
2719 ==========
2721 .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0,
2722 John P.Robertson, May 8, 2003, Page 7 - 11.
2723 https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
2724 """
2725 var, coeff, diop_type = classify_diop(eq, _dict=False)
2726 if diop_type == BinaryQuadratic.name:
2727 return _find_DN(var, coeff)
2730def _find_DN(var, coeff):
2732 x, y = var
2733 X, Y = symbols("X, Y", integer=True)
2734 A, B = _transformation_to_DN(var, coeff)
2736 u = (A*Matrix([X, Y]) + B)[0]
2737 v = (A*Matrix([X, Y]) + B)[1]
2738 eq = x**2*coeff[x**2] + x*y*coeff[x*y] + y**2*coeff[y**2] + x*coeff[x] + y*coeff[y] + coeff[1]
2740 simplified = _mexpand(eq.subs(zip((x, y), (u, v))))
2742 coeff = simplified.as_coefficients_dict()
2744 return -coeff[Y**2]/coeff[X**2], -coeff[1]/coeff[X**2]
2747def check_param(x, y, a, params):
2748 """
2749 If there is a number modulo ``a`` such that ``x`` and ``y`` are both
2750 integers, then return a parametric representation for ``x`` and ``y``
2751 else return (None, None).
2753 Here ``x`` and ``y`` are functions of ``t``.
2754 """
2755 from sympy.simplify.simplify import clear_coefficients
2757 if x.is_number and not x.is_Integer:
2758 return DiophantineSolutionSet([x, y], parameters=params)
2760 if y.is_number and not y.is_Integer:
2761 return DiophantineSolutionSet([x, y], parameters=params)
2763 m, n = symbols("m, n", integer=True)
2764 c, p = (m*x + n*y).as_content_primitive()
2765 if a % c.q:
2766 return DiophantineSolutionSet([x, y], parameters=params)
2768 # clear_coefficients(mx + b, R)[1] -> (R - b)/m
2769 eq = clear_coefficients(x, m)[1] - clear_coefficients(y, n)[1]
2770 junk, eq = eq.as_content_primitive()
2772 return _diop_solve(eq, params=params)
2775def diop_ternary_quadratic(eq, parameterize=False):
2776 """
2777 Solves the general quadratic ternary form,
2778 `ax^2 + by^2 + cz^2 + fxy + gyz + hxz = 0`.
2780 Returns a tuple `(x, y, z)` which is a base solution for the above
2781 equation. If there are no solutions, `(None, None, None)` is returned.
2783 Usage
2784 =====
2786 ``diop_ternary_quadratic(eq)``: Return a tuple containing a basic solution
2787 to ``eq``.
2789 Details
2790 =======
2792 ``eq`` should be an homogeneous expression of degree two in three variables
2793 and it is assumed to be zero.
2795 Examples
2796 ========
2798 >>> from sympy.abc import x, y, z
2799 >>> from sympy.solvers.diophantine.diophantine import diop_ternary_quadratic
2800 >>> diop_ternary_quadratic(x**2 + 3*y**2 - z**2)
2801 (1, 0, 1)
2802 >>> diop_ternary_quadratic(4*x**2 + 5*y**2 - z**2)
2803 (1, 0, 2)
2804 >>> diop_ternary_quadratic(45*x**2 - 7*y**2 - 8*x*y - z**2)
2805 (28, 45, 105)
2806 >>> diop_ternary_quadratic(x**2 - 49*y**2 - z**2 + 13*z*y -8*x*y)
2807 (9, 1, 5)
2808 """
2809 var, coeff, diop_type = classify_diop(eq, _dict=False)
2811 if diop_type in (
2812 HomogeneousTernaryQuadratic.name,
2813 HomogeneousTernaryQuadraticNormal.name):
2814 sol = _diop_ternary_quadratic(var, coeff)
2815 if len(sol) > 0:
2816 x_0, y_0, z_0 = list(sol)[0]
2817 else:
2818 x_0, y_0, z_0 = None, None, None
2820 if parameterize:
2821 return _parametrize_ternary_quadratic(
2822 (x_0, y_0, z_0), var, coeff)
2823 return x_0, y_0, z_0
2826def _diop_ternary_quadratic(_var, coeff):
2827 eq = sum([i*coeff[i] for i in coeff])
2828 if HomogeneousTernaryQuadratic(eq).matches():
2829 return HomogeneousTernaryQuadratic(eq, free_symbols=_var).solve()
2830 elif HomogeneousTernaryQuadraticNormal(eq).matches():
2831 return HomogeneousTernaryQuadraticNormal(eq, free_symbols=_var).solve()
2834def transformation_to_normal(eq):
2835 """
2836 Returns the transformation Matrix that converts a general ternary
2837 quadratic equation ``eq`` (`ax^2 + by^2 + cz^2 + dxy + eyz + fxz`)
2838 to a form without cross terms: `ax^2 + by^2 + cz^2 = 0`. This is
2839 not used in solving ternary quadratics; it is only implemented for
2840 the sake of completeness.
2841 """
2842 var, coeff, diop_type = classify_diop(eq, _dict=False)
2844 if diop_type in (
2845 "homogeneous_ternary_quadratic",
2846 "homogeneous_ternary_quadratic_normal"):
2847 return _transformation_to_normal(var, coeff)
2850def _transformation_to_normal(var, coeff):
2852 _var = list(var) # copy
2853 x, y, z = var
2855 if not any(coeff[i**2] for i in var):
2856 # https://math.stackexchange.com/questions/448051/transform-quadratic-ternary-form-to-normal-form/448065#448065
2857 a = coeff[x*y]
2858 b = coeff[y*z]
2859 c = coeff[x*z]
2860 swap = False
2861 if not a: # b can't be 0 or else there aren't 3 vars
2862 swap = True
2863 a, b = b, a
2864 T = Matrix(((1, 1, -b/a), (1, -1, -c/a), (0, 0, 1)))
2865 if swap:
2866 T.row_swap(0, 1)
2867 T.col_swap(0, 1)
2868 return T
2870 if coeff[x**2] == 0:
2871 # If the coefficient of x is zero change the variables
2872 if coeff[y**2] == 0:
2873 _var[0], _var[2] = var[2], var[0]
2874 T = _transformation_to_normal(_var, coeff)
2875 T.row_swap(0, 2)
2876 T.col_swap(0, 2)
2877 return T
2879 else:
2880 _var[0], _var[1] = var[1], var[0]
2881 T = _transformation_to_normal(_var, coeff)
2882 T.row_swap(0, 1)
2883 T.col_swap(0, 1)
2884 return T
2886 # Apply the transformation x --> X - (B*Y + C*Z)/(2*A)
2887 if coeff[x*y] != 0 or coeff[x*z] != 0:
2888 A = coeff[x**2]
2889 B = coeff[x*y]
2890 C = coeff[x*z]
2891 D = coeff[y**2]
2892 E = coeff[y*z]
2893 F = coeff[z**2]
2895 _coeff = {}
2897 _coeff[x**2] = 4*A**2
2898 _coeff[y**2] = 4*A*D - B**2
2899 _coeff[z**2] = 4*A*F - C**2
2900 _coeff[y*z] = 4*A*E - 2*B*C
2901 _coeff[x*y] = 0
2902 _coeff[x*z] = 0
2904 T_0 = _transformation_to_normal(_var, _coeff)
2905 return Matrix(3, 3, [1, S(-B)/(2*A), S(-C)/(2*A), 0, 1, 0, 0, 0, 1])*T_0
2907 elif coeff[y*z] != 0:
2908 if coeff[y**2] == 0:
2909 if coeff[z**2] == 0:
2910 # Equations of the form A*x**2 + E*yz = 0.
2911 # Apply transformation y -> Y + Z ans z -> Y - Z
2912 return Matrix(3, 3, [1, 0, 0, 0, 1, 1, 0, 1, -1])
2914 else:
2915 # Ax**2 + E*y*z + F*z**2 = 0
2916 _var[0], _var[2] = var[2], var[0]
2917 T = _transformation_to_normal(_var, coeff)
2918 T.row_swap(0, 2)
2919 T.col_swap(0, 2)
2920 return T
2922 else:
2923 # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, F may be zero
2924 _var[0], _var[1] = var[1], var[0]
2925 T = _transformation_to_normal(_var, coeff)
2926 T.row_swap(0, 1)
2927 T.col_swap(0, 1)
2928 return T
2930 else:
2931 return Matrix.eye(3)
2934def parametrize_ternary_quadratic(eq):
2935 """
2936 Returns the parametrized general solution for the ternary quadratic
2937 equation ``eq`` which has the form
2938 `ax^2 + by^2 + cz^2 + fxy + gyz + hxz = 0`.
2940 Examples
2941 ========
2943 >>> from sympy import Tuple, ordered
2944 >>> from sympy.abc import x, y, z
2945 >>> from sympy.solvers.diophantine.diophantine import parametrize_ternary_quadratic
2947 The parametrized solution may be returned with three parameters:
2949 >>> parametrize_ternary_quadratic(2*x**2 + y**2 - 2*z**2)
2950 (p**2 - 2*q**2, -2*p**2 + 4*p*q - 4*p*r - 4*q**2, p**2 - 4*p*q + 2*q**2 - 4*q*r)
2952 There might also be only two parameters:
2954 >>> parametrize_ternary_quadratic(4*x**2 + 2*y**2 - 3*z**2)
2955 (2*p**2 - 3*q**2, -4*p**2 + 12*p*q - 6*q**2, 4*p**2 - 8*p*q + 6*q**2)
2957 Notes
2958 =====
2960 Consider ``p`` and ``q`` in the previous 2-parameter
2961 solution and observe that more than one solution can be represented
2962 by a given pair of parameters. If `p` and ``q`` are not coprime, this is
2963 trivially true since the common factor will also be a common factor of the
2964 solution values. But it may also be true even when ``p`` and
2965 ``q`` are coprime:
2967 >>> sol = Tuple(*_)
2968 >>> p, q = ordered(sol.free_symbols)
2969 >>> sol.subs([(p, 3), (q, 2)])
2970 (6, 12, 12)
2971 >>> sol.subs([(q, 1), (p, 1)])
2972 (-1, 2, 2)
2973 >>> sol.subs([(q, 0), (p, 1)])
2974 (2, -4, 4)
2975 >>> sol.subs([(q, 1), (p, 0)])
2976 (-3, -6, 6)
2978 Except for sign and a common factor, these are equivalent to
2979 the solution of (1, 2, 2).
2981 References
2982 ==========
2984 .. [1] The algorithmic resolution of Diophantine equations, Nigel P. Smart,
2985 London Mathematical Society Student Texts 41, Cambridge University
2986 Press, Cambridge, 1998.
2988 """
2989 var, coeff, diop_type = classify_diop(eq, _dict=False)
2991 if diop_type in (
2992 "homogeneous_ternary_quadratic",
2993 "homogeneous_ternary_quadratic_normal"):
2994 x_0, y_0, z_0 = list(_diop_ternary_quadratic(var, coeff))[0]
2995 return _parametrize_ternary_quadratic(
2996 (x_0, y_0, z_0), var, coeff)
2999def _parametrize_ternary_quadratic(solution, _var, coeff):
3000 # called for a*x**2 + b*y**2 + c*z**2 + d*x*y + e*y*z + f*x*z = 0
3001 assert 1 not in coeff
3003 x_0, y_0, z_0 = solution
3005 v = list(_var) # copy
3007 if x_0 is None:
3008 return (None, None, None)
3010 if solution.count(0) >= 2:
3011 # if there are 2 zeros the equation reduces
3012 # to k*X**2 == 0 where X is x, y, or z so X must
3013 # be zero, too. So there is only the trivial
3014 # solution.
3015 return (None, None, None)
3017 if x_0 == 0:
3018 v[0], v[1] = v[1], v[0]
3019 y_p, x_p, z_p = _parametrize_ternary_quadratic(
3020 (y_0, x_0, z_0), v, coeff)
3021 return x_p, y_p, z_p
3023 x, y, z = v
3024 r, p, q = symbols("r, p, q", integer=True)
3026 eq = sum(k*v for k, v in coeff.items())
3027 eq_1 = _mexpand(eq.subs(zip(
3028 (x, y, z), (r*x_0, r*y_0 + p, r*z_0 + q))))
3029 A, B = eq_1.as_independent(r, as_Add=True)
3032 x = A*x_0
3033 y = (A*y_0 - _mexpand(B/r*p))
3034 z = (A*z_0 - _mexpand(B/r*q))
3036 return _remove_gcd(x, y, z)
3039def diop_ternary_quadratic_normal(eq, parameterize=False):
3040 """
3041 Solves the quadratic ternary diophantine equation,
3042 `ax^2 + by^2 + cz^2 = 0`.
3044 Explanation
3045 ===========
3047 Here the coefficients `a`, `b`, and `c` should be non zero. Otherwise the
3048 equation will be a quadratic binary or univariate equation. If solvable,
3049 returns a tuple `(x, y, z)` that satisfies the given equation. If the
3050 equation does not have integer solutions, `(None, None, None)` is returned.
3052 Usage
3053 =====
3055 ``diop_ternary_quadratic_normal(eq)``: where ``eq`` is an equation of the form
3056 `ax^2 + by^2 + cz^2 = 0`.
3058 Examples
3059 ========
3061 >>> from sympy.abc import x, y, z
3062 >>> from sympy.solvers.diophantine.diophantine import diop_ternary_quadratic_normal
3063 >>> diop_ternary_quadratic_normal(x**2 + 3*y**2 - z**2)
3064 (1, 0, 1)
3065 >>> diop_ternary_quadratic_normal(4*x**2 + 5*y**2 - z**2)
3066 (1, 0, 2)
3067 >>> diop_ternary_quadratic_normal(34*x**2 - 3*y**2 - 301*z**2)
3068 (4, 9, 1)
3069 """
3070 var, coeff, diop_type = classify_diop(eq, _dict=False)
3071 if diop_type == HomogeneousTernaryQuadraticNormal.name:
3072 sol = _diop_ternary_quadratic_normal(var, coeff)
3073 if len(sol) > 0:
3074 x_0, y_0, z_0 = list(sol)[0]
3075 else:
3076 x_0, y_0, z_0 = None, None, None
3077 if parameterize:
3078 return _parametrize_ternary_quadratic(
3079 (x_0, y_0, z_0), var, coeff)
3080 return x_0, y_0, z_0
3083def _diop_ternary_quadratic_normal(var, coeff):
3084 eq = sum([i * coeff[i] for i in coeff])
3085 return HomogeneousTernaryQuadraticNormal(eq, free_symbols=var).solve()
3088def sqf_normal(a, b, c, steps=False):
3089 """
3090 Return `a', b', c'`, the coefficients of the square-free normal
3091 form of `ax^2 + by^2 + cz^2 = 0`, where `a', b', c'` are pairwise
3092 prime. If `steps` is True then also return three tuples:
3093 `sq`, `sqf`, and `(a', b', c')` where `sq` contains the square
3094 factors of `a`, `b` and `c` after removing the `gcd(a, b, c)`;
3095 `sqf` contains the values of `a`, `b` and `c` after removing
3096 both the `gcd(a, b, c)` and the square factors.
3098 The solutions for `ax^2 + by^2 + cz^2 = 0` can be
3099 recovered from the solutions of `a'x^2 + b'y^2 + c'z^2 = 0`.
3101 Examples
3102 ========
3104 >>> from sympy.solvers.diophantine.diophantine import sqf_normal
3105 >>> sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11)
3106 (11, 1, 5)
3107 >>> sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11, True)
3108 ((3, 1, 7), (5, 55, 11), (11, 1, 5))
3110 References
3111 ==========
3113 .. [1] Legendre's Theorem, Legrange's Descent,
3114 https://public.csusm.edu/aitken_html/notes/legendre.pdf
3117 See Also
3118 ========
3120 reconstruct()
3121 """
3122 ABC = _remove_gcd(a, b, c)
3123 sq = tuple(square_factor(i) for i in ABC)
3124 sqf = A, B, C = tuple([i//j**2 for i,j in zip(ABC, sq)])
3125 pc = igcd(A, B)
3126 A /= pc
3127 B /= pc
3128 pa = igcd(B, C)
3129 B /= pa
3130 C /= pa
3131 pb = igcd(A, C)
3132 A /= pb
3133 B /= pb
3135 A *= pa
3136 B *= pb
3137 C *= pc
3139 if steps:
3140 return (sq, sqf, (A, B, C))
3141 else:
3142 return A, B, C
3145def square_factor(a):
3146 r"""
3147 Returns an integer `c` s.t. `a = c^2k, \ c,k \in Z`. Here `k` is square
3148 free. `a` can be given as an integer or a dictionary of factors.
3150 Examples
3151 ========
3153 >>> from sympy.solvers.diophantine.diophantine import square_factor
3154 >>> square_factor(24)
3155 2
3156 >>> square_factor(-36*3)
3157 6
3158 >>> square_factor(1)
3159 1
3160 >>> square_factor({3: 2, 2: 1, -1: 1}) # -18
3161 3
3163 See Also
3164 ========
3165 sympy.ntheory.factor_.core
3166 """
3167 f = a if isinstance(a, dict) else factorint(a)
3168 return Mul(*[p**(e//2) for p, e in f.items()])
3171def reconstruct(A, B, z):
3172 """
3173 Reconstruct the `z` value of an equivalent solution of `ax^2 + by^2 + cz^2`
3174 from the `z` value of a solution of the square-free normal form of the
3175 equation, `a'*x^2 + b'*y^2 + c'*z^2`, where `a'`, `b'` and `c'` are square
3176 free and `gcd(a', b', c') == 1`.
3177 """
3178 f = factorint(igcd(A, B))
3179 for p, e in f.items():
3180 if e != 1:
3181 raise ValueError('a and b should be square-free')
3182 z *= p
3183 return z
3186def ldescent(A, B):
3187 """
3188 Return a non-trivial solution to `w^2 = Ax^2 + By^2` using
3189 Lagrange's method; return None if there is no such solution.
3190 .
3192 Here, `A \\neq 0` and `B \\neq 0` and `A` and `B` are square free. Output a
3193 tuple `(w_0, x_0, y_0)` which is a solution to the above equation.
3195 Examples
3196 ========
3198 >>> from sympy.solvers.diophantine.diophantine import ldescent
3199 >>> ldescent(1, 1) # w^2 = x^2 + y^2
3200 (1, 1, 0)
3201 >>> ldescent(4, -7) # w^2 = 4x^2 - 7y^2
3202 (2, -1, 0)
3204 This means that `x = -1, y = 0` and `w = 2` is a solution to the equation
3205 `w^2 = 4x^2 - 7y^2`
3207 >>> ldescent(5, -1) # w^2 = 5x^2 - y^2
3208 (2, 1, -1)
3210 References
3211 ==========
3213 .. [1] The algorithmic resolution of Diophantine equations, Nigel P. Smart,
3214 London Mathematical Society Student Texts 41, Cambridge University
3215 Press, Cambridge, 1998.
3216 .. [2] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
3217 [online], Available:
3218 https://nottingham-repository.worktribe.com/output/1023265/efficient-solution-of-rational-conics
3219 """
3220 if abs(A) > abs(B):
3221 w, y, x = ldescent(B, A)
3222 return w, x, y
3224 if A == 1:
3225 return (1, 1, 0)
3227 if B == 1:
3228 return (1, 0, 1)
3230 if B == -1: # and A == -1
3231 return
3233 r = sqrt_mod(A, B)
3235 Q = (r**2 - A) // B
3237 if Q == 0:
3238 B_0 = 1
3239 d = 0
3240 else:
3241 div = divisors(Q)
3242 B_0 = None
3244 for i in div:
3245 sQ, _exact = integer_nthroot(abs(Q) // i, 2)
3246 if _exact:
3247 B_0, d = sign(Q)*i, sQ
3248 break
3250 if B_0 is not None:
3251 W, X, Y = ldescent(A, B_0)
3252 return _remove_gcd((-A*X + r*W), (r*X - W), Y*(B_0*d))
3255def descent(A, B):
3256 """
3257 Returns a non-trivial solution, (x, y, z), to `x^2 = Ay^2 + Bz^2`
3258 using Lagrange's descent method with lattice-reduction. `A` and `B`
3259 are assumed to be valid for such a solution to exist.
3261 This is faster than the normal Lagrange's descent algorithm because
3262 the Gaussian reduction is used.
3264 Examples
3265 ========
3267 >>> from sympy.solvers.diophantine.diophantine import descent
3268 >>> descent(3, 1) # x**2 = 3*y**2 + z**2
3269 (1, 0, 1)
3271 `(x, y, z) = (1, 0, 1)` is a solution to the above equation.
3273 >>> descent(41, -113)
3274 (-16, -3, 1)
3276 References
3277 ==========
3279 .. [1] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
3280 Mathematics of Computation, Volume 00, Number 0.
3281 """
3282 if abs(A) > abs(B):
3283 x, y, z = descent(B, A)
3284 return x, z, y
3286 if B == 1:
3287 return (1, 0, 1)
3288 if A == 1:
3289 return (1, 1, 0)
3290 if B == -A:
3291 return (0, 1, 1)
3292 if B == A:
3293 x, z, y = descent(-1, A)
3294 return (A*y, z, x)
3296 w = sqrt_mod(A, B)
3297 x_0, z_0 = gaussian_reduce(w, A, B)
3299 t = (x_0**2 - A*z_0**2) // B
3300 t_2 = square_factor(t)
3301 t_1 = t // t_2**2
3303 x_1, z_1, y_1 = descent(A, t_1)
3305 return _remove_gcd(x_0*x_1 + A*z_0*z_1, z_0*x_1 + x_0*z_1, t_1*t_2*y_1)
3308def gaussian_reduce(w, a, b):
3309 r"""
3310 Returns a reduced solution `(x, z)` to the congruence
3311 `X^2 - aZ^2 \equiv 0 \ (mod \ b)` so that `x^2 + |a|z^2` is minimal.
3313 Details
3314 =======
3316 Here ``w`` is a solution of the congruence `x^2 \equiv a \ (mod \ b)`
3318 References
3319 ==========
3321 .. [1] Gaussian lattice Reduction [online]. Available:
3322 https://web.archive.org/web/20201021115213/http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=404
3323 .. [2] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
3324 Mathematics of Computation, Volume 00, Number 0.
3325 """
3326 u = (0, 1)
3327 v = (1, 0)
3329 if dot(u, v, w, a, b) < 0:
3330 v = (-v[0], -v[1])
3332 if norm(u, w, a, b) < norm(v, w, a, b):
3333 u, v = v, u
3335 while norm(u, w, a, b) > norm(v, w, a, b):
3336 k = dot(u, v, w, a, b) // dot(v, v, w, a, b)
3337 u, v = v, (u[0]- k*v[0], u[1]- k*v[1])
3339 u, v = v, u
3341 if dot(u, v, w, a, b) < dot(v, v, w, a, b)/2 or norm((u[0]-v[0], u[1]-v[1]), w, a, b) > norm(v, w, a, b):
3342 c = v
3343 else:
3344 c = (u[0] - v[0], u[1] - v[1])
3346 return c[0]*w + b*c[1], c[0]
3349def dot(u, v, w, a, b):
3350 r"""
3351 Returns a special dot product of the vectors `u = (u_{1}, u_{2})` and
3352 `v = (v_{1}, v_{2})` which is defined in order to reduce solution of
3353 the congruence equation `X^2 - aZ^2 \equiv 0 \ (mod \ b)`.
3354 """
3355 u_1, u_2 = u
3356 v_1, v_2 = v
3357 return (w*u_1 + b*u_2)*(w*v_1 + b*v_2) + abs(a)*u_1*v_1
3360def norm(u, w, a, b):
3361 r"""
3362 Returns the norm of the vector `u = (u_{1}, u_{2})` under the dot product
3363 defined by `u \cdot v = (wu_{1} + bu_{2})(w*v_{1} + bv_{2}) + |a|*u_{1}*v_{1}`
3364 where `u = (u_{1}, u_{2})` and `v = (v_{1}, v_{2})`.
3365 """
3366 u_1, u_2 = u
3367 return sqrt(dot((u_1, u_2), (u_1, u_2), w, a, b))
3370def holzer(x, y, z, a, b, c):
3371 r"""
3372 Simplify the solution `(x, y, z)` of the equation
3373 `ax^2 + by^2 = cz^2` with `a, b, c > 0` and `z^2 \geq \mid ab \mid` to
3374 a new reduced solution `(x', y', z')` such that `z'^2 \leq \mid ab \mid`.
3376 The algorithm is an interpretation of Mordell's reduction as described
3377 on page 8 of Cremona and Rusin's paper [1]_ and the work of Mordell in
3378 reference [2]_.
3380 References
3381 ==========
3383 .. [1] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
3384 Mathematics of Computation, Volume 00, Number 0.
3385 .. [2] Diophantine Equations, L. J. Mordell, page 48.
3387 """
3389 if _odd(c):
3390 k = 2*c
3391 else:
3392 k = c//2
3394 small = a*b*c
3395 step = 0
3396 while True:
3397 t1, t2, t3 = a*x**2, b*y**2, c*z**2
3398 # check that it's a solution
3399 if t1 + t2 != t3:
3400 if step == 0:
3401 raise ValueError('bad starting solution')
3402 break
3403 x_0, y_0, z_0 = x, y, z
3404 if max(t1, t2, t3) <= small:
3405 # Holzer condition
3406 break
3408 uv = u, v = base_solution_linear(k, y_0, -x_0)
3409 if None in uv:
3410 break
3412 p, q = -(a*u*x_0 + b*v*y_0), c*z_0
3413 r = Rational(p, q)
3414 if _even(c):
3415 w = _nint_or_floor(p, q)
3416 assert abs(w - r) <= S.Half
3417 else:
3418 w = p//q # floor
3419 if _odd(a*u + b*v + c*w):
3420 w += 1
3421 assert abs(w - r) <= S.One
3423 A = (a*u**2 + b*v**2 + c*w**2)
3424 B = (a*u*x_0 + b*v*y_0 + c*w*z_0)
3425 x = Rational(x_0*A - 2*u*B, k)
3426 y = Rational(y_0*A - 2*v*B, k)
3427 z = Rational(z_0*A - 2*w*B, k)
3428 assert all(i.is_Integer for i in (x, y, z))
3429 step += 1
3431 return tuple([int(i) for i in (x_0, y_0, z_0)])
3434def diop_general_pythagorean(eq, param=symbols("m", integer=True)):
3435 """
3436 Solves the general pythagorean equation,
3437 `a_{1}^2x_{1}^2 + a_{2}^2x_{2}^2 + . . . + a_{n}^2x_{n}^2 - a_{n + 1}^2x_{n + 1}^2 = 0`.
3439 Returns a tuple which contains a parametrized solution to the equation,
3440 sorted in the same order as the input variables.
3442 Usage
3443 =====
3445 ``diop_general_pythagorean(eq, param)``: where ``eq`` is a general
3446 pythagorean equation which is assumed to be zero and ``param`` is the base
3447 parameter used to construct other parameters by subscripting.
3449 Examples
3450 ========
3452 >>> from sympy.solvers.diophantine.diophantine import diop_general_pythagorean
3453 >>> from sympy.abc import a, b, c, d, e
3454 >>> diop_general_pythagorean(a**2 + b**2 + c**2 - d**2)
3455 (m1**2 + m2**2 - m3**2, 2*m1*m3, 2*m2*m3, m1**2 + m2**2 + m3**2)
3456 >>> diop_general_pythagorean(9*a**2 - 4*b**2 + 16*c**2 + 25*d**2 + e**2)
3457 (10*m1**2 + 10*m2**2 + 10*m3**2 - 10*m4**2, 15*m1**2 + 15*m2**2 + 15*m3**2 + 15*m4**2, 15*m1*m4, 12*m2*m4, 60*m3*m4)
3458 """
3459 var, coeff, diop_type = classify_diop(eq, _dict=False)
3461 if diop_type == GeneralPythagorean.name:
3462 if param is None:
3463 params = None
3464 else:
3465 params = symbols('%s1:%i' % (param, len(var)), integer=True)
3466 return list(GeneralPythagorean(eq).solve(parameters=params))[0]
3469def diop_general_sum_of_squares(eq, limit=1):
3470 r"""
3471 Solves the equation `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`.
3473 Returns at most ``limit`` number of solutions.
3475 Usage
3476 =====
3478 ``general_sum_of_squares(eq, limit)`` : Here ``eq`` is an expression which
3479 is assumed to be zero. Also, ``eq`` should be in the form,
3480 `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`.
3482 Details
3483 =======
3485 When `n = 3` if `k = 4^a(8m + 7)` for some `a, m \in Z` then there will be
3486 no solutions. Refer to [1]_ for more details.
3488 Examples
3489 ========
3491 >>> from sympy.solvers.diophantine.diophantine import diop_general_sum_of_squares
3492 >>> from sympy.abc import a, b, c, d, e
3493 >>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345)
3494 {(15, 22, 22, 24, 24)}
3496 Reference
3497 =========
3499 .. [1] Representing an integer as a sum of three squares, [online],
3500 Available:
3501 https://www.proofwiki.org/wiki/Integer_as_Sum_of_Three_Squares
3502 """
3503 var, coeff, diop_type = classify_diop(eq, _dict=False)
3505 if diop_type == GeneralSumOfSquares.name:
3506 return set(GeneralSumOfSquares(eq).solve(limit=limit))
3509def diop_general_sum_of_even_powers(eq, limit=1):
3510 """
3511 Solves the equation `x_{1}^e + x_{2}^e + . . . + x_{n}^e - k = 0`
3512 where `e` is an even, integer power.
3514 Returns at most ``limit`` number of solutions.
3516 Usage
3517 =====
3519 ``general_sum_of_even_powers(eq, limit)`` : Here ``eq`` is an expression which
3520 is assumed to be zero. Also, ``eq`` should be in the form,
3521 `x_{1}^e + x_{2}^e + . . . + x_{n}^e - k = 0`.
3523 Examples
3524 ========
3526 >>> from sympy.solvers.diophantine.diophantine import diop_general_sum_of_even_powers
3527 >>> from sympy.abc import a, b
3528 >>> diop_general_sum_of_even_powers(a**4 + b**4 - (2**4 + 3**4))
3529 {(2, 3)}
3531 See Also
3532 ========
3534 power_representation
3535 """
3536 var, coeff, diop_type = classify_diop(eq, _dict=False)
3538 if diop_type == GeneralSumOfEvenPowers.name:
3539 return set(GeneralSumOfEvenPowers(eq).solve(limit=limit))
3542## Functions below this comment can be more suitably grouped under
3543## an Additive number theory module rather than the Diophantine
3544## equation module.
3547def partition(n, k=None, zeros=False):
3548 """
3549 Returns a generator that can be used to generate partitions of an integer
3550 `n`.
3552 Explanation
3553 ===========
3555 A partition of `n` is a set of positive integers which add up to `n`. For
3556 example, partitions of 3 are 3, 1 + 2, 1 + 1 + 1. A partition is returned
3557 as a tuple. If ``k`` equals None, then all possible partitions are returned
3558 irrespective of their size, otherwise only the partitions of size ``k`` are
3559 returned. If the ``zero`` parameter is set to True then a suitable
3560 number of zeros are added at the end of every partition of size less than
3561 ``k``.
3563 ``zero`` parameter is considered only if ``k`` is not None. When the
3564 partitions are over, the last `next()` call throws the ``StopIteration``
3565 exception, so this function should always be used inside a try - except
3566 block.
3568 Details
3569 =======
3571 ``partition(n, k)``: Here ``n`` is a positive integer and ``k`` is the size
3572 of the partition which is also positive integer.
3574 Examples
3575 ========
3577 >>> from sympy.solvers.diophantine.diophantine import partition
3578 >>> f = partition(5)
3579 >>> next(f)
3580 (1, 1, 1, 1, 1)
3581 >>> next(f)
3582 (1, 1, 1, 2)
3583 >>> g = partition(5, 3)
3584 >>> next(g)
3585 (1, 1, 3)
3586 >>> next(g)
3587 (1, 2, 2)
3588 >>> g = partition(5, 3, zeros=True)
3589 >>> next(g)
3590 (0, 0, 5)
3592 """
3593 if not zeros or k is None:
3594 for i in ordered_partitions(n, k):
3595 yield tuple(i)
3596 else:
3597 for m in range(1, k + 1):
3598 for i in ordered_partitions(n, m):
3599 i = tuple(i)
3600 yield (0,)*(k - len(i)) + i
3603def prime_as_sum_of_two_squares(p):
3604 """
3605 Represent a prime `p` as a unique sum of two squares; this can
3606 only be done if the prime is congruent to 1 mod 4.
3608 Examples
3609 ========
3611 >>> from sympy.solvers.diophantine.diophantine import prime_as_sum_of_two_squares
3612 >>> prime_as_sum_of_two_squares(7) # can't be done
3613 >>> prime_as_sum_of_two_squares(5)
3614 (1, 2)
3616 Reference
3617 =========
3619 .. [1] Representing a number as a sum of four squares, [online],
3620 Available: https://schorn.ch/lagrange.html
3622 See Also
3623 ========
3624 sum_of_squares()
3625 """
3626 if not p % 4 == 1:
3627 return
3629 if p % 8 == 5:
3630 b = 2
3631 else:
3632 b = 3
3634 while pow(b, (p - 1) // 2, p) == 1:
3635 b = nextprime(b)
3637 b = pow(b, (p - 1) // 4, p)
3638 a = p
3640 while b**2 > p:
3641 a, b = b, a % b
3643 return (int(a % b), int(b)) # convert from long
3646def sum_of_three_squares(n):
3647 r"""
3648 Returns a 3-tuple $(a, b, c)$ such that $a^2 + b^2 + c^2 = n$ and
3649 $a, b, c \geq 0$.
3651 Returns None if $n = 4^a(8m + 7)$ for some `a, m \in \mathbb{Z}`. See
3652 [1]_ for more details.
3654 Usage
3655 =====
3657 ``sum_of_three_squares(n)``: Here ``n`` is a non-negative integer.
3659 Examples
3660 ========
3662 >>> from sympy.solvers.diophantine.diophantine import sum_of_three_squares
3663 >>> sum_of_three_squares(44542)
3664 (18, 37, 207)
3666 References
3667 ==========
3669 .. [1] Representing a number as a sum of three squares, [online],
3670 Available: https://schorn.ch/lagrange.html
3672 See Also
3673 ========
3675 sum_of_squares()
3676 """
3677 special = {1:(1, 0, 0), 2:(1, 1, 0), 3:(1, 1, 1), 10: (1, 3, 0), 34: (3, 3, 4), 58:(3, 7, 0),
3678 85:(6, 7, 0), 130:(3, 11, 0), 214:(3, 6, 13), 226:(8, 9, 9), 370:(8, 9, 15),
3679 526:(6, 7, 21), 706:(15, 15, 16), 730:(1, 27, 0), 1414:(6, 17, 33), 1906:(13, 21, 36),
3680 2986: (21, 32, 39), 9634: (56, 57, 57)}
3682 v = 0
3684 if n == 0:
3685 return (0, 0, 0)
3687 v = multiplicity(4, n)
3688 n //= 4**v
3690 if n % 8 == 7:
3691 return
3693 if n in special.keys():
3694 x, y, z = special[n]
3695 return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
3697 s, _exact = integer_nthroot(n, 2)
3699 if _exact:
3700 return (2**v*s, 0, 0)
3702 x = None
3704 if n % 8 == 3:
3705 s = s if _odd(s) else s - 1
3707 for x in range(s, -1, -2):
3708 N = (n - x**2) // 2
3709 if isprime(N):
3710 y, z = prime_as_sum_of_two_squares(N)
3711 return _sorted_tuple(2**v*x, 2**v*(y + z), 2**v*abs(y - z))
3712 return
3714 if n % 8 in (2, 6):
3715 s = s if _odd(s) else s - 1
3716 else:
3717 s = s - 1 if _odd(s) else s
3719 for x in range(s, -1, -2):
3720 N = n - x**2
3721 if isprime(N):
3722 y, z = prime_as_sum_of_two_squares(N)
3723 return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
3726def sum_of_four_squares(n):
3727 r"""
3728 Returns a 4-tuple `(a, b, c, d)` such that `a^2 + b^2 + c^2 + d^2 = n`.
3730 Here `a, b, c, d \geq 0`.
3732 Usage
3733 =====
3735 ``sum_of_four_squares(n)``: Here ``n`` is a non-negative integer.
3737 Examples
3738 ========
3740 >>> from sympy.solvers.diophantine.diophantine import sum_of_four_squares
3741 >>> sum_of_four_squares(3456)
3742 (8, 8, 32, 48)
3743 >>> sum_of_four_squares(1294585930293)
3744 (0, 1234, 2161, 1137796)
3746 References
3747 ==========
3749 .. [1] Representing a number as a sum of four squares, [online],
3750 Available: https://schorn.ch/lagrange.html
3752 See Also
3753 ========
3755 sum_of_squares()
3756 """
3757 if n == 0:
3758 return (0, 0, 0, 0)
3760 v = multiplicity(4, n)
3761 n //= 4**v
3763 if n % 8 == 7:
3764 d = 2
3765 n = n - 4
3766 elif n % 8 in (2, 6):
3767 d = 1
3768 n = n - 1
3769 else:
3770 d = 0
3772 x, y, z = sum_of_three_squares(n)
3774 return _sorted_tuple(2**v*d, 2**v*x, 2**v*y, 2**v*z)
3777def power_representation(n, p, k, zeros=False):
3778 r"""
3779 Returns a generator for finding k-tuples of integers,
3780 `(n_{1}, n_{2}, . . . n_{k})`, such that
3781 `n = n_{1}^p + n_{2}^p + . . . n_{k}^p`.
3783 Usage
3784 =====
3786 ``power_representation(n, p, k, zeros)``: Represent non-negative number
3787 ``n`` as a sum of ``k`` ``p``\ th powers. If ``zeros`` is true, then the
3788 solutions is allowed to contain zeros.
3790 Examples
3791 ========
3793 >>> from sympy.solvers.diophantine.diophantine import power_representation
3795 Represent 1729 as a sum of two cubes:
3797 >>> f = power_representation(1729, 3, 2)
3798 >>> next(f)
3799 (9, 10)
3800 >>> next(f)
3801 (1, 12)
3803 If the flag `zeros` is True, the solution may contain tuples with
3804 zeros; any such solutions will be generated after the solutions
3805 without zeros:
3807 >>> list(power_representation(125, 2, 3, zeros=True))
3808 [(5, 6, 8), (3, 4, 10), (0, 5, 10), (0, 2, 11)]
3810 For even `p` the `permute_sign` function can be used to get all
3811 signed values:
3813 >>> from sympy.utilities.iterables import permute_signs
3814 >>> list(permute_signs((1, 12)))
3815 [(1, 12), (-1, 12), (1, -12), (-1, -12)]
3817 All possible signed permutations can also be obtained:
3819 >>> from sympy.utilities.iterables import signed_permutations
3820 >>> list(signed_permutations((1, 12)))
3821 [(1, 12), (-1, 12), (1, -12), (-1, -12), (12, 1), (-12, 1), (12, -1), (-12, -1)]
3822 """
3823 n, p, k = [as_int(i) for i in (n, p, k)]
3825 if n < 0:
3826 if p % 2:
3827 for t in power_representation(-n, p, k, zeros):
3828 yield tuple(-i for i in t)
3829 return
3831 if p < 1 or k < 1:
3832 raise ValueError(filldedent('''
3833 Expecting positive integers for `(p, k)`, but got `(%s, %s)`'''
3834 % (p, k)))
3836 if n == 0:
3837 if zeros:
3838 yield (0,)*k
3839 return
3841 if k == 1:
3842 if p == 1:
3843 yield (n,)
3844 else:
3845 be = perfect_power(n)
3846 if be:
3847 b, e = be
3848 d, r = divmod(e, p)
3849 if not r:
3850 yield (b**d,)
3851 return
3853 if p == 1:
3854 for t in partition(n, k, zeros=zeros):
3855 yield t
3856 return
3858 if p == 2:
3859 feasible = _can_do_sum_of_squares(n, k)
3860 if not feasible:
3861 return
3862 if not zeros and n > 33 and k >= 5 and k <= n and n - k in (
3863 13, 10, 7, 5, 4, 2, 1):
3864 '''Todd G. Will, "When Is n^2 a Sum of k Squares?", [online].
3865 Available: https://www.maa.org/sites/default/files/Will-MMz-201037918.pdf'''
3866 return
3867 if feasible is not True: # it's prime and k == 2
3868 yield prime_as_sum_of_two_squares(n)
3869 return
3871 if k == 2 and p > 2:
3872 be = perfect_power(n)
3873 if be and be[1] % p == 0:
3874 return # Fermat: a**n + b**n = c**n has no solution for n > 2
3876 if n >= k:
3877 a = integer_nthroot(n - (k - 1), p)[0]
3878 for t in pow_rep_recursive(a, k, n, [], p):
3879 yield tuple(reversed(t))
3881 if zeros:
3882 a = integer_nthroot(n, p)[0]
3883 for i in range(1, k):
3884 for t in pow_rep_recursive(a, i, n, [], p):
3885 yield tuple(reversed(t + (0,)*(k - i)))
3888sum_of_powers = power_representation
3891def pow_rep_recursive(n_i, k, n_remaining, terms, p):
3892 # Invalid arguments
3893 if n_i <= 0 or k <= 0:
3894 return
3896 # No solutions may exist
3897 if n_remaining < k:
3898 return
3899 if k * pow(n_i, p) < n_remaining:
3900 return
3902 if k == 0 and n_remaining == 0:
3903 yield tuple(terms)
3905 elif k == 1:
3906 # next_term^p must equal to n_remaining
3907 next_term, exact = integer_nthroot(n_remaining, p)
3908 if exact and next_term <= n_i:
3909 yield tuple(terms + [next_term])
3910 return
3912 else:
3913 # TODO: Fall back to diop_DN when k = 2
3914 if n_i >= 1 and k > 0:
3915 for next_term in range(1, n_i + 1):
3916 residual = n_remaining - pow(next_term, p)
3917 if residual < 0:
3918 break
3919 yield from pow_rep_recursive(next_term, k - 1, residual, terms + [next_term], p)
3922def sum_of_squares(n, k, zeros=False):
3923 """Return a generator that yields the k-tuples of nonnegative
3924 values, the squares of which sum to n. If zeros is False (default)
3925 then the solution will not contain zeros. The nonnegative
3926 elements of a tuple are sorted.
3928 * If k == 1 and n is square, (n,) is returned.
3930 * If k == 2 then n can only be written as a sum of squares if
3931 every prime in the factorization of n that has the form
3932 4*k + 3 has an even multiplicity. If n is prime then
3933 it can only be written as a sum of two squares if it is
3934 in the form 4*k + 1.
3936 * if k == 3 then n can be written as a sum of squares if it does
3937 not have the form 4**m*(8*k + 7).
3939 * all integers can be written as the sum of 4 squares.
3941 * if k > 4 then n can be partitioned and each partition can
3942 be written as a sum of 4 squares; if n is not evenly divisible
3943 by 4 then n can be written as a sum of squares only if the
3944 an additional partition can be written as sum of squares.
3945 For example, if k = 6 then n is partitioned into two parts,
3946 the first being written as a sum of 4 squares and the second
3947 being written as a sum of 2 squares -- which can only be
3948 done if the condition above for k = 2 can be met, so this will
3949 automatically reject certain partitions of n.
3951 Examples
3952 ========
3954 >>> from sympy.solvers.diophantine.diophantine import sum_of_squares
3955 >>> list(sum_of_squares(25, 2))
3956 [(3, 4)]
3957 >>> list(sum_of_squares(25, 2, True))
3958 [(3, 4), (0, 5)]
3959 >>> list(sum_of_squares(25, 4))
3960 [(1, 2, 2, 4)]
3962 See Also
3963 ========
3965 sympy.utilities.iterables.signed_permutations
3966 """
3967 yield from power_representation(n, 2, k, zeros)
3970def _can_do_sum_of_squares(n, k):
3971 """Return True if n can be written as the sum of k squares,
3972 False if it cannot, or 1 if ``k == 2`` and ``n`` is prime (in which
3973 case it *can* be written as a sum of two squares). A False
3974 is returned only if it cannot be written as ``k``-squares, even
3975 if 0s are allowed.
3976 """
3977 if k < 1:
3978 return False
3979 if n < 0:
3980 return False
3981 if n == 0:
3982 return True
3983 if k == 1:
3984 return is_square(n)
3985 if k == 2:
3986 if n in (1, 2):
3987 return True
3988 if isprime(n):
3989 if n % 4 == 1:
3990 return 1 # signal that it was prime
3991 return False
3992 else:
3993 f = factorint(n)
3994 for p, m in f.items():
3995 # we can proceed iff no prime factor in the form 4*k + 3
3996 # has an odd multiplicity
3997 if (p % 4 == 3) and m % 2:
3998 return False
3999 return True
4000 if k == 3:
4001 if (n//4**multiplicity(4, n)) % 8 == 7:
4002 return False
4003 # every number can be written as a sum of 4 squares; for k > 4 partitions
4004 # can be 0
4005 return True