%I #22 Jun 02 2017 22:25:28
%S 1093,3511,1194649,3837523,12327121,1305751357,4194412639,13473543253,
%T 43280521831,1427186233201,4584493014427,14726582775529,
%U 47305610361283,151957912148641,5010850864768711,16096154973653197,51705032124882319,166089997978464613
%N Numbers n > 1 where all prime factors are Wieferich primes, i.e., terms of A001220.
%C The prime terms are Wieferich primes.
%C All "Wieferich pseudoprimes", if any exist, are in the sequence (see second comment in A240719).
%H PrimeGrid, <a href="http://prpnet.primegrid.com:13000/">Wieferich prime search - PRPNet server statistics</a>
%e 4194412639 = 1093^2 * 3511. All prime factors are Wieferich primes, so 4194412639 is a term of the sequence.
%t Take[#, 19] &@ Rest@ Sort@ Map[1093^First@ # 3511^Last@ # &, Tuples[Range[0, 6], 2]] (* _Michael De Vlieger_, Mar 24 2016 *)
%o (PARI) is(n) = if(n==1, return(0)); my(f=factor(n)[, 1]); for(k=1, #f, if(Mod(2, f[k]^2)^(f[k]-1)!=1, return(0))); return(1)
%o (PARI) /* The following program is significantly faster; valid up to (p^x * q^y) < b, where b is the upper search bound for Wieferich primes (approximately 5*10^17 as of Mar 23 2016, see PrimeGrid PRPNet server statistics) */
%o my(p=1093, q=3511, v=vector(0), w=vector(1)); for(x=0, 4, for(y=0, 4, w[1]=p^x*q^y; v=concat(v, w))); vecextract(vecsort(v,,8), "2..25")
%Y Cf. A001220, A240719.
%K nonn
%O 1,1
%A _Felix Fröhlich_, Mar 23 2016