A364479 Happy palindromic primes.

7, 313, 383, 11311, 15451, 30103, 30803, 35053, 36263, 71317, 74047, 94349, 94649, 95959, 98689, 1221221, 1257521, 1262621, 1281821, 1311131, 1444441, 1551551, 1594951, 1597951, 1658561, 1703071, 1737371, 1764671, 1829281, 1924291, 1957591, 1970791, 1981891, 1988891, 3001003

Offset: 1

Author

Simon R Blow, Jul 26 2023

Keywords

• nonn
• base

Comments

All terms in the sequence are prime numbers that read the same backward and forward. Each term is a happy number, meaning it converges to 1 under the process of repeatedly summing the squares of its digits. The sequence is a subset of both the palindromic prime numbers (A002385) and the happy numbers (A007770).

There are no happy palindromic primes with an even number of digits: every palindromic number with an even number of digits is divisible by 11, so 11 itself is the only palindromic prime, and it is not a happy number.

Examples

			313 is a term as it is palindromic (can be reversed), is a prime and is happy: 3^2 + 1^2 + 3^2 = 19, 1^2 + 9^2 = 82, 8^2 + 2^2 = 68, 6^2 + 8^2 = 100, 1^2 + 0^2 + 0^2 = 1.
		

Links

• Chai Wah Wu, Table of n, a(n) for n = 1..10000
• Eric Weisstein's World of Mathematics, Happy Number.
• Eric Weisstein's World of Mathematics, Palindromic Prime.

Crossrefs

Intersection of A002385 and A007770.

Programs

Mathematica

happyQ[n_] := NestWhile[Plus @@ (IntegerDigits[#]^2) &, n, UnsameQ, All] == 1; Select[Prime[Range[220000]], PalindromeQ[#] && happyQ[#] &] (* Amiram Eldar, Jul 28 2023 *)

Python

def is_prime(num):
    return num > 1 and all(num % i != 0 for i in range(2, int(num ** 0.5) + 1))
def is_palindrome(num):
    return str(num) == str(num)[::-1]
def is_happy(num):
    while num != 1 and num != 4:
        num = sum(int(digit) ** 2 for digit in str(num))
    return num == 1
happy_palindromic_primes = [num for num in range(1, 10000000) if is_prime(num) and is_palindrome(num) and is_happy(num)]
print(happy_palindromic_primes)

Python

from itertools import islice
from sympy import isprime
def A364479_gen(): # generator of terms
    n = 1
    while True:
        for z in (1,3,5,7,9):
            for y in range(z*n, (z+1)*n):
                k, m = y//10, 0
                while k >= 10:
                    k, r = divmod(k, 10)
                    m = 10*m + r
                if isprime(a:=y*n + 10*m + k):
                    b = a
                    while b not in {1,37,58,89,145,42,20,4,16}:
                        b = sum((0, 1, 4, 9, 16, 25, 36, 49, 64, 81)[ord(d)-48] for d in str(b))
                    if b == 1:
                        yield a
        n *= 10
A364479_list = list(islice(A364479_gen(),20)) # Chai Wah Wu, Aug 02 2023