%I #34 Apr 07 2021 00:23:45

%S 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,

%T 35,36,40,42,44,45,48,49,50,54,55,56,60,63,64,66,70,72,75,77,80,81,84,

%U 88,90,96,98,99,100,101,105,108,110,112,120,121,125,126,128,131

%N Numbers all of whose prime factors are palindromes.

%C Multiplicative closure of A002385; A051038 and A046368 are subsequences. - _Reinhard Zumkeller_, Apr 11 2011

%H Ivan Neretin, <a href="/A033620/b033620.txt">Table of n, a(n) for n = 1..10550</a>

%H <a href="/index/Pri#prime_factors">Index entries sequences related to prime factors</a>

%F Sum_{n>=1} 1/a(n) = Product_{p in A002385} p/(p-1) = 5.0949... - _Amiram Eldar_, Sep 27 2020

%e 10 = 2 * 5 is a term since both 2 and 5 are palindromes.

%e 110 = 2 * 5 * 11 is a term since 2, 5 and 11 are palindromes.

%p N:= 5: # to get all terms of up to N digits

%p digrev:= proc(t) local L; L:= convert(t,base,10);

%p add(L[-i-1]*10^i,i=0..nops(L)-1);

%p end proc:

%p PPrimes:= [2,3,5,7,11]:

%p for d from 3 to N by 2 do

%p m:= (d-1)/2;

%p 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)]);

%p od:

%p PPrimes:= map(op,[PPrimes]):

%p M:= 10^N:

%p B:= Vector(M);

%p B[1]:= 1:

%p for p in PPrimes do

%p for k from 1 to floor(log[p](M)) do

%p R:= [$1..floor(M/p^k)];

%p B[p^k*R] := B[p^k*R] + B[R]

%p od

%p od:

%p select(t -> B[t] > 0, [$1..M]); # _Robert Israel_, Jul 05 2015

%p # alternative

%p isA033620:= proc(n)

%p for d in numtheory[factorset](n) do

%p if not isA002113(op(1,d)) then

%p return false;

%p end if;

%p end do;

%p true ;

%p end proc:

%p for n from 1 to 300 do

%p if isA033620(n) then

%p printf("%d,",n) ;

%p end if;

%p end do: # _R. J. Mathar_, Sep 09 2015

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

%o (Haskell)

%o a033620 n = a033620_list !! (n-1)

%o a033620_list = filter chi [1..] where

%o chi n = a136522 spf == 1 && (n' == 1 || chi n') where

%o n' = n `div` spf

%o spf = a020639 n -- cf. A020639

%o -- _Reinhard Zumkeller_, Apr 11 2011

%o (PARI) ispal(n)=n=digits(n);for(i=1,#n\2,if(n[i]!=n[#n+1-i],return(0)));1

%o is(n)=if(n<13,n>0,vecmin(apply(ispal,factor(n)[,1]))) \\ _Charles R Greathouse IV_, Feb 06 2013

%o (Python)

%o from sympy import isprime, primefactors

%o def pal(n): s = str(n); return s == s[::-1]

%o def ok(n): return all(pal(f) for f in primefactors(n))

%o print(list(filter(ok, range(1, 132)))) # _Michael S. Branicky_, Apr 06 2021

%Y Cf. A002113, A002385, A046368, A051038.

%K nonn,base,easy

%O 1,2

%A _N. J. A. Sloane_, May 17 1998