0

I have a triplet, for example (1806336, 1849600, 93636) and I need to solve the diophantine equation:

x^2-2w^2+2y^2-z^2-constant=0

Here the constant is 1806336 + 1849600 - 93636.

What I tried is:

from sympy.solvers.diophantine.diophantine import diop_general_pythagorean
from sympy.abc import x, w, y, z
tuple = [1806336, 1849600, 93636]
s = tuple[0]
t = tuple[1]
u = tuple[2]
diop_general_pythagorean(x**2 - 2*w**2 + 2*y**2 - z**2 - s-t+u)

which yields the exception

"This equation is not yet recognized or else has not been simplified sufficiently to put it in a form recognized by diop_classify()."

Without success, I also tried to avoid using the constant as a variable and wrote directly:

diop_general_pythagorean(x**2 - 2*w**2 + 2*y**2 - z**2 - 3562300)

The error remains. It would be great to have a chance solving such a diophantine equation directly using a constant which I may express as a variable.

azro
  • 53,056
  • 7
  • 34
  • 70
Eldar Sultanow
  • 205
  • 2
  • 9

2 Answers2

3

Maybe Z3 is an option?

from z3 import Ints, solve

x, y, z, w = Ints('x y z w')
s, t, u = (1806336, 1849600, 93636)
solve(x**2 - 2*w**2 + 2*y**2 - z**2 - s-t+u==0)

Gives [x = -368, w = -865, z = 14, y = 1569] as one possibility.

JohanC
  • 71,591
  • 8
  • 33
  • 66
  • Note that I copied your formula unchanged, although the way you wrote it, you might want `... - (s-t+u)` instead of `... - s-t+u`. – JohanC Dec 28 '21 at 18:38
  • Thank you. With your help I even could optimize the task. I am now using `solve(x*x-w*w==s, y*y-w*w==t, z*z-y*y==u, w!=0)`. Here is the [code](https://github.com/Sultanow/pythagorean/blob/main/pythagorean_search2.py). Meanwhile I implemented it in Mathematica as well using `FindInstance[ x*x - w*w == s*s && y*y - w*w == t*t && w != 0, {w, x, y}, Integers]`, which yields the same results, but drastically faster. The background is described here on [MSE](https://math.stackexchange.com/questions/4342283/for-integers-xyz-why-are-these-cases-impossible-for-mengolis-six-square-pr/4343676#4343676) – Eldar Sultanow Dec 28 '21 at 20:37
  • So in the MSE system 2 formulation, is this equivalent to seeking a tetrahedron comprised of 3 sides that are right triangles with integer-length sides (Pythagorean triples)? – smichr Dec 29 '21 at 19:20
2

You already suggest a way to approach this problem in a way that SymPy can help. SymPy can solve x**2 - y**2 - c where c is a constant. Your problem can be thought of in 3 parts with that form:

x**2 - w**2 - s + y**2 + w**2 - t + y**2 - z**2 + u = 0

Here is a way it might be done, using a smaller constant, 28, which I will break up into 15 + 8 + 5. I am going to limit output to positive values since the negatives are redundantly true because of each value being squared.

>>> from sympy import diophantine
>>> from sympy.abc import x, y
>>> pos = lambda eq: [i for i in diophantine(eq) if all(j>0 for j in i)]

Since I broke 28 into 3 terms with the same sign, I will solve x**2 - y**2 - c for those 3 values of c and look for a solution that fits the pattern:

>>> pos(x**2-y**2-15)
[(8, 7), (4, 1)]
>>> pos(x**2-y**2-8)
[(3, 1)]
>>> pos(x**2-y**2-5)
[(3, 2)]

So for (x,w),(y,w),(y,z) the points (4,1),(3,1),(3,2) work. So look for a solution between sub-problems for s and t that have a matching 2nd element and then a solution for u that has a first element that match either of the first elements in the s and t candidates.

smichr
  • 16,948
  • 2
  • 27
  • 34