%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