3
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • $\begingroup$ 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}} $\endgroup$ – flinty Oct 18 '20 at 14:36
6
$\begingroup$

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}}                *)
Improve this answer
$\endgroup$
5
  • $\begingroup$ I got a lot of error messages when running your last code, about 1/0. $\endgroup$ – Jan Oct 18 '20 at 16:08
  • $\begingroup$ @Jan those are warnings, not errors. You can suppress them if you want with Quiet@Do[If[Divisible[.... $\endgroup$ – Roman Oct 18 '20 at 16:17
  • $\begingroup$ +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. $\endgroup$ – Akku14 Oct 18 '20 at 18:24
  • $\begingroup$ @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. $\endgroup$ – Roman Oct 18 '20 at 19:27
  • $\begingroup$ @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 . $\endgroup$ – Akku14 Oct 19 '20 at 4:14
3
$\begingroup$

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 + 0``1] == 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.)

Improve this answer
$\endgroup$
3
$\begingroup$

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}}*)
Improve this answer
$\endgroup$
1
$\begingroup$

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]]

enter image description here

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,

enter image description here

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.