The excellent second solution by Roman, with R
slightly modified, produces
R = HornerForm[(a (3 + a) (-12 + 5 a) + 3 (9 - 14 b) b)/
((-1 + a) a (3 + a) - 3 (-1 + b) b)]
With[{s = 10^4}, Do[If[Divisible[a (3 + a) (-12 + 5 a) + 3 (9 - 14 b) b,
(-1 + a) a (3 + a) - 3 (-1 + b) b] && R >= 3, Sow[{a, b, R}]],
{a, s}, {b, s}] // Reap // Last // First]
(* {{3, 6, 24}, {5, 8, 244}, {5, 10, 31}, {5, 14, 19}, {9, 18, 177}, {9, 20, 46},
{12, 30, 45}, {32, 112, 139}, {33, 114, 573}, {35, 126, 220}, {45, 180, 553},
{47, 450, 16}, {48, 204, 129}, {63, 294, 3750}, {77, 396, 3889}, {116, 728, 46750},
{117, 2340, 15}, {159, 1166, 6826}, {240, 2156, 2098129}, {243, 2214, 576},
{357, 3906, 72807}, {372, 4154, 2509849}, {492, 6314, 398389}} *)
in about 350 seconds. I attempted to find faster approaches using various combinations of Tuples
, Table
, Cases
, and Select
, but the best I could do was
Flatten[Table[If[Divisible[a (-36 + a (3 + 5 a)) + (27 - 42 b) b,
a (-3 + a (2 + a)) + (3 - 3 b) b] && R > 2, {a, b, R}, Nothing, Nothing],
{a, 10000}, {b, 10000}], 1]
which produced the same results in the same amount of time.
The tutorial, DiophantineReduce discusses, among many other cases, "Equations with a Linear Variable", which this question is. Applying Reduce
Reduce[R == r && a > 0 && b > 0 && r > 2, {a, b, r}, Integers]
yields a lengthy result in less than a second, a portion of which is, in effect,
(* b > 1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3]) && r == R *)
(Not coincidentally, 1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3])
is the value of b
for which Denominator[R] == 0
.) Employing the inequality in my approach above,
Flatten[Table[If[Divisible[a (-36 + a (3 + 5 a)) + (27 - 42 b) b,
a (-3 + a (2 + a)) + (3 - 3 b) b] && R > 2, {a, b, R}, Nothing, Nothing], {a, 10000},
{b, Ceiling[1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3])], 10000}], 1]
reproduces the results given at the beginning of this answer in 15 seconds, a significant improvement. Applying this approach to a much larger domain (and using ParallelTable
on a six-processor PC) then yields
Flatten[ParallelTable[If[Divisible[a (-36 + a (3 + 5 a)) + (27 - 42 b) b,
a (-3 + a (2 + a)) + (3 - 3 b) b] && R > 2, {a, b, R}, Nothing, Nothing],
{a, 6000}, {b, Ceiling[1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3])],
300000}], 1]
(* {{3, 6, 24}, {5, 8, 244}, {5, 10, 31}, {5, 14, 19}, {9, 18, 177}, {9, 20, 46},
{12, 30, 45}, {32, 112, 139}, {33, 114, 573}, {35, 126, 220}, {45, 180, 553},
{47, 450, 16}, {48, 204, 129}, {63, 294, 3750}, {77, 396, 3889}, {116, 728, 46750},
{117, 2340, 15}, {159, 1166, 6826}, {240, 2156, 2098129}, {243, 2214, 576},
{357, 3906, 72807}, {372, 4154, 2509849}, {492, 6314, 398389}, {768, 12336, 1769},
{1266, 26028, 12553000}, {1545, 43860, 30}, {3792, 138336, 186},
{5973, 266574, 121035}} *)
in 1070 seconds. Here is a plot of b
vs a
.
Show[ListLogLogPlot[%[[2, All, ;;2]], PlotRange -> All, ImageSize -> Large, AxesLabel ->
{a, b}, LabelStyle -> {14, Bold, Black}], LogLogPlot[1/2 + Sqrt[3 - 12 a + 8 a^2 +
4 a^3]/(2 Sqrt[3]), {a, 1, 10000}, PlotRange -> All]]

Evidently, most of the points lie just above the inequality curve. This suggests that most, although not all, solutions can be obtained by searching just above the curve. For instance,
Flatten[ParallelTable[If[Divisible[a (-36 + a (3 + 5 a)) + (27 - 42 b) b,
a (-3 + a (2 + a)) + (3 - 3 b) b] && R > 2, {a, b, R}, Nothing, Nothing],
{a, 1000000}, {b, Ceiling[1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3])],
Ceiling[1/2 + Sqrt[3 - 12 a + 8 a^2 + 4 a^3]/(2 Sqrt[3])] + 100}], 1]
(* {{3, 6, 24}, {5, 8, 244}, {5, 10, 31}, {5, 14, 19}, {9, 18, 177}, {9, 20, 46},
{12, 30, 45}, {32, 112, 139}, {33, 114, 573}, {35, 126, 220}, {45, 180, 553},
{48, 204, 129}, {63, 294, 3750}, {77, 396, 3889}, {116, 728, 46750},
{159, 1166, 6826}, {240, 2156, 2098129}, {243, 2214, 576}, {357, 3906, 72807},
{372, 4154, 2509849}, {492, 6314, 398389}, {768, 12336, 1769},
{1266, 26028, 12553000}, {5973, 266574, 121035}, {12440, 801136, 1730566},
{43329, 5207358, 30979126197}, {44517, 5422980, 3270113811},
{137796, 29532312, 8075577424022}} *)
in 220 seconds. Plotted as before,

Addendum: Direct Solution with Reduce
Further review of Ref. 1 indicates that Reduce
can obtain integer zeros for bounded regions of {a, b}
, for instance,
SetSystemOptions["ReduceOptions" -> {"DiscreteSolutionBound" -> Infinity}];
SetSystemOptions["ReduceOptions" -> {"SieveMaxPoints" -> {10^3, 10^6}}];
Values@Solve[{r == R, 1000 >= a > 0, 1000 >= b > 0, r > 2}, {a, b, r},
Integers, Method -> Reduce]
yields the same sixteen results obtain by Roman in his answer, but over three times more slowly.
FindInstance
? For example:eqn = a (a + 3) (a (r - 5) + (12 - r))/9 == b (9 + b (-14 + r) - r)/3; constraints = And @@ {a >= 1, b >= 1, r >= 3}; FindInstance[eqn && constraints, {a, b, r}, PositiveIntegers]
gives{{a -> 5, b -> 10, r -> 31}}
$\endgroup$ – flinty Oct 18 '20 at 14:36