OFFSET
1,1
COMMENTS
The 2nd and 3rd repunit primes, 1111111111111111111 and 11111111111111111111111 are happy primes. - Thomas M. Green, Oct 23 2009
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
REFERENCES
R. K. Guy, Unsolved Problems Number Theory, Sect. E34.
LINKS
Nathaniel Johnston, Table of n, a(n) for n = 1..10000
Carlos Rivera, Puzzle 21. Happy primes, The Prime Puzzles and Problems Connection.
Eric Weisstein's World of Mathematics, Happy Number
Doctor Who, Episode 42
Wikipedia, Happy number
Wikipedia, Doctor Who, Episode 42
MATHEMATICA
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 *)
hpQ[p_]:=NestWhile[Total[IntegerDigits[#]^2]&, p, #!=1&, 1, 50]==1; Select[Prime[ Range[ 300]], hpQ] (* Harvey P. Dale, Jun 07 2022 *)
PROG
(PARI) has(n)=while(n>6, n=norml2(digits(n))); n==1
is(n)=has(n) && isprime(n) \\ Charles R Greathouse IV, Dec 14 2015
(Python)
from sympy import isprime
def swb(n): return sum(map(lambda x: x*x, map(int, str(n))))
def happy(bd):
while bd not in [1, 4]: bd = swb(bd) # iterate to fixed point or cycle
return bd == 1
def ok(n): return isprime(n) and happy(n)
def aupto(n): return [k for k in range(1, n+1) if ok(k)]
print(aupto(2012)) # Michael S. Branicky, Jul 13 2022
CROSSREFS
KEYWORD
nonn,easy,base
AUTHOR
EXTENSIONS
More terms from Patrick De Geest, Oct 15 1999
STATUS
approved