7
$\begingroup$

I know Mathematica is really powerful when it comes to functional operations like applying a list of tasks to a list of variables. Sometimes I feel like it isn't the most powerful tool when it comes to looking for a number in a large range. But that's what I'd like to use it for regardless. I can always use a While loop, but I'm not sure that is the most effective way of doing it. So I'm looking for a more elegant/effective way of doing this. Let me give you an example:

Question: Find the smallest positive integer whose sum of its digits squared is greater than 100.

My rather sloppy solution:

myInt = 0;
mySum = 0;
SumLimit = 100;
While[mySum < SumLimit,
 myDigits = IntegerDigits[myInt];
 mySum = Total[#^2 & /@ myDigits];
 myInt = myInt + 1;
]
myInt - 1

And the answer is 59.

But as you can clearly see, this method is quite slow. If I said 1000 instead of 100, it would take quite a long time to find. There are obviously better ways to solve this problem than brute force, but I'd basically like to have a more elegant implementation of the brute force method.

$\endgroup$
5
  • $\begingroup$ No, you can't use Range[100]. What would you use if I said 10^5 instead of 10^2? $\endgroup$
    – cartonn
    Jun 18, 2014 at 20:18
  • $\begingroup$ just tightening it up helps a little: ` While[Total@(IntegerDigits[++myInt]^2) < SumLimit] $\endgroup$
    – george2079
    Jun 18, 2014 at 20:35
  • $\begingroup$ @george2079 I like to have individual statements for the readability of the code. $\endgroup$
    – cartonn
    Jun 18, 2014 at 20:37
  • 1
    $\begingroup$ @Öskå I think it's safe to say most people are bad with Compile and that just shows that it's implementation and documentation are terrible @Kuba. $\endgroup$
    – RunnyKine
    Jun 18, 2014 at 21:03
  • $\begingroup$ Look at a related problem: Finding the largest integer that cannot be partitioned in a certain way $\endgroup$
    – Artes
    Jun 18, 2014 at 23:07

9 Answers 9

9
$\begingroup$

I think this works correctly:

ClearAll[min, doMin];
min[x_] := 
  doMin[x] // Reap // Last // Flatten // Reverse // FromDigits;
doMin[x_] :=
  With[
   {d = Range[9]^2},
   If[
    x > 81,
    Sow@ConstantArray[9, IntegerPart[x/81]]; doMin[Mod[x, 81]],
    Sow@Sqrt@Select[d, # >= x &, 1]]];

min[100] // AbsoluteTiming

{0.001005, 59}

min[1000000] // AbsoluteTiming

{0.024029, .... }

$\endgroup$
8
  • 1
    $\begingroup$ @mfvonh - I agree with the "59", but I don't agree with the "..." $\endgroup$
    – eldo
    Jun 18, 2014 at 21:12
  • $\begingroup$ @eldo Can you show me an instance where it produces the wrong answer? min[1000000] // IntegerDigits // #^2 & // Total == 1000009 $\endgroup$
    – mfvonh
    Jun 18, 2014 at 21:18
  • $\begingroup$ @mfvonh I think it was a joke.. $\endgroup$
    – Öskå
    Jun 18, 2014 at 21:23
  • $\begingroup$ @mfvonh - I'm sorry, but I don't see the answer on my screen for min[1000000] // AbsoluteTiming. I just see the timing and after that "...". An appeal court could reject the answer because of this :) $\endgroup$
    – eldo
    Jun 18, 2014 at 21:26
  • $\begingroup$ @eldo Oh I see, you just mean in the post? If so, I left the answer out because it is 12,346 digits long :P $\endgroup$
    – mfvonh
    Jun 18, 2014 at 21:29
5
$\begingroup$

I thinks this is short, easy, and fast solution:

x = 100;
n = FromDigits@
  Flatten@{Ceiling@Sqrt[x - 9^2 Floor@(x/9^2)], 
    ConstantArray[9, (Floor@(x/9^2))]}

59

Please try it and let me know.

$\endgroup$
8
  • $\begingroup$ @Mr.Wizard, what do you think? $\endgroup$ Jun 19, 2014 at 7:35
  • $\begingroup$ I consider your solution to be the best so far presented. It's even faster than mfvonh's method. I got 0.034 seconds for x = 10^7. Plus it seems to hold for all values >= 0. You can get an extra speed gain by defining w = Floor@(x/9^2) at the beginning :) $\endgroup$
    – eldo
    Jun 19, 2014 at 13:18
  • $\begingroup$ @eldo, thanks for your comments. $\endgroup$ Jun 19, 2014 at 15:53
  • $\begingroup$ My tentative victory vanishes :) $\endgroup$
    – mfvonh
    Jun 19, 2014 at 18:52
  • $\begingroup$ @mfvonh, I am always amazed by your knowledge of Mathematica. I wish I could be someday expert in Mathematica like you or like Mr. Wizard. thanks for your comment. $\endgroup$ Jun 19, 2014 at 19:36
3
$\begingroup$
    Minimize[{d0 + 10 d1 + 10^2 d2, d0^2 + d1^2 + d2^2 > 100, 
  d0 == 0 || d0 == 1 || d0 == 2 || d0 == 3 || d0 == 4 || d0 == 5 || 
   d0 == 6 || d0 == 7 || d0 == 8 || d0 == 9, 
  d1 == 0 || d1 == 1 || d1 == 2 || d1 == 3 || d1 == 4 || d1 == 5 || 
   d1 == 6 || d1 == 7 || d1 == 8 || d1 == 9, 
  d2 == 0 || d2 == 1 || d2 == 2 || d2 == 3 || d2 == 4 || d2 == 5 || 
   d2 == 6 || d2 == 7 || d2 == 8 || d2 == 9}, {d0, d1, d2}]
$\endgroup$
0
3
$\begingroup$

I'll propose two modifications to this problem [Read: I want to change it to make my proposed method look better]. First we up the constraint bound so that exhaustive search is more, well, exhausting. Then we weight the summands so that the method used by @Brian Megquier becomes more difficult to employ.

Specifically I raised the threshold from 100 to 10000 and I multiply the square of the jth digit by the jth prime. Again the objective is to find the smallest number subject to this threshold constraint. The method I used is 0-1 integer linear programming, mostly because it's all I know.

Here I set this up for a call to NMinimize. It recognizes ILPs and forks them over to some COIN-CLP library code.

squaresSumMin = 10000;
nvars = 15;
digits = Array[d, {nvars, 10}];
fvars = Flatten[digits];
c1 = Map[Total[#] == 1 &, digits];
c2 = Map[0 <= # <= 1 &, fvars];
c3 = {Prime[Range[nvars]].(digits.(Range[0, 9]^2)) >= squaresSumMin, 
   Element[fvars, Integers]};
obj = 10^Range[0, nvars - 1].(digits.Range[0, 9]);
constraints = Join[c1, c2, c3];

Timing[{min, vals} = NMinimize[{obj, constraints}, Flatten[digits]];]

(* Out[301]= {0.044000, Null} *)

Here is the winner.

Round[min]

(* Out[303]= 9899999989 *)

I actually checked this with a much slower run through (exact) Minimize, and it agrees after 20 minutes or so of deep soul searching.

A caveat is that this method will also have trouble should we go much higher with the threshold. I believe it is because the library code eventually has difficulty with machine precision arithmetic giving integers to close enough approximation. One could code an explicit branch-and-prune loop though I've not done so for this example.

$\endgroup$
2
$\begingroup$

Figure out how many nines,then fill the 1st number with the smallest non nine that works. This example is just to show the algorithm using pseudo code.

myInt=0
mySum=1000  \\or whatever
numNines=mySum/81

numNines=numNines-(numNines mod 1)  \\convert to int without rounding

for i=8 to 1; i--
  if (mySum mod 81) > i^2
    myInt=i

if numNines>=1
  for  j=0 to (numNines-1); j++
    myInt= (myInt*10) + 9
$\endgroup$
1
  • $\begingroup$ sure if you want to think about the problem instead of diving in and writing code (+1) $\endgroup$
    – george2079
    Jun 19, 2014 at 12:59
1
$\begingroup$

I could have written this as an oneliner, but I show it step for step:

x1 = Tuples[Range@9, 2];
x2 = First[#]^2 + Last[#]^2 & /@ x1;
p = First@Position[x2, x_ /; x > 100];
x1[[p]]

{{5, 9}}

Lots of possible improvements here, but it's late and I go.

$\endgroup$
1
$\begingroup$

If we do choose brute force this question is a duplicate of: Iterate until condition is met.

I would suggest starting with Select as I explained in my answer there. Also, try to vectorize operations (such as Power) as often as possible. Compare:

Select[Range@1*^6, Tr[IntegerDigits[#]^2] > 450 &, 1] // Timing
{1.544, {799999}}
myInt = 0;
mySum = 0;
SumLimit = 450;
While[mySum < SumLimit, myDigits = IntegerDigits[myInt];
  mySum = Total[#^2 & /@ myDigits];
  myInt = myInt + 1;] // Timing
myInt - 1
{5.975, Null}

799999

The range (Range@1*^6) was arbitrarily chosen; to make this open-ended I would use the block-based method I also described in that answer.

$\endgroup$
1
$\begingroup$

The idea here is among all the tuples with the same digits only the one sorted lowest needs to be checked. That sort/union operation is much faster than the total^2 operation so it is worthwhile:

 SumLimit = 450
 First@(FromDigits /@ 
       Select[Union@(Sort /@ Tuples[Range[0, 9], {Ceiling[SumLimit/81]}]),
           Total@(#^2) > SumLimit &])

{0.608404, 799999}

Faster still, here we avoid generating all the tuples to begin with.

 nd = Ceiling[SumLimit/81]+1;
 First@Sort[
     FromDigits /@ 
        Select[ #[[2]] & /@ 
           Nest[( Flatten[
              Map[(next = #[[1]] + 1;
                {next, #} & /@ ({#[[2]]}~Join~ Table[ Join[#[[2]] ,
                 ConstantArray[#[[1]], i ] ], {i, 1, 
                      nd - Length[#[[2]]]}])) & , #, 1], 1]) &  , {{0, {}}} , 10]  , 
                        Length[#] == nd && Total[#^2] > SumLimit&]] // Timing

{0.140401, 799999}

$\endgroup$
0
0
$\begingroup$

Did you try this? You may use the "select" built-in function

Mathematica Selecting

$\endgroup$
1
  • $\begingroup$ How is this different from the answer posted above by @Mr.Wizard? You just changed Tr to Total $\endgroup$
    – RunnyKine
    Jun 19, 2014 at 3:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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