A033620 Numbers all of whose prime factors are palindromes.

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 27, 28, 30, 32, 33, 35, 36, 40, 42, 44, 45, 48, 49, 50, 54, 55, 56, 60, 63, 64, 66, 70, 72, 75, 77, 80, 81, 84, 88, 90, 96, 98, 99, 100, 101, 105, 108, 110, 112, 120, 121, 125, 126, 128, 131

Offset: 1

Author

N. J. A. Sloane, May 17 1998

Keywords

• nonn
• base
• easy

Comments

Multiplicative closure of A002385; A051038 and A046368 are subsequences. - Reinhard Zumkeller, Apr 11 2011

Examples

			10 = 2 * 5 is a term since both 2 and 5 are palindromes.
110 = 2 * 5 * 11 is a term since 2, 5 and 11 are palindromes.
		

Links

• Ivan Neretin, Table of n, a(n) for n = 1..10550
• Index entries sequences related to prime factors

Crossrefs

Cf. A002113, A002385, A046368, A051038.

Programs

Haskell

a033620 n = a033620_list !! (n-1)
a033620_list = filter chi [1..] where
   chi n = a136522 spf == 1 && (n' == 1 || chi n') where
      n' = n `div` spf
      spf = a020639 n  -- cf. A020639
-- Reinhard Zumkeller, Apr 11 2011

Maple

N:= 5: # to get all terms of up to N digits
digrev:= proc(t) local L; L:= convert(t,base,10);
add(L[-i-1]*10^i,i=0..nops(L)-1);
end proc:
PPrimes:= [2,3,5,7,11]:
for d from 3 to N by 2 do
    m:= (d-1)/2;
    PPrimes:= PPrimes, select(isprime,[seq(seq(n*10^(m+1)+y*10^m+digrev(n), y=0..9), n=10^(m-1)..10^m-1)]);
od:
PPrimes:= map(op,[PPrimes]):
M:= 10^N:
B:= Vector(M);
B[1]:= 1:
for p in PPrimes do
  for k from 1 to floor(log[p](M)) do
     R:= [$1..floor(M/p^k)];
     B[p^k*R] := B[p^k*R] + B[R]
  od
od:
select(t -> B[t] > 0, [$1..M]); # Robert Israel, Jul 05 2015
# alternative
isA033620:= proc(n)
    for d in numtheory[factorset](n) do
        if not isA002113(op(1,d)) then
            return false;
        end if;
    end do;
    true ;
end proc:
for n from 1 to 300 do
    if isA033620(n) then
        printf("%d,",n) ;
    end if;
end do: # R. J. Mathar, Sep 09 2015

Mathematica

palQ[n_]:=Reverse[x=IntegerDigits[n]]==x; Select[Range[131],And@@palQ/@First/@FactorInteger[#]&] (* Jayanta Basu, Jun 05 2013 *)

PARI

ispal(n)=n=digits(n);for(i=1,#n\2,if(n[i]!=n[#n+1-i],return(0)));1
is(n)=if(n<13,n>0,vecmin(apply(ispal,factor(n)[,1]))) \\ Charles R Greathouse IV, Feb 06 2013

Python

from sympy import isprime, primefactors
def pal(n): s = str(n); return s == s[::-1]
def ok(n): return all(pal(f) for f in primefactors(n))
print(list(filter(ok, range(1, 132)))) # Michael S. Branicky, Apr 06 2021

Formula

Sum_{n>=1} 1/a(n) = Product_{p in A002385} p/(p-1) = 5.0949... - Amiram Eldar, Sep 27 2020