# Solving an equation in natural numbers

I am trying to solve the following equation in the Natural Numbers, with the condition $$a\ge1$$, $$b\ge1$$, and $$r\ge3$$:

$$\frac{a(a + 3)(a(r - 5) + (12 - r))}{9}=\frac{b (9 + b (-14 + r) - r)}{3}\tag1$$

The method I know use is, that I solve the equation for $$b$$ and I got:

$$b=\displaystyle\frac{1}{6} \left(\sqrt{3\cdot\frac{4 a (a+3) (r-14) (a (r-5)-r+12)+3 (r-9)^2}{(r-14)^2}}+\frac{15}{r-14}+3\right)\tag2$$

Now, I used Mathematica to check when the function under the square root is a perfect square, with the following code:

ParallelTable[
If[IntegerQ[
FullSimplify[
Sqrt[3*((
4 a (3 + a) (12 + a (-5 + r) - r) (-14 + r) +
3 (-9 + r)^2)/(-14 + r)^2)]]], {a, r}, Nothing], {a, 1,
10^5}, {r, 3, 10^5}] //. {} -> Nothing


And the solutions I got, I put in equation $$(1)$$ to check if I can find a solution to the original problem.

This method takes a very very long time, but I am not knowing if there is a faster and smarter way to program this. Can you help me with this. Thanks a lot in advance.

• What's wrong with just using 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}} Oct 18 '20 at 14:36

Much faster to solve for $$r$$ and check that it's an integer: (it's always rational; no square-roots involved)

Solve[(a(a+3)(a(r-5)+(12-r)))/9 == (b(9+b(-14+r)-r))/3, r] // FullSimplify

(*    {{r -> (a(3+a)(-12+5a)+3(9-14b)b)/((-1+a)a(3+a)-3(-1+b)b)}}    *)

R = (a(3+a)(-12+5a)+3(9-14b)b)/((-1+a)a(3+a)-3(-1+b)b);
With[{s = 10^3},
Do[If[IntegerQ[R] && 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}}                *)


Even faster: instead of calculating $$r$$ and checking if it's an integer, we can just check if its numerator is divisible by its denominator:

R = (a(3+a)(-12+5a)+3(9-14b)b)/((-1+a)a(3+a)-3(-1+b)b);
With[{s = 10^3},
Do[If[Divisible[a(3+a)(-12+5a)+3(9-14b)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}}                *)

• I got a lot of error messages when running your last code, about 1/0. Oct 18 '20 at 16:08
• @Jan those are warnings, not errors. You can suppress them if you want with Quiet@Do[If[Divisible[.... Oct 18 '20 at 16:17
• +1, Good idea @Roman You can make it even faster, if you regard a ContourPlot3D[ Evaluate@eqn, {a, 1, 2000}, {b, 1, 50000}, {r, 1000, 1000000}, PlotPoints -> 100]  . For r >1000 the b values strongly depend on a values for all greater r. You can calculate that relation e.g. for r == 1000 bb1000[a_] = b /. First@Solve[{1 < b, 1 < a, eqn /. r -> 1000}, b] // Simplify[#, a > 1] &  and need only test for b values +- 10 or so from that bb1000. Oct 18 '20 at 18:24
• @Akku14 yes indeed, most large solutions seem to have $b\approx a^{3/2}/\sqrt{3}$ as visible in your contour plot, which can be found by taking the limit $r\to\infty$ of the problem (setting the denominator in the expression for R to zero). There are exceptions, though, for example $(a,b,r)=(117,2340,15)$ and generally the solutions with small $r$ (which you've excluded in your plot). I would not be surprised if such exceptional solutions appear even for large $(a,b)$ and would therefore advise against your heuristic. Oct 18 '20 at 19:27
• @Roman, that's clear, that there are exceptional solutions for large (a,b), e.g. {1545, 43860, 30}  . My idea was to look for solutions with r >100 as i suggested above. In a second step You get the high a,b solutions, if you solve eqn for b, get br[a,r] and do your procedure for all a and r up to 100 (tab4 = With[{s = 2 10^3}, Do[If[IntegerQ[bbb = br[a, r]], Sow[{a, bbb, r}]], {a, s}, {r, 14, 100}] // Quiet // Reap // Last // First]) // AbsoluteTiming  . Oct 19 '20 at 4:14

Borrowing a fast perfect-square test from Fastest square number test, and shortening the length of the test case:

(* OP's *)
Table[
If[IntegerQ[
FullSimplify[
Sqrt[3*((4 a (3 + a) (12 + a (-5 + r) - r) (-14 + r) +
3 (-9 + r)^2)/(-14 + r)^2)]]], {a, r}, Nothing], {a, 1,
300}, {r, 3, 30}] // Flatten[#, 1] & // AbsoluteTiming

(*  {83.9498, {{5, 19}, {117, 15}, {252, 29}}}  *)

sQ[n_] := FractionalPart@Sqrt[n + 01] == 0;
Table[
If[IntegerQ[#] && sQ[#] &[
3*((4 a (3 + a) (12 + a (-5 + r) - r) (-14 + r) +
3 (-9 + r)^2)/(-14 + r)^2)], {a, r}, Nothing], {a, 1,
300}, {r, 3, 30}] // Flatten[#, 1] & // AbsoluteTiming

(*  {0.068718, {{5, 19}, {117, 15}, {252, 29}}}  *)


For the $$10^5 \times 10^5$$ search, the improved code will take on the order of 80000 seconds, but that's a lot less than the $$10^8$$ seconds that the OP's would take. (Divide by an appropriate factor if parallelized.)

Try NSolve with restricted parameter range 1<= a,b,r <=50

NSolve[{1/9 a (a + 3) (a (r - 5) + 12 - r) ==1/3 b (9 + b (-14 + r) - r) , 50 >= a >= 1, 50 >= b >= 1 ,50 > r >= 1}, {a, b, r}, Integers]
(**{{a -> 3, b -> 6, r -> 24},
{a -> 5, b -> 10, r -> 31},
{a -> 5,b -> 14, r -> 19},
{a -> 9, b -> 20, r -> 46},
{a -> 12, b -> 30,r -> 45}}*)


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.