|
%I #43 Jul 13 2022 03:42:13
%S 7,13,19,23,31,79,97,103,109,139,167,193,239,263,293,313,331,367,379,
%T 383,397,409,487,563,617,653,673,683,709,739,761,863,881,907,937,1009,
%U 1033,1039,1093,1151,1277,1303,1373,1427,1447,1481,1487,1511,1607,1663
%N Happy primes: primes that eventually reach 1 under iteration of "x -> sum of squares of digits of x".
%C The 2nd and 3rd repunit primes, 1111111111111111111 and 11111111111111111111111 are happy primes. - _Thomas M. Green_, Oct 23 2009
%C There are 200 terms up to 10^4, 1465 up to 10^5, 11144 up to 10^6, 91323 up to 10^7, 812371 up to 10^8, 7408754 up to 10^9, and 67982202 up to 10^10. These are consistent with b*prime(n) < a(n) < c*prime(n) with constants 0 < b < c. - _Charles R Greathouse IV_, Jan 06 2016
%D R. K. Guy, Unsolved Problems Number Theory, Sect. E34.
%H Nathaniel Johnston, <a href="/A035497/b035497.txt">Table of n, a(n) for n = 1..10000</a>
%H Carlos Rivera, <a href="http://www.primepuzzles.net/puzzles/puzz_021.htm">Puzzle 21. Happy primes</a>, The Prime Puzzles and Problems Connection.
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/HappyNumber.html">Happy Number</a>
%H Doctor Who, <a href="http://www.youtube.com/watch?v=ee2If8jSxUo">Episode 42</a>
%H Wikipedia, <a href="http://en.wikipedia.org/wiki/Happy_number">Happy number</a>
%H Wikipedia, <a href="http://en.wikipedia.org/wiki/42_(Doctor_Who)">Doctor Who, Episode 42</a>
%t g[n_] := Total[ IntegerDigits[n]^2]; fQ[n_] := NestWhileList[g@# &, n, UnsameQ, All][[-1]] == 1; Select[Prime@ Range@ 300, fQ@# &] (* _Robert G. Wilson v_, Jan 03 2013 *)
%t hpQ[p_]:=NestWhile[Total[IntegerDigits[#]^2]&,p,#!=1&,1,50]==1; Select[Prime[ Range[ 300]],hpQ] (* _Harvey P. Dale_, Jun 07 2022 *)
%o (PARI) has(n)=while(n>6, n=norml2(digits(n))); n==1
%o is(n)=has(n) && isprime(n) \\ _Charles R Greathouse IV_, Dec 14 2015
%o (Python)
%o from sympy import isprime
%o def swb(n): return sum(map(lambda x: x*x, map(int, str(n))))
%o def happy(bd):
%o while bd not in [1, 4]: bd = swb(bd) # iterate to fixed point or cycle
%o return bd == 1
%o def ok(n): return isprime(n) and happy(n)
%o def aupto(n): return [k for k in range(1, n+1) if ok(k)]
%o print(aupto(2012)) # _Michael S. Branicky_, Jul 13 2022
%Y Cf. A007770 (happy numbers), A046519.
%K nonn,easy,base
%O 1,1
%A _N. J. A. Sloane_
%E More terms from _Patrick De Geest_, Oct 15 1999
|
