1
$\begingroup$

I was working with Project Euler Problem 418 using Mathematica and got in trouble.

I wrote a function to find the unique factorization triple which minimizes c / a for integer n.

FactorizationTriple[n_] := 
 Block[{factor, base, exponent, l, e, varmat, vars, cons1, cons2, cons3, solve, a, b, c},
  factor = FactorInteger[n];
  base = factor[[All, 1]];
  exponent = factor[[All, 2]];
  l = Length[factor];
  varmat = Table[e[i, j], {i, 1, 3}, {j, 1, l}];
  vars = Flatten[varmat];
  cons1 = Thread[0 <= vars];
  cons2 = Thread[exponent == Total[varmat, {1}]];
  cons3 = {varmat[[1]].N[Log[base]] <= varmat[[2]].N[Log[base]] <= varmat[[3]].N[Log[base]]};
  solve = FindMinimum[{varmat[[3]].N[Log[base]] - varmat[[1]].N[Log[base]], Join[cons1, cons2, cons3], Element[vars, Integers]}, vars];
  a = Times @@ (Power @@@ ({base, Table[e[1, j], {j, 1, l}] /. solve[[2]]}\[Transpose]));
  b = Times @@ (Power @@@ ({base, Table[e[2, j], {j, 1, l}] /. solve[[2]]}\[Transpose]));
  c = Times @@ (Power @@@ ({base, Table[e[3, j], {j, 1, l}] /. solve[[2]]}\[Transpose]));
  {a, b, c}]

Since both f and cons in FindMinum are linear, I thought it uses Method->"LinearProgramming" and I expected it to return a global minimum.

FactorizationTriple does work when n = 165 or 100100 or 20!, but it can't give me the correct answer when n = 43!:

AbsoluteTiming[FactorizationTriple[43!]]
{1044.17, {392385912744443904, 392388272221065120, 392389380337500000}}

The correct answer is {a, b, c} = {392386762388275200, 392388272221065120, 392388530688000000}.

Questions:

  1. Should I use NMinimize or LinearProgramming instead? (I had some try but failed.)
  2. How to set the options in FindMinimum?
  3. How to improve the efficiency of FactorizationTriple? (It's too slow now.)
$\endgroup$
  • $\begingroup$ For starters, N[Log[base]] instances in cons3 convert exact values to inexact ones, causing loss of precision in the process. Have you tried without the enclosing N? (Secondly FindMinimum is not guaranteed to find a global minimum, Minimize does so.) $\endgroup$ – kirma Jan 16 at 15:25
  • 2
    $\begingroup$ @Kirma For ILP FindMinimum should find the global min, but as you suggest, this is subject to numeric issues (the library in question uses machine arithmetic to solve relaxed subproblems). Based on the sizes of the result integers, I rather suspect this is the underlying issue. $\endgroup$ – Daniel Lichtblau Jan 16 at 15:39
  • 1
    $\begingroup$ @DanielLichtblau I stand corrected. I was actually suspecting this would be the case, but in general it's not a safe assumption... :) $\endgroup$ – kirma Jan 16 at 15:40
  • $\begingroup$ @kirma I had tried without N before, but even if it calculate it for 1 hour, I can't get the result. $\endgroup$ – 李子涵 Jan 17 at 2:40

Your Answer

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

Browse other questions tagged or ask your own question.