# How to reduce the time to solve this system of equations?

I am trying to solve the system of equations. My code

  {x, y, z, a, b, m, n} /. Solve[{1/x + 1/y == 1/z, a/x  + b/y == m/n, 2 <= x <= 30, 2 <= y <= 30, 1 <= z <= 30, 2 <= a <= 10, 2 <= b <= 10, a > b, 2 <= m <= 10, 2 <= n <= 10, x > y, GCD[m, n] == 1, GCD[a, b, m] == 1, m <= n}, {x, y, z, a, b, m, n}, Integers]


For a long time, I can not get the result. How can I reduce the time?

• Try splitting the problem. Solve[{1/x + 1/y == 1/z, 2 <= y < x <= 30, 1 <= z <= 30}, {x, y, z}, Integers] returns a reasonable number of candidates, which you can then test with your other criteria. Jan 30, 2021 at 4:27

## 4 Answers

Inspired by @Eric Towers

Clear["*"];

cf = With[{eps = 10^-12.},
Compile[{},
Do[
If[And[GCD[m, n] == 1, Abs[a/x + b/y - m/n] < eps], Sow@{x, y, z, a, b, m, n}],
{x, 2, 30},
{y, 2, x - 1},
{z, 2, 30},
{a, If[Abs[1/x + 1/y - 1/z] > eps, Continue[]]; 2, 10},
{b, 2, a - 1},
{m, 2, 10},
{n, If[GCD[a, b, m] != 1, Continue[]]; m + 1, 10}
]]];

AbsoluteTiming[ans = Reap[cf[]][[2, 1]]]
Length[ans]


Output

{0.0080338, {{12,4,3,3,2,3,4},{12,4,3,4,2,5,6},{12,6,4,4,3,5,6},{12,6,4,5,2,3,4},{12,6,4,6,2,5,6},{15,10,6,3,2,2,5},{15,10,6,6,2,3,5},{15,10,6,6,3,7,10},{15,10,6,6,5,9,10},{15,10,6,7,2,2,3},{15,10,6,8,3,5,6},{15,10,6,9,2,4,5},{18,9,6,4,3,5,9},{18,9,6,6,2,5,9},{18,9,6,6,3,2,3},{18,9,6,6,4,7,9},{18,9,6,6,5,8,9},{18,9,6,7,4,5,6},{18,9,6,8,3,7,9},{18,9,6,9,3,5,6},{18,9,6,10,2,7,9},{18,9,6,10,3,8,9},{20,5,4,4,2,3,5},{20,5,4,4,3,4,5},{20,5,4,6,2,7,10},{20,5,4,7,2,3,4},{20,5,4,10,2,9,10},{24,8,6,3,2,3,8},{24,8,6,6,3,5,8},{24,8,6,6,4,3,4},{24,8,6,6,5,7,8},{24,8,6,7,3,2,3},{24,8,6,8,4,5,6},{24,8,6,9,2,5,8},{24,8,6,9,4,7,8},{24,12,8,5,2,3,8},{24,12,8,6,5,2,3},{24,12,8,7,4,5,8},{24,12,8,8,5,3,4},{24,12,8,8,6,5,6},{24,12,8,9,3,5,8},{24,12,8,9,6,7,8},{24,12,8,10,3,2,3},{24,12,8,10,4,3,4},{28,21,12,4,3,2,7},{28,21,12,8,3,3,7},{30,6,5,5,3,2,3},{30,6,5,5,4,5,6},{30,6,5,6,3,7,10},{30,6,5,7,4,9,10},{30,6,5,8,2,3,5},{30,6,5,9,3,4,5},{30,6,5,10,3,5,6},{30,15,10,5,2,3,10},{30,15,10,6,3,2,5},{30,15,10,8,5,3,5},{30,15,10,9,6,7,10},{30,15,10,9,8,5,6},{30,15,10,10,4,3,5},{30,15,10,10,5,2,3},{30,15,10,10,7,4,5},{30,20,12,6,2,3,10},{30,20,12,9,2,2,5},{30,20,12,9,8,7,10}}}

64

Previous answer

Clear["*"];

gcd = Function[{a0, b0},
Module[{a = a0, b = b0, t}, While[b != 0, {a, b} = {b, Mod[a, b]}];a]];

cf = With[{gcd = gcd, eps = 10^-12.},
Compile[{{x, _Integer}},
Module[{bag = InternalBag[Most@{0}]},
Do[If[And[Abs[1/x+1/y-1/z]<eps,Abs[a/x+b/y-m/n]<eps,gcd[m,n]==1,gcd[gcd[a,b],m]==1],
InternalStuffBag[bag, {x, y, z, a, b, m, n}, 1]],
{y, 2, x - 1}, {z, 2, 30}, {a, 2, 10}, {b, 2, a - 1}, {m, 2, 10}, {n, m + 1, 10}];
InternalBagPart[bag, All]],
RuntimeAttributes -> {Listable}, CompilationTarget -> "C"]];

AbsoluteTiming[ans = Join @@ (Partition[#, 7] & /@ cf[Range[2, 30]])]
Length[ans]


Verifying

And @@ Function[{x,y,z,a,b,m,n}, And[a>b,x>y,n>m,1/x+1/y==1/z,a/x+b/y==m/n,
GCD[m,n]==1,GCD[a,b,m]==1]] @@@ ans


True

• I want to know how you master the compiling skills, which is very complicated for me. Jan 30, 2021 at 6:31
• @chyanog ,very sophisticated solution. You got 42 solutions. With the method @J.M. proposes in his comment, you get 64 solutions verified. (sol3 = Solve[{1/xx + 1/yy == 1/zz, 2 <= yy < xx <= 30, 1 <= zz <= 30}, {xx, yy, zz}, Integers]) // AbsoluteTiming  . (sol5 = DeleteCases[Solve[{1/x + 1/y == 1/z, a/x + b/y == m/n, 2 <= a <= 10, 2 <= b <= 10, a > b, 2 <= m <= 10, 2 <= n <= 10, x > y, m <= n, GCD[m, n] == 1, GCD[a, b, m] == 1, x == xx, y == yy, z == zz} /. #, {x, y, z, a, b, m, n}, Integers] & /@ sol3, {}]) // AbsoluteTiming  Jan 30, 2021 at 12:42
• @Akku14 Thanks for pointing out, I ignored the floating point precision issue at first. Jan 30, 2021 at 13:35
Clear["Global*"]

system = And @@ {1/x + 1/y == 1/z, a/x + b/y == m/n, 2 <= x <= 30,
2 <= y <= 30, 1 <= z <= 30, 2 <= a <= 10, 2 <= b <= 10, a > b,
2 <= m <= 10, 2 <= n <= 10, x > y, GCD[m, n] == 1, GCD[a, b, m] == 1,
m <= n};


You can use FindInstance to identify a few solutions

(sol = FindInstance[system, {x, y, z, a, b, m, n}, Integers,
2]) // AbsoluteTiming

(* {307.661, {{x -> 30, y -> 15, z -> 10, a -> 9, b -> 6, m -> 7,
n -> 10}, {x -> 15, y -> 10, z -> 6, a -> 6, b -> 2, m -> 3, n -> 5}}} *)


Verifying,

system /. sol

(* {True, True} *)


You don't actually have that many cases:

• $$x$$ and $$y$$: $$\frac{1}{2} 30 \cdot 30 = 450$$.
• $$z$$: $$30$$.
• $$a$$ and $$b$$: $$\frac{1}{2} 9 \cdot 9 \approx 40$$.
• $$m$$ and $$n$$: $$\frac{1}{2} 9 \cdot 9 \approx 40$$.

So $$450 \cdot 30 \cdot 40 \cdot 40 = 21.6 \times 10^6$$ cases. And this doesn't discard cases that can be early rejected, for instance, rejecting those $$z$$ not satisfying $$1/x + 1/y = 1/z$$. Generating the cases directly and discarding early, we get something like

compute := Reap[
For[x = 2, x <= 30, x++,
For[y = 2, y < x, y++,
For[z = 1, z <= 30, z++,
If[1/x + 1/y == 1/z,
For[a = 2, a <= 10, a++,
For[b = 2, b < a, b++,
For[m = 2, m <= 10, m++,
If[GCD[a, b, m] == 1,
For[n = m, n <= 10, n++,
If[And[a/x + b/y == m/n, GCD[m, n] == 1],
Sow[{x, y, z, a, b, m, n}];
];
]
];
]
]
]
];
]
]
];
Clear[x, y, z, a, b, m, n];
][[2, 1]]


Then,

RepeatedTiming[
compute
]
(*  {0.130, {
{12, 4, 3, 3, 2, 3, 4}, {12, 4, 3, 4, 2, 5, 6},
{12, 6, 4, 4, 3, 5, 6}, {12, 6, 4, 5, 2, 3, 4},
{12, 6, 4, 6, 2, 5, 6}, {15, 10, 6, 3, 2, 2, 5},
{15, 10, 6, 6, 2, 3, 5}, {15, 10, 6, 6, 3, 7, 10},
{15, 10, 6, 6, 5, 9, 10}, {15, 10, 6, 7, 2, 2, 3},
{15, 10, 6, 8, 3, 5, 6}, {15, 10, 6, 9, 2, 4, 5},
{18, 9, 6, 4, 3, 5, 9}, {18, 9, 6, 6, 2, 5, 9},
{18, 9, 6, 6, 3, 2, 3}, {18, 9, 6, 6, 4, 7, 9},
{18, 9, 6, 6, 5, 8, 9}, {18, 9, 6, 7, 4, 5, 6},
{18, 9, 6, 8, 3, 7, 9}, {18, 9, 6, 9, 3, 5, 6},
{18, 9, 6, 10, 2, 7, 9}, {18, 9, 6, 10, 3, 8, 9},
{20, 5, 4, 4, 2, 3, 5}, {20, 5, 4, 4, 3, 4, 5},
{20, 5, 4, 6, 2, 7, 10}, {20, 5, 4, 7, 2, 3, 4},
{20, 5, 4, 10, 2, 9, 10}, {24, 8, 6, 3, 2, 3, 8},
{24, 8, 6, 6, 3, 5, 8}, {24, 8, 6, 6, 4, 3, 4},
{24, 8, 6, 6, 5, 7, 8}, {24, 8, 6, 7, 3, 2, 3},
{24, 8, 6, 8, 4, 5, 6}, {24, 8, 6, 9, 2, 5, 8},
{24, 8, 6, 9, 4, 7, 8}, {24, 12, 8, 5, 2, 3, 8},
{24, 12, 8, 6, 5, 2, 3}, {24, 12, 8, 7, 4, 5, 8},
{24, 12, 8, 8, 5, 3, 4}, {24, 12, 8, 8, 6, 5, 6},
{24, 12, 8, 9, 3, 5, 8}, {24, 12, 8, 9, 6, 7, 8},
{24, 12, 8, 10, 3, 2, 3}, {24, 12, 8, 10, 4, 3, 4},
{28, 21, 12, 4, 3, 2, 7}, {28, 21, 12, 8, 3, 3, 7},
{30, 6, 5, 5, 3, 2, 3}, {30, 6, 5, 5, 4, 5, 6},
{30, 6, 5, 6, 3, 7, 10}, {30, 6, 5, 7, 4, 9, 10},
{30, 6, 5, 8, 2, 3, 5}, {30, 6, 5, 9, 3, 4, 5},
{30, 6, 5, 10, 3, 5, 6}, {30, 15, 10, 5, 2, 3, 10},
{30, 15, 10, 6, 3, 2, 5}, {30, 15, 10, 8, 5, 3, 5},
{30, 15, 10, 9, 6, 7, 10}, {30, 15, 10, 9, 8, 5, 6},
{30, 15, 10, 10, 4, 3, 5}, {30, 15, 10, 10, 5, 2, 3},
{30, 15, 10, 10, 7, 4, 5}, {30, 20, 12, 6, 2, 3, 10},
{30, 20, 12, 9, 2, 2, 5}, {30, 20, 12, 9, 8, 7, 10}
}}  *)


IF we're curious how many cases are actually fully generated (so actually accounting for the early rejections), we can add a counter and modify the inner For loop.

count = 0;
Reap[
For[x = 2, x <= 30, x++,
...
For[n = m, n <= 10, n++,
count++;
If[And[a/x + b/y == m/n, GCD[m, n] == 1],
Sow[{x, y, z, a, b, m, n}];
];
]
...
];
Clear[x, y, z, a, b, m, n];
][[2, 1]];
count

(* 15816 *)


So several cases are rejected by the Ifs.

Alternatively, we can stage Solves similar to the structure above: first generate all possible x,y,z triples, then try to solve for a,b,m,ns for each of those.

Clear[x, y, z, a, b, m, n];
extend[xyz_] := ({x, y, z, a, b, m, n} /.
Join[xyz, #]) & /@
Solve[{a/x + b/y == m/n, 2 <= a <= 10, 2 <= b <= 10, a > b,
2 <= m <= 10, 2 <= n <= 10, GCD[m, n] == 1, GCD[a, b, m] == 1,
m <= n} /. xyz, {a, b, m, n}, Integers]

Join @@ (extend /@
Solve[{1/x + 1/y == 1/z, 2 <= x <= 30, 2 <= y <= 30, 1 <= z <= 30,
x > y}, {x, y, z}, Integers])


extend is a function taking a triple of the form, for instance, {x -> 6, y -> 3, z -> 2}, one element of the list output by the Solve over {x,y,z}. It uses this to find a, b, m, n that satisfy the rest of the constraints. It's slower than the For loops,

RepeatedTiming[
Join @@ (extend /@
Solve[{1/x + 1/y == 1/z, 2 <= x <= 30, 2 <= y <= 30, 1 <= z <= 30,
x > y}, {x, y, z}, Integers])]

(* {1.153, { ... }}  *)


but faster than "a long time".

• The last solution is more or less what I had in mind in my comment. (+1) Jan 31, 2021 at 4:26

If you want speed for larger upper bounds on the variables $$\{x,y,z,a,b,m,n\}$$, then the following is much faster than even the compiled solution by @chyanog (+1).

Rewrite the original equation $$1/x+1/y-1/z=0$$ as $$(x-z)(y-z)=z^2$$, after forming a common denominator $$x y z>0$$, and multiplying through by it. If $$r*s$$ is a factorization of a given $$z^2$$ on the right-hand side, then the solutions are $$\{x,y\}=\{r,s\}+z$$. This is from making a correspondence between the two factors on the left-hand side of the equation and those of $$r*s$$ on the right.

Generate divisor pairs $$r*s=z^2$$ with the larger one first.

DivisorPairs[n_] := Thread[{Reverse[#], #}][[ ;; Ceiling[Length[#]/2]]] &[Divisors[n]]


Run through the range of $$z$$, generating solution triples $$\{x,y,z\}$$. Note, the maximum $$z$$ is only half the maximum of $$x$$ or $$y$$.

xyzSolve[max_] :=
Block[{xyz},
Flatten[Table[
Pick[#, UnitStep[max - #[[All, 1]]], 1] &[
Map[Join[#, {z}] &, z + Most[DivisorPairs[z^2]]]],
{z, 2, max/2 - 1}], 1]]


For example, there are 12 solutions to $$1/x+1/y=1/z$$, $$x>y$$, $$2\le x\le 30$$, $$2\le y\le 30$$, $$1\le z\le 30$$.

xyzSolve[30]


{{6, 3, 2}, {12, 4, 3}, {20, 5, 4}, {12, 6, 4}, {30, 6, 5}, {24, 8, 6}, {18, 9, 6}, {15, 10, 6}, {24, 12, 8}, {30, 15, 10}, {30, 20, 12}, {28, 21, 12}}

For any solution $$\{x,y\}$$ of $$1/x+1/y=1/z$$, and given $$\{m,n\}$$, the second equation $$a/x+b/y=m/n$$ becomes a simple linear Diophantine equation $$a (n y)+b (n x)=m x y$$ for the unknowns $$\{a,b\}$$. Allow for $$b_{min}\le b < a\le a_{max}$$, and solve via an extended GCD.

abSolve[u_Integer, v_Integer, n_Integer, bmin_, amax_] :=
Block[{g, r, s, k, ab},
{g, {r, s}} = ExtendedGCD[u, v];
ab = Table[{r*n + v*k, s*n - u*k}/g, {k, Ceiling[-r*n/v], Quotient[s*n, u]}];
Select[Pick[ab, Denominator[ab[[All, 1]]], 1], amax >= #[[1]] > #[[2]] >= bmin &]
]


Given $$n$$, use the following RelativePrimes[n] to find allowed $$m\le n$$.

RelativePrimes[n_Integer] := Pick[#, GCD[n, #], 1] &[Range[n]]


Given a solution $$\{x,y,z\}$$, find corresponding values $$\{a,b,m,n\}$$, assuming $$n\ge 2$$ and $$m\ge 2$$.

getSolution[{x_, y_, z_}, bmin_, amax_, nmax_] :=
Pick[#, GCD @@@ (#[[All, {4, 5, 6}]]), 1] &[Flatten[Table[
Map[Join[{x, y, z}, #, {m, n}] &, abSolve[n y, n x, m x y, bmin, amax]],
{n, 2, nmax}, {m, Rest[RelativePrimes[n]]}], 2]]


Solve the original problem in about 0.01 s.

RepeatedTiming[
Length[Flatten[Map[getSolution[#, 2, 10, 10] &, xyzSolve[30]], 1]]
]


{0.0094, 64}

If the upper limit is increased to 100 from 30, then this code takes about 6.5 seconds to find the 19841 valid solutions. The compiled solution takes about 85 times longer.

RepeatedTiming[
Length[Flatten[Map[getSolution[#, 2, 100, 100] &, xyzSolve[100]], 1]]
]


{6.32, 19841}