41
$\begingroup$

This being the Q number 2K in the site, and this being the day we got the confirmation of mathematica.se graduating soon, I think a celebration question is in order.

So...

What is the fastest way to compute the happy prime number 2000 in Mathematica?

Edit

Here are the Timing[ ] results so far:

 {"JM",      {5.610, 137653}}
 {"Leonid",  {5.109, {12814, 137653}}}
 {"wxffles", {4.11, {12814, 137653}}}
 {"Rojo",    {0.765, {12814, 137653}}}
 {"Rojo1",   {0.547, {12814, 137653}}}
Improve this question
$\endgroup$
7
  • 8
    $\begingroup$ You should offer a 2K bounty to celebrate $\endgroup$ – Rojo Jun 20 '12 at 23:17
  • $\begingroup$ @Rojo A bounty can't be started until a few days after the Q was posted $\endgroup$ – Dr. belisarius Jun 20 '12 at 23:31
  • 2
    $\begingroup$ We can extend the celebration :) $\endgroup$ – Rojo Jun 20 '12 at 23:32
  • $\begingroup$ You tested the boringer simpler version of my edit? $\endgroup$ – Rojo Jun 21 '12 at 0:41
  • $\begingroup$ Not yet, testing right now $\endgroup$ – Dr. belisarius Jun 21 '12 at 0:44

10 Answers 10

Active Oldest Votes
25
$\begingroup$

This answer should be read upside down, since the last edit has the fastest, neatest and shortest answer

Module[{$guard = True},

happyQ[i_] /; $guard := Block[{$guard = False, appeared},
   appeared[_] = False;
   happyQ[i]
   ]
 ]

e : happyQ[_?appeared] := e = False;

happyQ[1] = True;

e : happyQ[i_] := e = (appeared[i] = True; happyQ[#.#&@IntegerDigits[i]])

Now, taking this from @LeonidShiffrin

happyPrimeN[n_] := Module[{m = 0, pctr = 0},
   While[m < n, If[happyQ@Prime[++pctr], m++]];
   {pctr, Prime[pctr]}];

EDIT

Ok, this was cool, but if you don't mind wasting a little memory and not resetting appeared, it becomes simple and less cool

appeared[_] = False;
happyQ[1] = True;
happyQ[_?appeared] = False;
e : happyQ[i_] := e = (appeared[i] = True; happyQ[#.# &@IntegerDigits[i]])

EDIT2

Slightly faster but I like it twice as much

happyQ[1] = True;
e : happyQ[i_] := (e = False; e = happyQ[#.# &@IntegerDigits[i]])

or perhaps to make it slightly shorter and a little bit more memory efficient, reducing the recursion tree's height

happyQ[1] = True;
e : happyQ[i_] := e = happyQ[e = False; #.# &@IntegerDigits[i]]
Improve this answer
$\endgroup$
10
  • $\begingroup$ Chagracia! Para vos también @belisarius. +1 to the question! $\endgroup$ – Rojo Jun 20 '12 at 23:57
  • 5
    $\begingroup$ Speedy! And in the best traditions, obfuscated too. $\endgroup$ – wxffles Jun 21 '12 at 0:03
  • 1
    $\begingroup$ This is very cool!+1 $\endgroup$ – Leonid Shifrin Jun 21 '12 at 0:07
  • 1
    $\begingroup$ You got an extra ^2 after IntegerDigits in your first solution. I just realized my solution is exactly the same as your except it's written in a much more boring way. $\endgroup$ – Szabolcs Jun 21 '12 at 14:50
  • $\begingroup$ Edit 2 is awesome! $\endgroup$ – Simon Woods Jun 29 '12 at 21:59
13
$\begingroup$

This probably counts as cheating, since it uses the fact that all unhappy numbers end up in a cycle including 4. But I like it for simplicity...

happyQ[1]=True;
happyQ[4]=False;
happyQ[n_]:=happyQ[n]=happyQ[#.#&@IntegerDigits[n]]

This works with Leonid's happyPrimeN function.

Improve this answer
$\endgroup$
4
  • $\begingroup$ +1 and happy 2K, Simon. I'll include the timing later (need to work now) $\endgroup$ – Dr. belisarius Jun 21 '12 at 13:56
  • $\begingroup$ Can you prove that this is the case? EDIT There's a proof on Wikipedia. $\endgroup$ – Szabolcs Jun 21 '12 at 14:20
  • $\begingroup$ @Szabolcs: there's a proof here. $\endgroup$ – J. M.'s torpor♦ Jun 21 '12 at 14:23
  • $\begingroup$ Hehe, can't fight knowledge with brute force. +1! $\endgroup$ – Rojo Jun 22 '12 at 13:06
12
$\begingroup$

Simple top-level solution

Here is a simplistic completely top-level code:

ClearAll[happyQ];
happyQ[n_] :=
  Block[{appeared},
    appeared[_] = False;
    Take[
       NestWhileList[
          Total[IntegerDigits[#]^2] &,
          n,
          (! appeared[#] && (appeared[#] = True)) &
       ], -2] == {1, 1}];

Clear[happyPrimeN];
happyPrimeN[n_] :=
  Module[{m = 0, pctr = 0},
    While[m < n, If[happyQ@Prime[++pctr], m++]];
    {pctr, Prime[pctr]}
  ];

Using this, we get for example:

happyPrimeN/@Range[5]

(* {{4,7},{6,13},{8,19},{9,23},{11,31}}  *)

And for 2000th happy prime, we have:

happyPrimeN[2000] // AbsoluteTiming

(* {1.5693359, {12814, 137653}}  *)

which is not particularly fast, but probably ok. I am sure that there are faster solutions though.

Java solution with memoization

One thing I want to mention here: I had about 10 iterations of this one before I finally optimized it, and when I did, I looked closer at @Rojo's solution and found that I just arrived to a Java port of it. So, while I did it independently, I just want to stress that the following code does not contain new or better ideas than those used by @Rojo for his beautiful solution.

Ok, so:

  1. Load the Java reloader

  2. Compile the following class:

    JCompileLoad@"import java.util.*;

    public class HappyPrimes{   
        public Map<Integer,Boolean> happy = new HashMap<Integer,Boolean>(10000);
        private int max;

        public HappyPrimes(int max){
             this.max = max;
             happy.put(1,true); 
        }

        public  int getDigitsSqSum(int num){
            int result = 0;
            while(num>0){
                int dig = num % 10;
                result+=dig*dig;
                num /=10;
            }       
            return  result;
        }

        private boolean isHappy(int num){
            if(happy.containsKey(num)){
                return happy.get(num);
            }
            happy.put(num,false);
            boolean result  = isHappy(getDigitsSqSum(num));
            happy.put(num,result);
            return result;
        }

        public  int[] currentMaxHappyPrime(int[] primes,
                           int startPrime, int currentMax){
            int done = 0;               
            int i = 0;      
            for( ; i< primes.length ; i++){     
                 if(isHappy(primes[i])&& ++currentMax == max){
                      done = 1;
                      break;
                 }                                          
            }
            startPrime+=i;
            return new int[]{startPrime,currentMax,done};
        } 
    }";

  3. The "top-level" function follows:

    ClearAll[happyPrimeNJ];
happyPrimeNJ[n_, chunk_: 5000] :=
  JavaBlock[
     With[{o = JavaNew["HappyPrimes", n]},
       {#, Prime[#]} &@(First[#] + 1) &@
          NestWhile[
             o@currentMaxHappyPrime[
                 Prime[Range[First@# + 1, First@# + chunk]], #[[1]], #[[2]]
             ] &,
             {0, 0, 0},
             Last@# != 1 &]
     ]
  ];

What happens here is that I use Mathematica to generate primes in chunks. I send those to Java and count the number of happy primes in a given chunk. When I get enough, I stop and return the prime index. At intermediate steps, I return a list of 3 numbers: current total number of processed primes, current total number of happy primes among those, and a flag telling me whether or not I should continue.

Here is how we use it:

happyPrimeNJ[50000]//AbsoluteTiming

{1.2324219,{365523,5263169}}

My benchmarks show that it is systematically several times (up to 10) faster than @Rojo's version, but we don't see a dramatic speed-up as in some other cases, since @Rojo very cleverly uses the language, and Mathematica hash tables (DownValues) are pretty efficient. Also, for (relatively) small numbers of happy primes (such as 2000), the speedup is not so apparent since there is a constant overhead of Java calls which is of the same order as the total time it takes to process those.

Summary and conclusions

The first method I presented is relatively slow, being the most straightforward. The second one, based on Java, is fast. However, it does not really compete with the elegant and terse solution of @Rojo, and moreover, is more or less a direct port of it to Java (even though I arrived at it mostly independently).

Improve this answer
$\endgroup$
6
  • $\begingroup$ +1 happy 2K, Leonid! As more answers come in, I'll publish benchmarking figures from my machine. $\endgroup$ – Dr. belisarius Jun 20 '12 at 22:58
  • 1
    $\begingroup$ @belisarius Thanks, and happy 2K for you too! And, I am not done with this yet :) $\endgroup$ – Leonid Shifrin Jun 20 '12 at 22:59
  • $\begingroup$ Due to some problems in my JVM setup, I can't test your last version. As this is now CW, you may go on and edit the scoring section of the Q with the relevant info. Sorry! $\endgroup$ – Dr. belisarius Jun 21 '12 at 2:11
  • $\begingroup$ @belisarius No problem :) I will edit the scoring section when I get some time. B.t.w., what platform are you in? The reloader should work out of the box, because it uses the JVM which comes bundled with Mathematica. So if it doesn't, I want to know about it, so I can see what the problem is, and fix it. $\endgroup$ – Leonid Shifrin Jun 21 '12 at 14:22
  • $\begingroup$ My problem with the JVM is just because I mixed up the bundled one with my own while playing a little. Nothing serious, but I have to reinstall a few things before both of them can work properly :( $\endgroup$ – Dr. belisarius Jun 21 '12 at 14:32
11
$\begingroup$ 
Clear[happyPrimeN];
happyPrimeN[2000] = 137653;

happyPrimeN[2000] // AbsoluteTiming

{0., 137653}

But seriously, here's a memoised, recursive happyQ that can be used with Leonid's happyPrimeN

Clear[sos, happyQ];
sos[k_Integer] := sos[k] = #.# &[IntegerDigits[k]];
happyQ[k_Integer] := happyQ[k] = happyQ[k, {}];
happyQ[1, history_List] := True;
happyQ[k_Integer, history_List] := 
   With[{h = sos[k]}, If[MemberQ[history, h], False, happyQ[h, Prepend[history, h]]]];

happyPrimeN[2000] // AbsoluteTiming
happyPrimeN[2000] // AbsoluteTiming

{1.4531250, {12814, 137653}}

{0.0468750, {12814, 137653}}

Improve this answer
$\endgroup$
5
  • $\begingroup$ Note: Total[IntegerDigits[#]^2] is faster than my silly Plus @@ (#^2 & /@ IntegerDigits[k]), but I can't steal all of Leonid's code. I always forget about Listable. $\endgroup$ – wxffles Jun 20 '12 at 23:52
  • $\begingroup$ C'mon, don't be ashamed! This is for fun! $\endgroup$ – Dr. belisarius Jun 20 '12 at 23:54
  • 1
    $\begingroup$ If you really wanted an alternative: sos[k_Integer] := Norm[IntegerDigits[k]]^2. $\endgroup$ – J. M.'s torpor♦ Jun 21 '12 at 8:39
  • 1
    $\begingroup$ You just made me look for an alternative to Total[IntegerDigits[#]^2], and since what I found seems to be faster than that one, I used it myself, hehe. But I hereby, on the record, grant you the exclusive license to use it #.#&[IntegerDigits[#]] $\endgroup$ – Rojo Jun 21 '12 at 12:54
  • $\begingroup$ @Rojo I was wondering when someone would get around to #.#&[..] $\endgroup$ – Brett Champion Jun 21 '12 at 15:27
10
$\begingroup$

Here's my take:

(* Brent's algorithm for cycle detection *)
happyQ[start_Integer] := Module[{cyc, f, hare, pow, tortoise},
       f = Total[IntegerDigits[#]^2] &;
       cyc = pow = 1;
       tortoise = start; hare = f[start];
       While[tortoise =!= hare,
             If[pow == cyc,
                tortoise = hare; pow *= 2; cyc = 0;];
             hare = f[hare];
             cyc++];
       cyc === 1]

happyPrimeN[1] = 7;
happyPrimeN[n_Integer] := happyPrimeN[n] = 
  Block[{$RecursionLimit = Infinity}, 
   NestWhile[NextPrime, happyPrimeN[n - 1], (! happyQ[#] &), {2, 1}]]
Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ @bel: confetti and alcohol is in order indeed. :) $\endgroup$ – J. M.'s torpor♦ Jun 21 '12 at 0:44
9
$\begingroup$

Here is my method (not sure it counts):

happyPrimeN = Import["http://oeis.org/A035497/b035497.txt", "table"][[#, 2]] &

happyPrimeN[2000] // AbsoluteTiming

Out[14]= {0.7490428, 137653}
Improve this answer
$\endgroup$
1
  • $\begingroup$ Hehe ... I have a veeeery slow Internet connection :D $\endgroup$ – Dr. belisarius Jun 21 '12 at 2:54
6
$\begingroup$

a short functional style solution:

HappyQ[n_Integer?Positive] := NestWhile[
    Total[IntegerDigits[#]^2] &, n,
    Unequal,
    All
  ] == 1
NextHappyPrime[n_Integer?Positive] := NestWhile[
    NextPrime,
    NextPrime[n],
    Composition[Not, HappyQ]
  ]
HappyPrimeN[n_Integer?Positive] := Nest[NextHappyPrime, 7, n - 1]

HappyPrimeN[2000]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Elegant! +1 and happy 2K whatever that means $\endgroup$ – Rojo Jun 29 '12 at 21:06
5
$\begingroup$

EDIT: After reading the other answers, this seems to be the same as Rojo's method, except written in a less interesting way.

Here's my shot at the problem. I didn't look at the other solutions to keep this more fun (it's community wiki anyway).

happyQ;

Begin["happyQ`"];

happyQ[num_] :=
 Block[{seenQ},
  seenQ[_] = False;
  isHappy[num]
  ]

propagate = Total[IntegerDigits[#]^2] &;

isHappy[1] = True;

isHappy[num_] :=
 If[seenQ[num],
  False,
  seenQ[num] = True;
  isHappy[num] = isHappy[propagate[num]]
  ]

End[];


happyPrimeN[n_Integer] :=
 Block[{$RecursionLimit = Infinity, count, p},
  count = 0;
  p = 1;
  While[count < 2000,
   p = NextPrime[p];
   While[! happyQ[p], p = NextPrime[p]];
   count++
   ];
  p
 ]

I get a timing of 0.8 seconds on this machine.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Okay, looking at other solutions, nothing new here. But it was fun :-) $\endgroup$ – Szabolcs Jun 21 '12 at 14:20
  • $\begingroup$ +1, interesting scoping. $\endgroup$ – rcollyer Jun 21 '12 at 14:21
  • $\begingroup$ No matter how it is written, it is a cool method. +1. $\endgroup$ – Leonid Shifrin Jun 21 '12 at 16:24
3
$\begingroup$

I'm pretty late to the party but here's my attempt:

It involves grabbing a dynamic amount of primes at a time and testing.

happyQ[1] = True;
happyQ[4] = False;
happyQ[n_] := happyQ[n] = happyQ[Total[IntegerDigits[n]^2]]

happyPrime[n_] := iHappyPrime[n, 0]

iHappyPrime[n_, s_] := With[{primes = Select[Range[s, NextPrime[s, n]], PrimeQ]},
    With[{count = Count[primes, v_ /; happyQ[v]]},
        If[count < n,
            iHappyPrime[n - count, Last[primes] + 1],
            Select[primes, happyQ][[n]]
        ]
    ]
]

In[1]:= happyPrime[2000] // Timing

Out[1]= {0.479137, 137653}
Improve this answer
$\endgroup$
1
$\begingroup$

I'm very late here, but why not?

First I steal the list of happy numbers below 1000 from wikipedia:

happy1000 = {1, 7, 10, 13, 19, 23, 28, 31, 32, 44, 49, 68, 70, 79, 82,
   86, 91, 94, 97, 100, 103, 109, 129, 130, 133, 139, 167, 176, 188, 
  190, 192, 193, 203, 208, 219, 226, 230, 236, 239, 262, 263, 280, 
  291, 293, 301, 302, 310, 313, 319, 320, 326, 329, 331, 338, 356, 
  362, 365, 367, 368, 376, 379, 383, 386, 391, 392, 397, 404, 409, 
  440, 446, 464, 469, 478, 487, 490, 496, 536, 556, 563, 565, 566, 
  608, 617, 622, 623, 632, 635, 637, 638, 644, 649, 653, 655, 656, 
  665, 671, 673, 680, 683, 694, 700, 709, 716, 736, 739, 748, 761, 
  763, 784, 790, 793, 802, 806, 818, 820, 833, 836, 847, 860, 863, 
  874, 881, 888, 899, 901, 904, 907, 910, 912, 913, 921, 923, 931, 
  932, 937, 940, 946, 964, 970, 973, 989, 998, 1000};

Then I set DownValues for happyQ

Do[happyQ[n] = True, {n, happy1000}];
Do[happyQ[n] = False, {n, Complement[Range[1000], happy1000]}];
happyQ[n_Integer] := happyQ[#.#&@IntegerDigits[n]]

And finally I output the list of 2 and the first 2000 happy primes:

NestList[NestWhile[NextPrime, NextPrime[#], Not@*happyQ] &, 2, 2000]

or more briefly I do

Nest[NestWhile[NextPrime, NextPrime[#], Not@*happyQ] &, 2, 2000] // AbsoluteTiming

{0.216036, 137653}

Undoubtedly, since three years ago execution speed must have improved.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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