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

-
@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[1000000] // IntegerDigits // #^2 & // Total == 1000009 – mfvonh 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[1000000] // 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 – mfvonh 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? – Algohi 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
@eldo, thanks for your comments. – Algohi Jun 19 '14 at 15:53
My tentative victory vanishes :) – mfvonh 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. – Algohi 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}]
-

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) – george2079 Jun 19 '14 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.

-

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[ #[[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}

-

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 – RunnyKine Jun 19 '14 at 3:38