# Find the smallest positive integer that satisfies a certain condition

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


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.

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

I think this works correctly:

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

min // AbsoluteTiming


{0.001005, 59}

min // AbsoluteTiming


{0.024029, .... }

• @mfvonh - I agree with the "59", but I don't agree with the "..."
– eldo
Jun 18 '14 at 21:12
• @eldo Can you show me an instance where it produces the wrong answer? min // IntegerDigits // #^2 & // Total == 1000009 Jun 18 '14 at 21:18
• @mfvonh I think it was a joke..
– Öskå
Jun 18 '14 at 21:23
• @mfvonh - I'm sorry, but I don't see the answer on my screen for min // AbsoluteTiming. I just see the timing and after that "...". An appeal court could reject the answer because of this :)
– eldo
Jun 18 '14 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 Jun 18 '14 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.

• @Mr.Wizard, what do you think? Jun 19 '14 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 '14 at 13:18
• My tentative victory vanishes :) Jun 19 '14 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. Jun 19 '14 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}]


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= {0.044000, Null} *)


Here is the winner.

Round[min]

(* Out= 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.

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

• sure if you want to think about the problem instead of diving in and writing code (+1) Jun 19 '14 at 12:59

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.

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.

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[ #[] & /@
Nest[( Flatten[
Map[(next = #[] + 1;
{next, #} & /@ ({#[]}~Join~ Table[ Join[#[] ,
ConstantArray[#[], i ] ], {i, 1,
nd - Length[#[]]}])) & , #, 1], 1]) &  , {{0, {}}} , 10]  ,
Length[#] == nd && Total[#^2] > SumLimit&]] // Timing


{0.140401, 799999}

Did you try this? You may use the "select" built-in function • How is this different from the answer posted above by @Mr.Wizard? You just changed Tr to Total` Jun 19 '14 at 3:38