A057612 Numbers that are both Mersenne numbers (A001348) and lucky numbers (A000959).
3, 7, 31, 127, 8191, 131071, 524287, 8388607
Offset: 1
Extensions
More terms from Mark Weston (mweston(AT)uvic.ca), Oct 10 2001
This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
seq(numtheory:-order(10, 2^ithprime(i)-1),i=1..20); # Robert Israel, May 25 2020
f[n_] := Block[{ds = Divisors[n - 1]}, p = Position[ PowerMod[10, ds, n], 1]; If[p == {}, Length[ RealDigits[1/n][[1, 1]]], Take[ds, p[[1, 1]]][[ -1]]]]; Table[ f[2^Prime[n] - 1], {n, 11}]
Corresponding to the initial terms 2, 3, 5, 7, 13, 17, 19, 31 ... we get the Mersenne primes 2^2 - 1 = 3, 2^3 - 1 = 7, 2^5 - 1 = 31, 127, 8191, 131071, 524287, 2147483647, ... (see A000668).
MersennePrimeExponent[Range[48]] (* Eric W. Weisstein, Jul 17 2017; updated Oct 21 2024 *)
isA000043(n) = isprime(2^n-1) \\ Michael B. Porter, Oct 28 2009
is(n)=my(h=Mod(2,2^n-1)); for(i=1, n-2, h=2*h^2-1); h==0||n==2 \\ Lucas-Lehmer test for exponent e. - Joerg Arndt, Jan 16 2011, and Charles R Greathouse IV, Jun 05 2013 forprime(e=2,5000,if(is(e),print1(e,", "))); /* terms < 5000 */
from sympy import isprime, prime
for n in range(1,100):
if isprime(2**prime(n)-1):
print(prime(n), end=', ') # Stefano Spezia, Dec 06 2018
A000668:=Filtered(List(Filtered([1..600], IsPrime),i->2^i-1),IsPrime); # Muniru A Asiru, Oct 01 2017
A000668 := proc(n) local i; i := 2^(ithprime(n))-1: if (isprime(i)) then return i fi: end: seq(A000668(n), n=1..31); # Jani Melik, Feb 09 2011 # Alternate: seq(numtheory:-mersenne([i]),i=1..26); # Robert Israel, Jul 13 2014
2^Array[MersennePrimeExponent, 18] - 1 (* Jean-François Alcover, Feb 17 2018, Mersenne primes with less than 1000 digits *) 2^MersennePrimeExponent[Range[18]] - 1 (* Eric W. Weisstein, Sep 04 2021 *)
forprime(p=2,1e5,if(ispseudoprime(2^p-1),print1(2^p-1", "))) \\ Charles R Greathouse IV, Jul 15 2011
LL(e) = my(n, h); n = 2^e-1; h = Mod(2, n); for (k=1, e-2, h=2*h*h-1); return(0==h) \\ after Joerg Arndt in A000043 forprime(p=1, , if(LL(p), print1(p, ", "))) \\ Felix Fröhlich, Feb 17 2018
from sympy import isprime, primerange print([2**n-1 for n in primerange(1, 1001) if isprime(2**n-1)]) # Karl V. Keller, Jr., Jul 16 2020
a(2)=12 because 2^12 - 1 = 4095 = 5*(3^2)*7*13 is divisible by a square.
[n: n in [1..250] | not IsSquarefree(2^n-1)]; // Vincenzo Librandi, Jul 14 2015
N:= 250:
B:= Vector(N):
for n from 1 to N do
if B[n] <> 1 then
F:= ifactors(2^n-1,easy)[2];
if max(seq(t[2],t=F)) > 1 or (hastype(F,symbol)
and not numtheory:-issqrfree(2^n-1)) then
B[[seq(n*k,k=1..floor(N/n))]]:= 1;
fi
fi;
od:
select(t -> B[t]=1, [$1..N]); # Robert Israel, Nov 17 2015
Select[Range[240], !SquareFreeQ[2^#-1]&] (* Vladimir Joseph Stephan Orlovsky, Mar 18 2011 *)
default(factor_add_primes,1); is(n)=my(f=factor(n>>valuation(n,2))[,1],N,o); for(i=1,#f,if(n%(f[i]-1) == 0, return(1))); N=2^n-1; fordiv(n,d,f=factor(2^d-1)[,1]; for(i=1,#f, if(d==n,return(!issquarefree(N))); o=valuation(N,f[i]); if(o>1, return(1)); N/=f[i]^o)) \\ Charles R Greathouse IV, Feb 02 2014
is(n)=!issquarefree(2^n-1) \\ Charles R Greathouse IV, Feb 04 2014
a(4) = 2 because 2^4 - 1 = 15 = 3*5.
From _Gus Wiseman_, Jul 04 2019: (Start)
The sequence of Mersenne numbers together with their prime indices begins:
1: {}
3: {2}
7: {4}
15: {2,3}
31: {11}
63: {2,2,4}
127: {31}
255: {2,3,7}
511: {4,21}
1023: {2,5,11}
2047: {9,24}
4095: {2,2,3,4,6}
8191: {1028}
16383: {2,14,31}
32767: {4,11,36}
65535: {2,3,7,55}
131071: {12251}
262143: {2,2,2,4,8,21}
524287: {43390}
1048575: {2,3,3,5,11,13}
(End)
a[q_] := Module[{x, n}, x=FactorInteger[2^n-1]; n=Length[x]; Sum[Table[x[i][2], {i, n}][j], {j, n}]]
a[n_Integer] := PrimeOmega[2^n - 1]; Table[a[n], {n,200}] (* Vladimir Joseph Stephan Orlovsky, Jul 22 2011 *)
a(n)=bigomega(2^n-1) \\ Charles R Greathouse IV, Apr 01 2013
The fourth term is 23 (10111 in binary), since no prime less than 23 has exactly 4 1's in its binary representation.
a061712 n = fromJust $ find ((== n) . a000120) a000040_list -- Reinhard Zumkeller, Feb 10 2013
with(combstruct): a:=proc(n) local m,is,s,t,r; if n=1 then return 2 fi; r:=+infinity; for m from 0 to 100 do is := iterstructs(Combination(n-2+m),size=n-2); while not finished(is) do s := nextstruct(is); t := 2^(n-1+m)+1+add(2^i,i=s); # print(s,t); if isprime(t) then r:=min(t,r) fi; od; if r<+infinity then return r fi; od; return 0; end: seq(a(n),n=1..60); # Max Alekseyev, Aug 03 2005
Do[k = 1; While[ Count[ IntegerDigits[ Prime[k], 2], 1] != n, k++ ]; Print[ Prime[k]], {n, 1, 30} ]
(* Second program: *)
a[n_] := Module[{m, s, k, p}, For[m=0, True, m++, s = {1, Sequence @@ #, 1} & /@ Permutations[Join[Table[1, {n-2}], Table[0, {m}]]] // Sort; For[k=1, k <= Length[ s], k++, p = FromDigits[s[[k]], 2]; If[PrimeQ[p], Print["a(", n, ") = ", p]; Return[p]]]]]; a[1] = 2; Array[a, 100] (* Jean-François Alcover, Mar 16 2015 *)
Module[{hw=Table[{n,DigitCount[n,2,1]},{n,Prime[Range[250*10^6]]}]}, Table[ SelectFirst[hw,#[[2]]==k&],{k,31}]][[All,1]] (* Requires Mathematica version 10 or later *) (* Harvey P. Dale, Jan 01 2019 *)
a(n)=forprime(p=2, , if (hammingweight(p) == n, return(p));); \\ Michel Marcus, Mar 16 2015
from itertools import combinations from sympy import isprime def A061712(n): l, k = n-1, 2**n while True: for d in combinations(range(l-1,-1,-1),l-n+1): m = k-1 - sum(2**(e) for e in d) if isprime(m): return m l += 1 k *= 2 # Chai Wah Wu, Sep 02 2021
Element 2 = 4 because Mersenne2 = (2^3)-1 = 7; 7 is the 4th prime.
Array[PrimePi[2^MersennePrimeExponent[#] - 1] &, 8] (* Michael De Vlieger, Apr 21 2019 *)
LL(e) = if(e==2, return(1)); my(n, h); n = 2^e-1; h = Mod(2, n); for (k=1, e-2, h=2*h*h-1); return(0==h) \\ after Joerg Arndt in A000043 forprime(p=1, , if(LL(p), print1(primepi(2^p-1), ", "))) \\ Felix Fröhlich, Apr 19 2019
a(3) = 2 because 2^3 + 1 = 9 = 3*3.
a[n_] := Module[{x=FactorInteger[2^n+1]}, Sum[x[[i]][[2]], {i, Length[x]}]]
A054992[n_Integer] := PrimeOmega[2^n + 1]; Table[A054992[n], {n,103}] (* Vladimir Joseph Stephan Orlovsky, Jul 22 2011 *)
a(n)=bigomega(2^n+1) \\ Charles R Greathouse IV, Apr 29 2015
A137576[n_] := Module[{t}, (t = MultiplicativeOrder[2, 2 n + 1])* DivisorSum[2 n + 1, EulerPhi[#]/MultiplicativeOrder[2, #]&] - t + 1]; okQ[n_] := n > 1 && CompositeQ[n] && n == A137576[(n-1)/2]; Reap[For[k = 2, k < 2*10^6, k++, If[okQ[k], Print[k]; Sow[k]]]][[2, 1]] (* Jean-François Alcover, Jan 11 2019, from PARI *)
f(n)=my(t); sumdiv(2*n+1, d, eulerphi(d)/(t=znorder(Mod(2, d))))*t-t+1; \\ A137576 isok(n) = (n>1) && !isprime(n) && (n == f((n-1)/2)); \\ Michel Marcus, Oct 05 2018
Comments