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

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 = Internal`Bag[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], 
       Internal`StuffBag[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}];
     Internal`BagPart[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


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".