OFFSET
1,2
COMMENTS
A204455(7*a(n)) = 7, and only for these numbers. - Wolfdieter Lang, Feb 04 2012
LINKS
Reinhard Zumkeller, Table of n, a(n) for n = 1..10000 (first 100 terms from Vincenzo Librandi)
Vaclav Kotesovec, Graph - the asymptotic ratio (250000 terms)
FORMULA
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
Sum_{n>=1} 1/a(n) = (2*7)/((2-1)*(7-1)) = 7/3. - Amiram Eldar, Sep 22 2020
a(n) ~ exp(sqrt(2*log(2)*log(7)*n)) / sqrt(14). - Vaclav Kotesovec, Sep 22 2020
MATHEMATICA
fQ[n_] := PowerMod[14, n, n]==0; Select[Range[30000], fQ] (* Vincenzo Librandi, Feb 04 2012 *)
PROG
(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
(PARI) isA003591(n)=n>>=valuation(n, 2); ispower(n, , &n); n==1||n==7 \\ Charles R Greathouse IV, Jun 28 2011
(Magma) [n: n in [1..26000] | PrimeDivisors(n) subset [2, 7]]; // Bruno Berselli, Sep 24 2012
(Haskell)
import Data.Set (singleton, deleteFindMin, insert)
a003591 n = a003591_list !! (n-1)
a003591_list = f $ singleton 1 where
f s = y : f (insert (2 * y) $ insert (7 * y) s')
where (y, s') = deleteFindMin s
-- Reinhard Zumkeller, May 16 2015
(GAP) Filtered([1..30000], n->PowerMod(14, n, n)=0); # Muniru A Asiru, Mar 19 2019
(Python)
from sympy import integer_log
def A003591(n):
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x): return n+x-sum((x//7**i).bit_length() for i in range(integer_log(x, 7)[0]+1))
return bisection(f, n, n) # Chai Wah Wu, Sep 16 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
STATUS
approved