9
$\begingroup$

I am trying to solve a simple 'sum and product' problem. Basically, the sum of four prices and the product of those four prices are numerically equal (you must, of course, disregard the units). The value of the sum/product is given. Here is my attempt in MMA 11.3.

prices = {a, b, c, d}
pints = prices \[Element] Integers
ppos = prices > 0
ptotal = Total[prices] == 711
pprod = Times @@ prices == 100^3 711
FindInstance[{pints, ppos, ptotal, pprod}, prices]

This spins and spins without returning with an answer.

This Python code utilizing the Z3 Solver module will return within a minute with the answer.

import z3    
def Find711PricesZ3(target = 7.11):
    a, b, c, d = z3.Ints('a b c d')
    s = z3.Solver()
    s.add((a + b + c + d) == target * 100)
    s.add((a * b * c * d) == target * 100**4)
    s.add(a > -1, b > -1, c > -1, d > -1)
    solvable = s.check()
    print(solvable)
    if solvable == z3.sat:
        m = s.model()
        print(m)

Then you run Find711PricesZ3() and voila! Less than a minute later you have an answer.

What is wrong with my Mathematica formulation?

Also, I was trying to Monitor[] the FindInstance[] command with Monitor[...,prices] and various other alternatives, but none of those worked. Is it possible to see which values of prices Mathematica is working on to know that the process is working?

Thanks so much for your assistance.

$\endgroup$
10
  • $\begingroup$ @MichaelE2 Do you not see my Mathematica code in the post? I see it on my end. It is above the Python code. Let me know if you still don't see it and I will try to figure out why. $\endgroup$
    – JoeMarino
    Jun 24 at 22:52
  • $\begingroup$ Sorry, Joe, for some reason my brain parsed it as python. <face palm> $\endgroup$
    – Michael E2
    Jun 24 at 22:54
  • $\begingroup$ FindInstance does some symbolic solving and not just a brute-force search. This shows the (first) Reduce command it calls: Trace[ FindInstance[{pints, ppos, ptotal, pprod}, prices], r_Reduce :> Return[Hold[r], Trace], TraceInternal -> True] $\endgroup$
    – Michael E2
    Jun 24 at 22:58
  • $\begingroup$ @MichaelE2 Does Trace need to finish evaluating the expression before it returns with the first Reduce command it calls? This Trace seems to be spinning for me again, not coming back with an answer, as if it is trying to complete the FindInstance. Maybe Trace only returns once expression has been fully evaluated? $\endgroup$
    – JoeMarino
    Jun 24 at 23:18
  • $\begingroup$ The code I posted uses Return to break out of Trace before the Reduce command is executed. It's very quick for me (less than 0.01 sec. according to AbsoluteTiming). Did you change any of your variable definitions since posting your Q? I assume you copy-pasted my code, but if you typed it in, check it. The Hold[r] is important because that keeps the Reduce command from executing. $\endgroup$
    – Michael E2
    Jun 24 at 23:25

5 Answers 5

13
$\begingroup$

Mathematica can find all the solutions.

Reduce[{a + b + c + d == 711, a*b*c*d == 100^3 711}, {a, b, c, 
  d}, PositiveIntegers]

%// ToRules // List
{{a -> 120, b -> 125, c -> 150, d -> 316}, {a -> 120, b -> 125, 
  c -> 316, d -> 150}, {a -> 120, b -> 150, c -> 125, 
  d -> 316}, {a -> 120, b -> 150, c -> 316, d -> 125}, {a -> 120, 
  b -> 316, c -> 125, d -> 150}, {a -> 120, b -> 316, c -> 150, 
  d -> 125}, {a -> 125, b -> 120, c -> 150, d -> 316}, {a -> 125, 
  b -> 120, c -> 316, d -> 150}, {a -> 125, b -> 150, c -> 120, 
  d -> 316}, {a -> 125, b -> 150, c -> 316, d -> 120}, {a -> 125, 
  b -> 316, c -> 120, d -> 150}, {a -> 125, b -> 316, c -> 150, 
  d -> 120}, {a -> 150, b -> 120, c -> 125, d -> 316}, {a -> 150, 
  b -> 120, c -> 316, d -> 125}, {a -> 150, b -> 125, c -> 120, 
  d -> 316}, {a -> 150, b -> 125, c -> 316, d -> 120}, {a -> 150, 
  b -> 316, c -> 120, d -> 125}, {a -> 150, b -> 316, c -> 125, 
  d -> 120}, {a -> 316, b -> 120, c -> 125, d -> 150}, {a -> 316, 
  b -> 120, c -> 150, d -> 125}, {a -> 316, b -> 125, c -> 120, 
  d -> 150}, {a -> 316, b -> 125, c -> 150, d -> 120}, {a -> 316, 
  b -> 150, c -> 120, d -> 125}, {a -> 316, b -> 150, c -> 125, 
  d -> 120}}
$\endgroup$
7
  • $\begingroup$ +1, I did try Reduce, but used Integers and it was taking long time so stopped it. did not think of trying domain PositiveIntegers Mathematica beats Python again :) $\endgroup$
    – Nasser
    Jun 25 at 0:44
  • 1
    $\begingroup$ Are all the solutions equivalent? $\endgroup$
    – Michael E2
    Jun 25 at 0:49
  • $\begingroup$ @MichaelE2 Yes, all of them are equivalent upto permutation. $\endgroup$
    – cvgmt
    Jun 25 at 0:55
  • $\begingroup$ Thanks,I thought maybe my eyes missed something. How long did it take? $\endgroup$
    – Michael E2
    Jun 25 at 0:56
  • $\begingroup$ @MichaelE2 About 10min $\endgroup$
    – cvgmt
    Jun 25 at 0:57
10
$\begingroup$

Second answer

Here's a fast way to solve the problem, working from a different angle based on the factorization of 100^3 * 711, which obviously figures in pprod:

Extract[#, Position[Total[#, {5}], 711]] &@
   Outer[Times,
    2^Select[
      IntegerPartitions[6, {4}, Range[0, Floor[Log[2, 711]]]],
      Times @@ # == 0 &], (* at least one odd price to get 711 *)
    3^Flatten[
      Permutations /@ 
       IntegerPartitions[2, {4}, Range[0, Floor[Log[3, 711]]]], 1],
    5^Flatten[
      Permutations /@ 
       IntegerPartitions[6, {4}, Range[0, Floor[Log[5, 711]]]], 1],
    79^Flatten[
      Permutations /@ 
       IntegerPartitions[1, {4}, Range[0, Floor[Log[79, 711]]]], 1],
    1
    ] // Flatten // AbsoluteTiming

(*  {0.017399, {120, 316, 150, 125}}  *)

First answer

Oh, cool, while I was doing something else, this finished:

prices = {a, b, c, d};
pints = prices \[Element] Integers;
ppos = And @@ Thread[prices > 0];
ptotal = Total[prices] == 711;
pprod = Times @@ prices == 100^3 711;
FindInstance[{pints, ppos, ptotal, pprod, 
  0 < a < 500 && 0 < b < 500 && 0 < c < 500 && 
   0 < d < 500}, prices, Integers]

(*  {{a -> 120, b -> 125, c -> 150, d -> 316}}  *)

All I know is that it took a long time. I neglected to time it. Less than 15 minutes, I think.

$\endgroup$
2
  • $\begingroup$ Wow! This second answer is magnificent! I am currently digesting it. Not often can one realize a 6000x increase in performance so soon. Is the reason the factor 2 doesn't get permuted because one of the factors can sit still, while all the others permute around it? That will still generate all possibilities? $\endgroup$
    – JoeMarino
    Jun 25 at 15:37
  • $\begingroup$ @JoeMarino Re permuting 2: Take any solution {a, b, c, d}. You can reorder it so that the 2-power factors appear in one of the orders given by the Select[..] code. I chose 2 because its partitions gave the most permutations. One could choose any of the prime factors. It just changes how many redundant solutions you get (redundant up to reordering). $\endgroup$
    – Michael E2
    Jun 25 at 15:54
6
$\begingroup$

OK. So, the original code I posted was not faulty. I didn't have a code problem, I had an old computer problem. As I watched in Task Manager, when Mathematica would run the original code, it would gobble up RAM, and at some point reach the limit of physical RAM on the machine. At that point, Mathematica would begin transferring some of what was in RAM to disk (noted by a high disk % and a reduction of RAM being used). This would repeat until terminated (I couldn't wait for this to finish; when you start using hard disk instead of RAM, you are in for a long wait).

As @MichaelE2 and @cvgmt have noted, on a modern machine, the first three iterations of the code posted will run and complete in about 10 or 15 minutes. FindInstance[] is capable of finding more than one solution, and this robustness makes it slower than the Python code, which is definitely only designed to find just one solution.

We have improvements to this, of course. MichaelE2 created a version that utilizes the factorization of the product constraint. It runs quickly. Since it can help a new Mathematica user (and I am kinda new to Mathematica), I'd like to break down his code for the less experienced.

The prime factors of 711 million are 2, 3, 5, and 79.

FactorInteger[711 100^3]

{{2, 6}, {3, 2}, {5, 6}, {79, 1}}

Now each power above (6, 2, 6 and 1) is partitioned, limited by the number of times that the base will divide 711. Only partitions of length 4 are needed, since we have four prices.

IntegerPartitions[6, {4}, Range[0, Floor[Log[2, 711]]]]
IntegerPartitions[2, {4}, Range[0, Floor[Log[3, 711]]]]
IntegerPartitions[6, {4}, Range[0, Floor[Log[5, 711]]]]
IntegerPartitions[1, {4}, Range[0, Floor[Log[79, 711]]]]

{{6, 0, 0, 0}, {5, 1, 0, 0}, {4, 2, 0, 0}, {4, 1, 1, 0}, {3, 3, 0,
0}, {3, 2, 1, 0} ....

Now, we need to create permutations for three of these four (see MichaelE2's note above for why). We'll take a look at the factor 5.

Permutations /@ 
  IntegerPartitions[6, {4}, Range[0, Floor[Log[5, 711]]]] // 
 Flatten[#, 1] &

{{4, 2, 0, 0}, {4, 0, 2, 0}, {4, 0, 0, 2}, {2, 4, 0, 0}, {2, 0, 4,
0}, {2, 0, 0, 4}, {0, 4, 2, 0}, {0, 4, 0, 2}, {0, 2, 4, 0},...

Now, we raise the base factor to each of these exponents. Then Outer[] does its job of making every possible combination (grab one from each list and apply the given function, here Times[]).

Outer[Times, 
 2^Select[IntegerPartitions[6, {4}, Range[0, Floor[Log[2, 711]]]], 
   Times @@ # == 0 &],(*at least one odd price to get 711*)
 3^Flatten[
   Permutations /@ 
    IntegerPartitions[2, {4}, Range[0, Floor[Log[3, 711]]]], 1], 
 5^Flatten[
   Permutations /@ 
    IntegerPartitions[6, {4}, Range[0, Floor[Log[5, 711]]]], 1], 
 79^Flatten[
   Permutations /@ 
    IntegerPartitions[1, {4}, Range[0, Floor[Log[79, 711]]]], 1], 1]

{{{{{28440000, 25, 1, 1}, {360000, 1975, 1, 1}, {360000, 25, 79, 1}, {360000, 25, 1, 79}},...

Outer[] has first grabbed 2^{6, 0, 0, 0}, 3^{2, 0, 0, 0}, 5^{4, 2, 0, 0}, and 79^{1, 0, 0, 0} and created the products of 2^i * 3^i * 5^i * 79^i for i = 1 to 4. It will then change to 79^{0, 1, 0, 0}, and repeat. It does this for all combos.

Now we only have to execute Extract[#, Position[Total[#, {5}], 711]] & on the result of all this. Each item is totaled, and if that item totals 711, that location is fed to Extract[] to retrieve the four prices from the original, un-totaled list that are the solution. Basically, for each set of four prices created by Outer[], if the total is 711, retrieve the prices.

I know when I first started looking at compact Mathematica code, it could be daunting. I hope this extended look at MichaelE2's solution helps the newer users out there who stumble across this later.

I am still trying to digest chyanog's solution. I still need to decipher the first command executed, which is a lot of the secret sauce. When I figure it out, I'll try to post back. I can vouch that the code does work. If you have repeated applications to make, the compiled function will be speedy.

EDIT

OK. I think I get it now. Here is the first command from chyanog's solution.

Through[{Floor@*MinValue, Ceiling@*MaxValue}[{#, a + b + c + d == 711,
      a*b*c*d == 100^3 711, 0 < a <= b <= c <= d}, {a, b, c, d}]] & /@ {a, b, c, d}

So, Through[] is a new one for me. So is the Composition[] function. The composition, using the shorthand @*, is nothing more than Floor[MinValue[...]] and Ceiling[MaxValue[...]. What Through[] is going to do is take these two compositions and apply them both to whatever comes next. That means we end up taking the arguments

{#, a + b + c + d == 711,a*b*c*d == 100^3 711, 0 < a <= b <= c <= d}, {a, b, c, d}]] &

and feed them to both Floor[MinValue[...]] and Ceiling[MaxValue[...] (this ends up being the expression and constraints that those functions want). What MinValue[] and MaxValue[] will do is look at the constraints it was given and provide a range of values that a variable can take and still satisfy the constraints.

The constraints were formulated as a function (note the # placeholder and the & terminator), so we can Map[] it (using /@) onto each variable separately. We then create the pairs in the output by executing this pseudocode:

For each var in {a, b, c, d}, find the minimum and maximum values that the variable takes on while satisfying the constraints.

The first pair in the output, {74, 131}, is created from these two commands:

Floor[MinValue[{a, a + b + c + d == 711,a*b*c*d == 100^3 711, 0 < a <= b <= c <= d}, {a, b, c, d}]]

and

Ceiling[MaxValue[{a, a + b + c + d == 711,a*b*c*d == 100^3 711, 0 < a <= b <= c <= d}, {a, b, c, d}]]

From there, I feel the code speaks for itself.

END EDIT

Thanks, everyone, for the valuable contributions. If I have interpreted any of the code incorrectly, hopefully someone will let us know.

$\endgroup$
5
$\begingroup$

Without loss of generality, assuming $a\leq b\leq c\leq d$ , we can determine the range of $a,b,c,d$

Through[{Floor@*MinValue, Ceiling@*MaxValue}[{#, a + b + c + d == 711,
      a*b*c*d == 100^3 711, 0 < a <= b <= c <= d}, {a, b, c, d}]] & /@ {a, b, c, d}

{{74, 131}, {107, 213}, {130, 248}, {212, 320}}

$c, d$ can be represented by $a, b$

Solve[{a + b + c + d == 711, a*b*c*d == 100^3 711}, {c, d}][[1]] /. 
  Sqrt[_] :> Sqrt[delta] // Simplify

$\left\{c\to \frac{1}{2} \left(-\frac{\sqrt{\text{delta}}}{a b}-a-b+711\right),d\to \frac{1}{2} \left(\frac{\sqrt{\text{delta}}}{a b}-a-b+711\right)\right\}$

So we can get the answer by enumerating $a$ and $b$

Compile[{}, Module[{delta, c, d, tol = 1.*^-10},
    Do[delta = a b (a b (a + b - 711)^2 - 4 100^3 711);
     If[delta >= 0,
      c = 1/2 (711 - a - b - Sqrt[delta]/(a b)); 
      d = 1/2 (711 - a - b + Sqrt[delta]/(a b));
      If[a <= b && Abs[Round@c - c] < tol && Abs[Round@d - d] < tol &&
         130 < c < 248 && 212 < d < 320, Print[Round@{a, b, c, d}]];
      ]
     , {a, 74., 131}, {b, 107., 213}]
    ]][] // AbsoluteTiming

{120, 125, 150, 316}
{0.0021569, Null}

If you need all possibilities
Permutations[{120, 125, 150, 316}]

{{120,125,150,316},{120,125,316,150},{120,150,125,316},{120,150,316,125},{120,316,125,150},{120,316,150,125},{125,120,150,316},{125,120,316,150},{125,150,120,316},{125,150,316,120},{125,316,120,150},{125,316,150,120},{150,120,125,316},{150,120,316,125},{150,125,120,316},{150,125,316,120},{150,316,120,125},{150,316,125,120},{316,120,125,150},{316,120,150,125},{316,125,120,150},{316,125,150,120},{316,150,120,125},{316,150,125,120}}

$\endgroup$
4
  • $\begingroup$ Thank you so much for your contribution. I will definitely study this solution. I have a long way to go to master when to compile and when not to compile. Also, the underlying mathematics here will take me a second to pick up. Who knew such a 'simple' puzzle could get such advanced treatment? $\endgroup$
    – JoeMarino
    Jun 25 at 15:40
  • $\begingroup$ Wait, the range for 'd' was determined to be [212, 320], yet the Do loop uses [21, 320]. I presume that is a typo? $\endgroup$
    – JoeMarino
    Jun 25 at 17:04
  • $\begingroup$ @JoeMarino This is a typo. $\endgroup$
    – chyanog
    Jun 26 at 8:29
  • $\begingroup$ If all three steps are timed together, the total time is 0.217672 sec., and the first step is rather an essential step in the solution. $\endgroup$
    – Michael E2
    Jun 26 at 15:57
3
$\begingroup$

NMinimize can do the job in 4 seconds.

( For some parameters, you have to play around with WorkingPrecision, MaxIterations and PrecisionGoal. )

prices = {a, b, c, d};
pints = # \[Element] Integers & /@ prices;
ppos = 0 < # < 711 & /@ prices;
ptotal = Total[prices] == 711;
pprod = Times @@ prices == 100^3 711;

jj = Join[pints, ppos, {ptotal}, {pprod}]

(*   {a \[Element] Integers, b \[Element] Integers, c \[Element] Integers, 
 d \[Element] Integers, 0 < a < 711, 0 < b < 711, 0 < c < 711, 
 0 < d < 711, a + b + c + d == 711, a b c d == 711000000}   *)

(nmin = NMinimize[{a, jj}, prices, MaxIterations -> 500, 
    WorkingPrecision -> 30]) // Timing

(*   {4.234, {125.000000000000000000000000000, {a -> 125, b -> 316, 
   c -> 120, d -> 150}}}   *)

(nmin = NMinimize[{d, jj}, prices, MaxIterations -> 500, 
    WorkingPrecision -> 200, AccuracyGoal -> 1, 
    PrecisionGoal -> 1]) // Timing

(*   {4.14, {120.\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000, {a -> 125, 
   b -> 316, c -> 150, d -> 120}}}   *)
$\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.