Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
No, you can't use Range[100]. What would you use if I said 10^5 instead of 10^2? –  cartonn Jun 18 at 20:18
    
just tightening it up helps a little: ` While[Total@(IntegerDigits[++myInt]^2) < SumLimit] –  george2079 Jun 18 at 20:35
    
@george2079 I like to have individual statements for the readability of the code. –  cartonn Jun 18 at 20:37
1  
@Ö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. –  RunnyKine Jun 18 at 21:03
    
Look at a related problem: Finding the largest integer that cannot be partitioned in a certain way –  Artes Jun 18 at 23:07

9 Answers 9

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, .... }

share|improve this answer
1  
@mfvonh - I agree with the "59", but I don't agree with the "..." –  eldo Jun 18 at 21:12
    
@eldo Can you show me an instance where it produces the wrong answer? min[1000000] // IntegerDigits // #^2 & // Total == 1000009 –  mfvonh Jun 18 at 21:18
    
@mfvonh I think it was a joke.. –  Öskå Jun 18 at 21:23
    
@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 :) –  eldo Jun 18 at 21:26
    
@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 –  mfvonh Jun 18 at 21:29

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.

share|improve this answer
    
@Mr.Wizard, what do you think? –  Algohi Jun 19 at 7:35
    
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 :) –  eldo Jun 19 at 13:18
    
@eldo, thanks for your comments. –  Algohi Jun 19 at 15:53
    
My tentative victory vanishes :) –  mfvonh Jun 19 at 18:52
    
@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. –  Algohi Jun 19 at 19:36
    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}]
share|improve this answer

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
share|improve this answer
    
sure if you want to think about the problem instead of diving in and writing code (+1) –  george2079 Jun 19 at 12:59

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.

share|improve this answer

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.

share|improve this answer

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.

share|improve this answer

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}

share|improve this answer

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

Mathematica Selecting

share|improve this answer
    
How is this different from the answer posted above by @Mr.Wizard? You just changed Tr to Total –  RunnyKine Jun 19 at 3:38

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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