A056652 - OEIS

%I #26 Sep 28 2020 01:03:51

%S 3,7,9,21,27,31,49,63,81,93,127,147,189,217,243,279,343,381,441,567,

%T 651,729,837,889,961,1029,1143,1323,1519,1701,1953,2187,2401,2511,

%U 2667,2883,3087,3429,3937,3969,4557,5103,5859,6223,6561,6727,7203,7533,8001,8191,8649,9261

%N Integers > 1 whose prime divisors are all Mersenne primes (i.e., of the form (2^p - 1)).

%H Amiram Eldar, <a href="/A056652/b056652.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..471 from Michael De Vlieger)

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

%e 63 is included because the prime factorization of 63 is 3^2 * 7 = (2^2 -1)^2 *(2^3 -1).

%p isA000668 := proc(n)

%p     if n in [   3, 7, 31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111, 162259276829213363391578010288127, 170141183460469231731687303715884105727] then

%p         true;

%p     else

%p         false;

%p     end if;

%p end proc:

%p isA056652 := proc(n)

%p     local p;

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

%p         if not isA000668(p) then

%p             return false;

%p         end if;

%p     end do:

%p     true ;

%p end proc:

%p for n from 2 to 1000 do

%p     if isA056652(n) then

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

%p     end if;

%p end do: # _R. J. Mathar_, Feb 19 2017

%t Block[{nn = 10^4, s}, s = TakeWhile[Select[2^Prime@ Range@ 8 - 1, PrimeQ], # <= nn &]; Select[Range@ nn, AllTrue[FactorInteger[#][[All, 1]], MemberQ[s, #] &] &]] (* _Michael De Vlieger_, Sep 03 2017 *)

%o (PARI) isok(n) = {if (n==1, return (0)); my(f = factor(n)); for (k=1, #f~, if (! ((q=ispower(f[k, 1]+1,,&e)) && isprime(q) && (e==2)), return(0));); 1;} \\ _Michel Marcus_, Apr 25 2016

%Y Cf. A000668, A046528.

%K nonn

%O 1,1

%A _Leroy Quet_, Aug 09 2000

%E Offset corrected and more terms added by _Michel Marcus_, Apr 25 2016