Find the smallest positive integer that satisfies a certain condition

Question

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.

Answer 1

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

Answer 2

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.

Answer 3

    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}]

Answer 4

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.

Answer 5

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

Answer 6

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}

Answer 7

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.

Answer 8

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.

Answer 9

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