%I #56 Sep 18 2024 03:44:13
%S 1,2,4,7,8,14,16,28,32,49,56,64,98,112,128,196,224,256,343,392,448,
%T 512,686,784,896,1024,1372,1568,1792,2048,2401,2744,3136,3584,4096,
%U 4802,5488,6272,7168,8192,9604,10976,12544,14336,16384,16807,19208,21952,25088
%N Numbers of form 2^i*7^j, with i, j >= 0.
%C A204455(7*a(n)) = 7, and only for these numbers. - _Wolfdieter Lang_, Feb 04 2012
%H Reinhard Zumkeller, <a href="/A003591/b003591.txt">Table of n, a(n) for n = 1..10000</a> (first 100 terms from Vincenzo Librandi)
%H Vaclav Kotesovec, <a href="/A003591/a003591.jpg">Graph - the asymptotic ratio (250000 terms)</a>
%F The characteristic function of this sequence is given by Sum_{n >= 1} x^a(n) = Sum_{n >= 1} mu(14*n)*x^n/(1 - x^n), where mu(n) is the Möbius function A008683. Cf. with the formula of Hanna in A051037. - _Peter Bala_, Mar 18 2019
%F Sum_{n>=1} 1/a(n) = (2*7)/((2-1)*(7-1)) = 7/3. - _Amiram Eldar_, Sep 22 2020
%F a(n) ~ exp(sqrt(2*log(2)*log(7)*n)) / sqrt(14). - _Vaclav Kotesovec_, Sep 22 2020
%t fQ[n_] := PowerMod[14,n,n]==0; Select[Range[30000], fQ] (* _Vincenzo Librandi_, Feb 04 2012 *)
%o (PARI) list(lim)=my(v=List(),N);for(n=0,log(lim)\log(7),N=7^n;while(N<=lim,listput(v,N);N<<=1));vecsort(Vec(v)) \\ _Charles R Greathouse IV_, Jun 28 2011
%o (PARI) isA003591(n)=n>>=valuation(n,2);ispower(n,,&n);n==1||n==7 \\ _Charles R Greathouse IV_, Jun 28 2011
%o (Magma) [n: n in [1..26000] | PrimeDivisors(n) subset [2,7]]; // _Bruno Berselli_, Sep 24 2012
%o (Haskell)
%o import Data.Set (singleton, deleteFindMin, insert)
%o a003591 n = a003591_list !! (n-1)
%o a003591_list = f $ singleton 1 where
%o f s = y : f (insert (2 * y) $ insert (7 * y) s')
%o where (y, s') = deleteFindMin s
%o -- _Reinhard Zumkeller_, May 16 2015
%o (GAP) Filtered([1..30000],n->PowerMod(14,n,n)=0); # _Muniru A Asiru_, Mar 19 2019
%o (Python)
%o from sympy import integer_log
%o def A003591(n):
%o def bisection(f,kmin=0,kmax=1):
%o while f(kmax) > kmax: kmax <<= 1
%o while kmax-kmin > 1:
%o kmid = kmax+kmin>>1
%o if f(kmid) <= kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o return kmax
%o def f(x): return n+x-sum((x//7**i).bit_length() for i in range(integer_log(x,7)[0]+1))
%o return bisection(f,n,n) # _Chai Wah Wu_, Sep 16 2024
%Y Cf. A003586, A003592, A003593, A003594, A003595.
%K nonn
%O 1,2
%A _N. J. A. Sloane_