cp's OEIS Frontend

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.

Showing 1-10 of 38 results. Next

A372402 Position of 2310^n among 11-smooth numbers A051038.

Original entry on oeis.org

1, 283, 3847, 20996, 74228, 203084, 469053, 960396, 1797086, 3135610, 5173909, 8156188, 12377846, 18190320, 26005929, 36302854, 49629820, 66611231, 87951744, 114441450, 146960432, 186483973, 234087084, 290949702, 358361266, 437725888, 530566933, 638532124, 763398291, 907076258
Offset: 0

Views

Author

Michael De Vlieger, Jun 03 2024

Keywords

Comments

Also position of 2310^(n+1) in A147572.

Crossrefs

Programs

  • Mathematica
    Table[
      Sum[Floor@ Log[11, 2310^n/(2^i*3^j*5^k*7^m)] + 1,
        {i, 0, Log[2, 2310^n]},
        {j, 0, Log[3, 2310^n/2^i]},
        {k, 0, Log[5, 2310^n/(2^i*3^j)]},
        {m, 0, Log[7, 2310^n/(2^i*3^j*5^k)]}],
      {n, 0, 8}]
  • Python
    # uses imports/function in A372401
    print(list(islice(A372401gen(p=11), 7))) # Michael S. Branicky, Jun 05 2024
    
  • Python
    from sympy import integer_log, prevprime
    def A372402(n):
        def g(x,m): return sum((x//3**i).bit_length() for i in range(integer_log(x,3)[0]+1)) if m==3 else sum(g(x//(m**i),prevprime(m))for i in range(integer_log(x,m)[0]+1))
        return g(2310**n,11) # Chai Wah Wu, Sep 16 2024

Extensions

a(14)-a(18) from Michael S. Branicky, Jun 05 2024
More terms from David A. Corneth, Jun 05 2024

A323051 Numbers that cannot be written as a sum of two or fewer 11-smooth numbers (A051038).

Original entry on oeis.org

479, 958, 1151, 1319, 1437, 1559, 1679, 1916, 2302, 2351, 2395, 2638, 2874, 2999, 3013, 3071, 3118, 3353, 3358, 3453, 3671, 3737, 3769, 3832, 3911, 3957, 4199, 4309, 4311, 4604, 4677, 4702, 4703, 4751, 4790, 4919, 5037, 5057, 5269, 5276, 5389, 5443, 5519, 5597, 5683
Offset: 1

Views

Author

Carlos Alves, Jan 03 2019

Keywords

Comments

Similar to A323046 (3-smooth), A323049 (5-smooth) or A323050 (7-smooth).
This sequence is a subsequence of A323046, A323049, and A323050.
Notice that A045535(4) = a(1) = 479.

Crossrefs

See A323046 (3-smooth), A323049 (5-smooth) or A323050 (7-smooth). Cf. A051038, A045535 (or A062241).

Programs

  • Mathematica
    f[n_] := Union@Flatten@Table[2^a*3^b*5^c*7^d, {a, 0, Log2[n]}, {b, 0, Log[3, n/2^a]}, {c, 0, Log[5, n/(2^a*3^b)]}, {d, 0, Log[7, n/(2^a*3^b*5^c)]}];
    b = Block[{nn = 3000, s}, s = f[nn]; {0, 1}~Join~
        Select[Union@Flatten@Outer[Plus, s, s], # <= nn &]];
    Complement[Range[3000], b]

A003586 3-smooth numbers: numbers of the form 2^i*3^j with i, j >= 0.

Original entry on oeis.org

1, 2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32, 36, 48, 54, 64, 72, 81, 96, 108, 128, 144, 162, 192, 216, 243, 256, 288, 324, 384, 432, 486, 512, 576, 648, 729, 768, 864, 972, 1024, 1152, 1296, 1458, 1536, 1728, 1944, 2048, 2187, 2304, 2592, 2916, 3072, 3456, 3888
Offset: 1

Views

Author

Paul Zimmermann, Dec 11 1996

Keywords

Comments

This sequence is easily confused with A033845, which gives numbers of the form 2^i*3^j with i, j >= 1. Don't simply say "numbers of the form 2^i*3^j", but specify which sequence you mean. - N. J. A. Sloane, May 26 2024
These numbers were once called "harmonic numbers", see Lenstra links. - N. J. A. Sloane, Jul 03 2015
Successive numbers k such that phi(6k) = 2k. - Artur Jasinski, Nov 05 2008
Where record values greater than 1 occur in A088468: A160519(n) = A088468(a(n)). - Reinhard Zumkeller, May 16 2009
Also numbers that are divisible by neither 6k - 1 nor 6k + 1, for all k > 0. - Robert G. Wilson v, Oct 26 2010
Also numbers m such that the rooted tree with Matula-Goebel number m has m antichains. The Matula-Goebel number of a rooted tree can be defined in the following recursive manner: to the one-vertex tree there corresponds the number 1; to a tree T with root degree 1 there corresponds the t-th prime number, where t is the Matula-Goebel number of the tree obtained from T by deleting the edge emanating from the root; to a tree T with root degree m>=2 there corresponds the product of the Matula-Goebel numbers of the m branches of T. The vertices of a rooted tree can be regarded as a partially ordered set, where u<=v holds for two vertices u and v if and only if u lies on the unique path between v and the root. An antichain is a nonempty set of mutually incomparable vertices. Example: m=4 is in the sequence because the corresponding rooted tree is \/=ARB (R is the root) having 4 antichains (A, R, B, AB). - Emeric Deutsch, Jan 30 2012
A204455(3*a(n)) = 3, and only for these numbers. - Wolfdieter Lang, Feb 04 2012
The number of terms less than or equal to n is Sum_{i=0..floor(log_2(n))} floor(log_3(n/2^i) + 1), or Sum_{i=0..floor(log_3(n))} floor(log_2(n/3^i) + 1), which requires fewer terms to compute. - Robert G. Wilson v, Aug 17 2012
Named 3-friables in French. - Michel Marcus, Jul 17 2013
In the 14th century Levi Ben Gerson proved that the only pairs of terms which differ by 1 are (1,2), (2,3), (3,4), and (8,9); see A235365, A235366, A236210. - Jonathan Sondow, Jan 20 2014
Range of values of A000005(n) (and also A181819(n)) for cubefree numbers n. - Matthew Vandermast, May 14 2014
A036561 is a permutation of this sequence. - L. Edson Jeffery, Sep 22 2014
Also the sorted union of A000244 and A007694. - Lei Zhou, Apr 19 2017
The sum of the reciprocals of the 3-smooth numbers is equal to 3. Brief proof: 1 + 1/2 + 1/3 + 1/4 + 1/6 + 1/8 + 1/9 + ... = (Sum_{k>=0} 1/2^k) * (Sum_{m>=0} 1/3^m) = (1/(1-1/2)) * (1/(1-1/3)) = (2/(2-1)) * (3/(3-1)) = 3. - Bernard Schott, Feb 19 2019
Also those integers k for which, for every prime p > 3, p^(2k) - 1 == 0 (mod 24k). - Federico Provvedi, May 23 2022
For n>1, the exponents’ parity {parity(i), parity(j)} of one out of four consecutive terms is {odd, odd}. Therefore, for n>1, at least one out of every four consecutive terms is a Zumkeller number (A083207). If for the term whose parity is {even, odd}, even also means nonzero, then this term is also a Zumkeller number (as is the case with the last of the four consecutive terms 1296, 1458, 1536, 1728). - Ivan N. Ianakiev, Jul 10 2022
Except the initial terms 2, 3, 4, 8, 9 and 16, these are numbers k such that k^6 divides 6^k. Except the initial terms 2, 3, 4, 6, 8, 9, 16, 18 and 27, these are numbers k such that k^12 divides 12^k. - Mohammed Yaseen, Jul 21 2022
In music theory, a comma is a ratio, close to 1 (typically less than 1.04), between two natural numbers divisible by only small primes (typically single digit). In this sequence, a(131) / a(130) = 531441 / 524288 ~ 1.013643 is the Pythagorean comma (A221363), the difference between 12 perfect fifths and 7 octaves. - Hal M. Switkay, Mar 23 2025

References

  • J.-M. De Koninck & A. Mercier, 1001 Problèmes en Théorie Classique des Nombres, Problème 654 pp. 85, 287-8, Ellipses Paris 2004.
  • S. Ramanujan, Collected Papers, Ed. G. H. Hardy et al., Cambridge 1927; Chelsea, NY, 1962, p. xxiv.
  • R. Tijdeman, Some applications of Diophantine approximation, pp. 261-284 of Surveys in Number Theory (Urbana, May 21, 2000), ed. M. A. Bennett et al., Peters, 2003.

Crossrefs

Cf. A051037, A002473, A051038, A080197, A080681, A080682, A117221, A105420, A062051, A117222, A117220, A090184, A131096, A131097, A186711, A186712, A186771, A088468, A061987, A080683 (p-smooth numbers with other values of p), A025613 (a subsequence).
Cf. also A000244, A007694. - Lei Zhou, Apr 19 2017
Cf. A191475 (successive values of i), A191476 (successive values of j), A022330 (indices of the pure terms 2^i), A022331 (indices of the pure terms 3^j). - N. J. A. Sloane, May 26 2024
Cf. A221363.

Programs

  • Haskell
    import Data.Set (Set, singleton, insert, deleteFindMin)
    smooth :: Set Integer -> [Integer]
    smooth s = x : smooth (insert (3*x) $ insert (2*x) s')
      where (x, s') = deleteFindMin s
    a003586_list = smooth (singleton 1)
    a003586 n = a003586_list !! (n-1)
    -- Reinhard Zumkeller, Dec 16 2010
    
  • Magma
    [n: n in [1..4000] | PrimeDivisors(n) subset [2,3]]; // Bruno Berselli, Sep 24 2012
  • Maple
    A003586 := proc(n) option remember; if n = 1 then 1; else for a from procname(n-1)+1 do numtheory[factorset](a) minus {2,3} ; if % = {} then return a; end if; end do: end if; end proc: # R. J. Mathar, Feb 28 2011
    with(numtheory): for i from 1 to 23328 do if(i/phi(i)=3)then print(i/6) fi od; # Gary Detlefs, Jun 28 2011
  • Mathematica
    a[1] = 1; j = 1; k = 1; n = 100; For[k = 2, k <= n, k++, If[2*a[k - j] < 3^j, a[k] = 2*a[k - j], {a[k] = 3^j, j++}]]; Table[a[i], {i, 1, n}] (* Hai He (hai(AT)mathteach.net) and Gilbert Traub, Dec 28 2004 *)
    aa = {}; Do[If[EulerPhi[6 n] == 2 n, AppendTo[aa, n]], {n, 1, 1000}]; aa (* Artur Jasinski, Nov 05 2008 *)
    fQ[n_] := Union[ MemberQ[{1, 5}, # ] & /@ Union@ Mod[ Rest@ Divisors@ n, 6]] == {False}; fQ[1] = True; Select[ Range@ 4000, fQ] (* Robert G. Wilson v, Oct 26 2010 *)
    powerOfTwo = 12; Select[Nest[Union@Join[#, 2*#, 3*#] &, {1}, powerOfTwo-1], # < 2^powerOfTwo &] (* Robert G. Wilson v and T. D. Noe, Mar 03 2011 *)
    fQ[n_] := n == 3 EulerPhi@ n; Select[6 Range@ 4000, fQ]/6 (* Robert G. Wilson v, Jul 08 2011 *)
    mx = 4000; Sort@ Flatten@ Table[2^i*3^j, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}] (* Robert G. Wilson v, Aug 17 2012 *)
    f[n_] := Block[{p2, p3 = 3^Range[0, Floor@ Log[3, n] + 1]}, p2 = 2^Floor[Log[2, n/p3] + 1]; Min[ Select[ p2*p3, IntegerQ]]]; NestList[f, 1, 54] (* Robert G. Wilson v, Aug 22 2012 *)
    Select[Range@4000, Last@Map[First, FactorInteger@#] <= 3 &] (* Vincenzo Librandi, Aug 25 2016 *)
    Select[Range[4000],Max[FactorInteger[#][[All,1]]]<4&] (* Harvey P. Dale, Jan 11 2017 *)
  • PARI
    test(n)=for(p=2,3, while(n%p==0, n/=p)); n==1;
    for(n=1,4000,if(test(n),print1(n",")))
    
  • PARI
    list(lim)=my(v=List(),N);for(n=0,log(lim\1+.5)\log(3),N=3^n;while(N<=lim,listput(v,N);N<<=1));vecsort(Vec(v)) \\ Charles R Greathouse IV, Jun 28 2011
    
  • PARI
    is_A003586(n)=n<5||vecmax(factor(n,5)[, 1])<5 \\ M. F. Hasler, Jan 16 2015
    
  • PARI
    list(lim)=my(v=List(), N); for(n=0, logint(lim\=1,3), N=3^n; while(N<=lim, listput(v, N); N<<=1)); Set(v) \\ Charles R Greathouse IV, Jan 10 2018
    
  • Python
    from itertools import count, takewhile
    def aupto(lim):
        pows2 = list(takewhile(lambda x: xMichael S. Branicky, Jul 08 2022
    
  • Python
    from sympy import integer_log
    def A003586(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//3**i).bit_length() for i in range(integer_log(x,3)[0]+1))
        return bisection(f,n,n) # Chai Wah Wu, Sep 15 2024
    
  • Python
    # faster for initial segment of sequence
    import heapq
    from itertools import islice
    def A003586gen(): # generator of terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], [2, 3]
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(A003586gen(), 65))) # Michael S. Branicky, Sep 17 2024
    (C++) // Returns A003586 <= threshold without approximations nor sorting
    #include 
    std::forward_list A003586(const int threshold) {
        std::forward_list sequence;
        auto start_it = sequence.before_begin();
        for (int i = 1; i <= threshold; i *= 2) {
            for (int inc = 1; std::next(start_it) != sequence.end() && inc <= i; inc *= 3)
                ++start_it;
            auto it = start_it;
            for (int j = 1; i * j <= threshold; j *= 3) {
                sequence.emplace_after(it, i * j);
                for (int inc = 1; std::next(it) != sequence.end() && inc <= i; inc *= 2)
                    ++it;
            }
        }
        return sequence;
    } // Eben Gino Lester, Apr 17 2025
    
  • Sage
    def isA003586(n) :
        return not any(d != 2 and d != 3 for d in prime_divisors(n))
    @CachedFunction
    def A003586(n) :
        if n == 1 : return 1
        k = A003586(n-1) + 1
        while not isA003586(k) : k += 1
        return k
    [A003586(n) for n in (1..55)] # Peter Luschny, Jul 20 2012
    

Formula

An asymptotic formula for a(n) is roughly a(n) ~ 1/sqrt(6)*exp(sqrt(2*log(2)*log(3)*n)). - Benoit Cloitre, Nov 20 2001
A061987(n) = a(n + 1) - a(n), a(A084791(n)) = A084789(n), a(A084791(n) + 1) = A084790(n). - Reinhard Zumkeller, Jun 03 2003
Union of powers of 2 and 3 with n such that psi(n) = 2*n, where psi(n) = n*Product_(1 + 1/p) over all prime factors p of n = A001615(n). - Lekraj Beedassy, Sep 07 2004; corrected by Franklin T. Adams-Watters, Mar 19 2009
a(n) = 2^A022328(n)*3^A022329(n). - N. J. A. Sloane, Mar 19 2009
The characteristic function of this sequence is given by Sum_{n >= 1} x^a(n) = Sum_{n >= 1} moebius(6*n)*x^n/(1 - x^n). - Paul D. Hanna, Sep 18 2011
a(n) = A007694(n+1)/2. - Lei Zhou, Apr 19 2017

Extensions

Deleted claim that this sequence is union of 2^n (A000079) and 3^n (A000244) sequences -- this does not include the terms which are not pure powers. - Walter Roscello (wroscello(AT)comcast.net), Nov 16 2008

A002473 7-smooth numbers: positive numbers whose prime divisors are all <= 7.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, 48, 49, 50, 54, 56, 60, 63, 64, 70, 72, 75, 80, 81, 84, 90, 96, 98, 100, 105, 108, 112, 120, 125, 126, 128, 135, 140, 144, 147, 150, 160, 162, 168, 175, 180, 189, 192
Offset: 1

Views

Author

Keywords

Comments

Also called humble numbers; sometimes also called highly composite numbers, but this usually refers to A002182.
Successive numbers k such that phi(210k) = 48k. - Artur Jasinski, Nov 05 2008
The divisors of 10! (A161466) are a finite subsequence. - Reinhard Zumkeller, Jun 10 2009
Numbers n such that A198487(n) > 0 and A107698(n) > 0. - Jaroslav Krizek, Nov 04 2011
A262401(a(n)) = a(n). - Reinhard Zumkeller, Sep 25 2015
Numbers which are products of single-digit numbers. - N. J. A. Sloane, Jul 02 2017
Phi(a(n)) is 7-smooth. In fact, the Euler Phi function applied to p-smooth numbers, for any prime p, is p-smooth. - Richard Locke Peterson, May 09 2020
Also those integers k, such that, for every prime p > 5, p^(12k) - 1 == 0 (mod 5040k). - Federico Provvedi, Jun 06 2022
The nonprimes with this property are all terms except for 2, 3, 5 and 7, i.e.: (1, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, ...); the composite terms are all but the first one of this subsequence. ["Trivial" data provided mainly for search purpose.] - M. F. Hasler, Jun 06 2023

References

  • B. C. Berndt, Ramanujan's Notebooks Part IV, Springer-Verlag, see p. 52.
  • N. J. A. Sloane, A Handbook of Integer Sequences, Academic Press, 1973 (includes this sequence).
  • N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).

Crossrefs

Subsequence of A080672, complement of A068191. Subsequences: A003591, A003594, A003595, A195238, A059405.
Not the same as A063938. For p-smooth numbers with other values of p, see A003586, A051037, A051038, A080197, A080681, A080682, A080683.
Cf. A002182, A067374, A210679, A238985 (zeroless terms), A006530.
Cf. A262401.

Programs

  • Haskell
    import Data.Set (singleton, deleteFindMin, fromList, union)
    a002473 n = a002473_list !! (n-1)
    a002473_list = f $ singleton 1 where
       f s = x : f (s' `union` fromList (map (* x) [2,3,5,7]))
             where (x, s') = deleteFindMin s
    -- Reinhard Zumkeller, Mar 08 2014, Apr 02 2012, Apr 01 2012
    
  • Magma
    [n: n in [1..200] | PrimeDivisors(n) subset PrimesUpTo(7)]; // Bruno Berselli, Sep 24 2012
    
  • Mathematica
    Select[Range[250], Max[Transpose[FactorInteger[ # ]][[1]]]<=7&]
    aa = {}; Do[If[EulerPhi[210 n] == 48 n, AppendTo[aa, n]], {n, 1, 1200}]; aa (* Artur Jasinski, Nov 05 2008 *)
    mxExp = 8; Select[Union[Times @@@ Flatten[Table[Tuples[{2, 3, 5, 7}, n], {n, mxExp}], 1]], # <= 2^mxExp &] (* Harvey P. Dale, Aug 13 2012 *)
    mx = 200; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}] (* Robert G. Wilson v, Aug 17 2012 *)
  • PARI
    test(n)=m=n; forprime(p=2,7, while(m%p==0,m=m/p)); return(m==1)
    for(n=1,200,if(test(n),print1(n",")))
    
  • PARI
    is_A002473(n)=n<11||vecmax(factor(n,8)[,1])<8 \\ M. F. Hasler, Jan 16 2015
    
  • PARI
    list(lim)=my(v=List(),t); for(a=0,logint(lim\1,7), for(b=0,logint(lim\7^a,5), for(c=0,logint(lim\7^a\5^b,3), t=3^c*5^b*7^a; while(t<=lim, listput(v,t); t<<=1)))); Set(v) \\ Charles R Greathouse IV, Feb 22 2017
    
  • Python
    import heapq
    from itertools import islice
    from sympy import primerange
    def A002473gen(p=7): # generate all p-smooth terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(A002473gen(), 65))) # Michael S. Branicky, Nov 19 2022
    
  • Python
    from sympy import integer_log
    def A002473(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):
            c = n+x
            for i in range(integer_log(x,7)[0]+1):
                i7 = 7**i
                m = x//i7
                for j in range(integer_log(m,5)[0]+1):
                    j5 = 5**j
                    r = m//j5
                    for k in range(integer_log(r,3)[0]+1):
                        c -= (r//3**k).bit_length()
            return c
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024

Formula

A006530(a(n)) <= 7. - Reinhard Zumkeller, Apr 01 2012
Sum_{n>=1} 1/a(n) = Product_{primes p <= 7} p/(p-1) = (2*3*5*7)/(1*2*4*6) = 35/8. - Amiram Eldar, Sep 22 2020

Extensions

More terms from James Sellers, Dec 23 1999
Additional comments from Michel Lecomte, Jun 09 2007
Edited by M. F. Hasler, Jan 16 2015

A051037 5-smooth numbers, i.e., numbers whose prime divisors are all <= 5.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405
Offset: 1

Views

Author

Keywords

Comments

Sometimes called the Hamming sequence, since Hamming asked for an efficient algorithm to generate the list, in ascending order, of all numbers of the form 2^i*3^j*5^k for i,j,k >= 0. The problem was popularized by Edsger Dijkstra.
Numbers k such that 8*k = EulerPhi(30*k). - Artur Jasinski, Nov 05 2008
Where record values greater than 1 occur in A165704: A165705(n) = A165704(a(n)). - Reinhard Zumkeller, Sep 26 2009
Also called "harmonic whole numbers", see Howard and Longair, 1982, Table I, page 121. - Hugo Pfoertner, Jul 16 2020
Also called ugly numbers, although it is not clear why. - Gus Wiseman, May 21 2021
Some woody bamboo species have extraordinarily long and stable flowering intervals that belong to this sequence. The model by Veller, Nowak & Davis justifies this observation from the evolutionary point of view. - Andrey Zabolotskiy, Jun 27 2021
Also those integers k for which, for every prime p > 5, p^(4*k) - 1 == 0 (mod 240*k). - Federico Provvedi, May 23 2022
As noted in the comments to A085152, Størmer's theorem implies that the only pairs of consecutive integers that appear as consecutive terms of this sequence are (1,2), (2,3), (3,4), (4,5), (5,6), (8,9), (9,10), (15,16), (24,25), and (80,81). These all represent significant musical intervals. - Hal M. Switkay, Dec 05 2022

Examples

			From _Gus Wiseman_, May 21 2021: (Start)
The sequence of terms together with their prime indices begins:
      1: {}            25: {3,3}
      2: {1}           27: {2,2,2}
      3: {2}           30: {1,2,3}
      4: {1,1}         32: {1,1,1,1,1}
      5: {3}           36: {1,1,2,2}
      6: {1,2}         40: {1,1,1,3}
      8: {1,1,1}       45: {2,2,3}
      9: {2,2}         48: {1,1,1,1,2}
     10: {1,3}         50: {1,3,3}
     12: {1,1,2}       54: {1,2,2,2}
     15: {2,3}         60: {1,1,2,3}
     16: {1,1,1,1}     64: {1,1,1,1,1,1}
     18: {1,2,2}       72: {1,1,1,2,2}
     20: {1,1,3}       75: {2,3,3}
     24: {1,1,1,2}     80: {1,1,1,1,3}
(End)
		

Crossrefs

Subsequences: A003592, A003593, A051916 , A257997.
For p-smooth numbers with other values of p, see A003586, A002473, A051038, A080197, A080681, A080682, A080683.
The partitions with these Heinz numbers are counted by A001399.
The conjugate opposite is A033942, counted by A004250.
The opposite is A059485, counted by A004250.
The non-3-smooth case is A080193, counted by A069905.
The conjugate is A037144, counted by A001399.
The complement is A279622, counted by A035300.
Requiring the sum of prime indices to be even gives A344297.

Programs

  • Haskell
    import Data.Set (singleton, deleteFindMin, insert)
    a051037 n = a051037_list !! (n-1)
    a051037_list = f $ singleton 1 where
       f s = y : f (insert (5 * y) $ insert (3 * y) $ insert (2 * y) s')
                   where (y, s') = deleteFindMin s
    -- Reinhard Zumkeller, May 16 2015
    
  • Magma
    [n: n in [1..500] | PrimeDivisors(n) subset [2,3,5]]; // Bruno Berselli, Sep 24 2012
    
  • Maple
    A051037 := proc(n)
        option remember;
        local a;
        if n = 1 then
            1;
        else
            for a from procname(n-1)+1 do
                numtheory[factorset](a) minus {2, 3,5 } ;
                if % = {} then
                    return a;
                end if;
            end do:
        end if;
    end proc:
    seq(A051037(n),n=1..100) ; # R. J. Mathar, Nov 05 2017
  • Mathematica
    mx = 405; Sort@ Flatten@ Table[ 2^a*3^b*5^c, {a, 0, Log[2, mx]}, {b, 0, Log[3, mx/2^a]}, {c, 0, Log[5, mx/(2^a*3^b)]}] (* Or *)
    Select[ Range@ 405, Last@ Map[First, FactorInteger@ #] < 7 &] (* Robert G. Wilson v *)
    With[{nn=10},Select[Union[Times@@@Flatten[Table[Tuples[{2,3,5},n],{n,0,nn}],1]],#<=2^nn&]] (* Harvey P. Dale, Feb 28 2022 *)
  • PARI
    test(n)= {m=n; forprime(p=2,5, while(m%p==0,m=m/p)); return(m==1)}
    for(n=1,500,if(test(n),print1(n",")))
    
  • PARI
    a(n)=local(m); if(n<1,0,n=a(n-1); until(if(m=n, forprime(p=2,5, while(m%p==0,m/=p)); m==1),n++); n)
    
  • PARI
    list(lim)=my(v=List(),s,t); for(i=0,logint(lim\=1,5), t=5^i; for(j=0,logint(lim\t,3), s=t*3^j; while(s<=lim, listput(v,s); s<<=1))); Set(v) \\ Charles R Greathouse IV, Sep 21 2011; updated Sep 19 2016
    
  • PARI
    smooth(P:vec,lim)={ my(v=List([1]),nxt=vector(#P,i,1),indx,t);
    while(1, t=vecmin(vector(#P,i,v[nxt[i]]*P[i]),&indx);
    if(t>lim,break); if(t>v[#v],listput(v,t)); nxt[indx]++);
    Vec(v)
    };
    smooth([2,3,5], 1e4) \\ Charles R Greathouse IV, Dec 03 2013
    
  • PARI
    is_A051037(n)=n<7||vecmax(factor(n,6)[, 1])<7 \\ M. F. Hasler, Jan 16 2015
    
  • Python
    def isok(n):
      while n & 1 == 0: n >>= 1
      while n % 3 == 0: n //= 3
      while n % 5 == 0: n //= 5
      return n == 1 #  Darío Clavijo, Dec 30 2022
    
  • Python
    from sympy import integer_log
    def A051037(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):
            c = n+x
            for i in range(integer_log(x,5)[0]+1):
                for j in range(integer_log(y:=x//5**i,3)[0]+1):
                    c -= (y//3**j).bit_length()
            return c
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024
    
  • Python
    # faster for initial segment of sequence
    import heapq
    from itertools import islice
    def A051037gen(): # generator of terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], [2, 3, 5]
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                        heapq.heappush(h, v*p)
    print(list(islice(A051037gen(), 65))) # Michael S. Branicky, Sep 17 2024

Formula

Let s(n) = Card(k | a(k)Benoit Cloitre, Dec 30 2001
The characteristic function of this sequence is given by:
Sum_{n>=1} x^a(n) = Sum_{n>=1} -Möbius(30*n)*x^n/(1-x^n). - Paul D. Hanna, Sep 18 2011
a(n) = A143207(n) / 30. - Reinhard Zumkeller, Sep 13 2011
A204455(15*a(n)) = 15, and only for these numbers. - Wolfdieter Lang, Feb 04 2012
A006530(a(n)) <= 5. - Reinhard Zumkeller, May 16 2015
Sum_{n>=1} 1/a(n) = Product_{primes p <= 5} p/(p-1) = (2*3*5)/(1*2*4) = 15/4. - Amiram Eldar, Sep 22 2020

A080197 13-smooth numbers: numbers whose prime divisors are all <= 13.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48, 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81, 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120
Offset: 1

Views

Author

Klaus Brockhaus, Feb 10 2003

Keywords

Comments

Numbers of the form 2^r*3^s*5^t*7^u*11^v*13^w with r, s, t, u, v, w >= 0.

Examples

			33 = 3*11 and 39 = 3*13 are terms but 34 = 2*17 is not.
		

Crossrefs

Cf. A000079, A080196. For p-smooth numbers with other values of p, see A003586, A051037, A002473, A051038, A080681, A080682, A080683.

Programs

  • Magma
    [n: n in [1..150] | PrimeDivisors(n) subset PrimesUpTo(13)]; // Bruno Berselli, Sep 24 2012
    
  • Mathematica
    mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}] (* Robert G. Wilson v, Aug 17 2012 *)
  • PARI
    test(n)=m=n; forprime(p=2,13, while(m%p==0,m=m/p)); return(m==1)
    for(n=1,200,if(test(n),print1(n",")))
    
  • PARI
    is_A080197(n,p=13)=n<=p||vecmax(factor(n,p+1)[,1])<=p \\ M. F. Hasler, Jan 16 2015
    
  • PARI
    list(lim,p=13)=if(p==2, return(powers(2, logint(lim\1,2)))); my(v=[],q=precprime(p-1),t=1); for(e=0,logint(lim\=1,p), v=concat(v, list(lim\t,q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020
    
  • Python
    import heapq
    from itertools import islice
    from sympy import primerange
    def agen(p=13): # generate all p-smooth terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(agen(), 69))) # Michael S. Branicky, Nov 20 2022
    
  • Python
    from sympy import integer_log, prevprime
    def A080197(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 g(x,m): return sum((x//3**i).bit_length() for i in range(integer_log(x,3)[0]+1)) if m==3 else sum(g(x//(m**i),prevprime(m))for i in range(integer_log(x,m)[0]+1))
        def f(x): return n+x-g(x,13)
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024

Formula

Sum_{n>=1} 1/a(n) = Product_{primes p <= 13} p/(p-1) = (2*3*5*7*11*13)/(1*2*4*6*10*12) = 1001/192. - Amiram Eldar, Sep 22 2020

A033620 Numbers all of whose prime factors are palindromes.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 27, 28, 30, 32, 33, 35, 36, 40, 42, 44, 45, 48, 49, 50, 54, 55, 56, 60, 63, 64, 66, 70, 72, 75, 77, 80, 81, 84, 88, 90, 96, 98, 99, 100, 101, 105, 108, 110, 112, 120, 121, 125, 126, 128, 131
Offset: 1

Views

Author

N. J. A. Sloane, May 17 1998

Keywords

Comments

Multiplicative closure of A002385; A051038 and A046368 are subsequences. - Reinhard Zumkeller, Apr 11 2011

Examples

			10 = 2 * 5 is a term since both 2 and 5 are palindromes.
110 = 2 * 5 * 11 is a term since 2, 5 and 11 are palindromes.
		

Crossrefs

Programs

  • Haskell
    a033620 n = a033620_list !! (n-1)
    a033620_list = filter chi [1..] where
       chi n = a136522 spf == 1 && (n' == 1 || chi n') where
          n' = n `div` spf
          spf = a020639 n  -- cf. A020639
    -- Reinhard Zumkeller, Apr 11 2011
    
  • Maple
    N:= 5: # to get all terms of up to N digits
    digrev:= proc(t) local L; L:= convert(t,base,10);
    add(L[-i-1]*10^i,i=0..nops(L)-1);
    end proc:
    PPrimes:= [2,3,5,7,11]:
    for d from 3 to N by 2 do
        m:= (d-1)/2;
        PPrimes:= PPrimes, select(isprime,[seq(seq(n*10^(m+1)+y*10^m+digrev(n), y=0..9), n=10^(m-1)..10^m-1)]);
    od:
    PPrimes:= map(op,[PPrimes]):
    M:= 10^N:
    B:= Vector(M);
    B[1]:= 1:
    for p in PPrimes do
      for k from 1 to floor(log[p](M)) do
         R:= [$1..floor(M/p^k)];
         B[p^k*R] := B[p^k*R] + B[R]
      od
    od:
    select(t -> B[t] > 0, [$1..M]); # Robert Israel, Jul 05 2015
    # alternative
    isA033620:= proc(n)
        for d in numtheory[factorset](n) do
            if not isA002113(op(1,d)) then
                return false;
            end if;
        end do;
        true ;
    end proc:
    for n from 1 to 300 do
        if isA033620(n) then
            printf("%d,",n) ;
        end if;
    end do: # R. J. Mathar, Sep 09 2015
  • Mathematica
    palQ[n_]:=Reverse[x=IntegerDigits[n]]==x; Select[Range[131],And@@palQ/@First/@FactorInteger[#]&] (* Jayanta Basu, Jun 05 2013 *)
  • PARI
    ispal(n)=n=digits(n);for(i=1,#n\2,if(n[i]!=n[#n+1-i],return(0)));1
    is(n)=if(n<13,n>0,vecmin(apply(ispal,factor(n)[,1]))) \\ Charles R Greathouse IV, Feb 06 2013
    
  • Python
    from sympy import isprime, primefactors
    def pal(n): s = str(n); return s == s[::-1]
    def ok(n): return all(pal(f) for f in primefactors(n))
    print(list(filter(ok, range(1, 132)))) # Michael S. Branicky, Apr 06 2021

Formula

Sum_{n>=1} 1/a(n) = Product_{p in A002385} p/(p-1) = 5.0949... - Amiram Eldar, Sep 27 2020

A080682 19-smooth numbers: numbers whose prime divisors are all <= 19.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 38, 39, 40, 42, 44, 45, 48, 49, 50, 51, 52, 54, 55, 56, 57, 60, 63, 64, 65, 66, 68, 70, 72, 75, 76, 77, 78, 80, 81, 84, 85, 88, 90, 91, 95, 96, 98, 99, 100
Offset: 1

Views

Author

Cino Hilliard, Mar 02 2003

Keywords

Crossrefs

For p-smooth numbers with other values of p, see A003586, A051037, A002473, A051038, A080197, A080681, A080683.

Programs

  • Magma
    [n: n in [1..100] | PrimeDivisors(n) subset PrimesUpTo(19)]; // Bruno Berselli, Sep 24 2012
    
  • Mathematica
    mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n*17^o*19^p, {i, 0, Log[2,mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]},{l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}, {o, 0, Log[17, mx/(2^i*3^j*5^k*7^l*11^m*13^n)]}, {p, 0, Log[19, mx/(2^i*3^j*5^k*7^l*11^m*13^n*17^o)]}] (* Robert G. Wilson v, Jan 19 2016 *)
    Select[Range[100],Max[FactorInteger[#][[All,1]]]<20&] (* Harvey P. Dale, Sep 20 2018 *)
  • PARI
    test(n)= {m=n; forprime(p=2,19, while(m%p==0,m=m/p)); return(m==1)}
    for(n=1,200,if(test(n),print1(n",")))
    
  • PARI
    list(lim,p=19)=if(p==2, return(powers(2, logint(lim\1,2)))); my(v=[],q=precprime(p-1),t=1); for(e=0,logint(lim\=1,p), v=concat(v, list(lim\t,q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020
    
  • Python
    import heapq
    from itertools import islice
    from sympy import primerange
    def agen(p=19): # generate all p-smooth terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(agen(), 72))) # Michael S. Branicky, Nov 20 2022
    
  • Python
    from sympy import integer_log
    def A080682(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 g(x,m): return sum((x//3**i).bit_length() for i in range(integer_log(x,3)[0]+1)) if m==3 else sum(g(x//(m**i),prevprime(m))for i in range(integer_log(x,m)[0]+1))
        def f(x): return n+x-g(x,19)
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024

Formula

Sum_{n>=1} 1/a(n) = Product_{primes p <= 19} p/(p-1) = (2*3*5*7*11*13*17*19)/(1*2*4*6*10*12*16*18) = 323323/55296. - Amiram Eldar, Sep 22 2020

A080683 23-smooth numbers: numbers whose prime divisors are all <= 23.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 48, 49, 50, 51, 52, 54, 55, 56, 57, 60, 63, 64, 65, 66, 68, 69, 70, 72, 75, 76, 77, 78, 80, 81, 84, 85, 88, 90, 91, 92, 95
Offset: 1

Views

Author

Cino Hilliard, Mar 02 2003

Keywords

Comments

Coincides for the first 111 terms with A174228 (divisors of 24!). - Bruno Berselli, Sep 24 2012

Crossrefs

For p-smooth numbers with other values of p, see A003586, A051037, A002473, A051038, A080197, A080681, A080682.

Programs

  • Magma
    [n: n in [1..100] | PrimeDivisors(n) subset PrimesUpTo(23)]; // Bruno Berselli, Sep 24 2012
    
  • Maple
    select(t -> max(numtheory:-factorset(t)) <= 23, [$1..1000]); # Robert Israel, Jan 22 2016
  • Mathematica
    mx = 100; Sort@ Flatten@ Table[ 2^a*3^b*5^c*7^d*11^e*13^f*17^g*19^h*23^i, {a, 0, Log[2, mx]}, {b, 0, Log[3, mx/2^a]}, {c, 0, Log[5, mx/(2^a*3^b)]}, {d, 0, Log[7, mx/(2^a*3^b*5^c)]}, {e, 0, Log[11, mx/(2^a*3^b*5^c*7^d)]}, {f, 0, Log[13, mx/(2^a*3^b*5^c*7^d*11^e)]}, {g, 0, Log[17, mx/(2^a*3^b*5^c*7^d*11^e*13^f)]}, {h, 0, Log[19, mx/(2^a*3^b*5^c*7^d*11^e*13^f*17^g)]}, {i, 0, Log[23, mx/(2^a*3^b*5^c*7^d*11^e*13^f*17^g*19^h)]}] (* Robert G. Wilson v, Jan 19 2016 *)
  • PARI
    test(n)=m=n; forprime(p=2,23, while(m%p==0,m=m/p)); return(m==1)
    for(n=1,100,if(test(n),print1(n",")))
    
  • PARI
    list(lim,p=23)=if(p==2, return(powers(2, logint(lim\1,2)))); my(v=[],q=precprime(p-1),t=1); for(e=0,logint(lim\=1,p), v=concat(v, list(lim\t,q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020
    
  • Python
    import heapq
    from itertools import islice
    from sympy import primerange
    def agen(p=23): # generate all p-smooth terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(agen(), 72))) # Michael S. Branicky, Nov 20 2022
    
  • Python
    from sympy import integer_log, prevprime
    def A080683(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 g(x,m): return sum((x//3**i).bit_length() for i in range(integer_log(x,3)[0]+1)) if m==3 else sum(g(x//(m**i),prevprime(m))for i in range(integer_log(x,m)[0]+1))
        def f(x): return n+x-g(x,23)
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024

Formula

Sum_{n>=1} 1/a(n) = Product_{primes p <= 23} p/(p-1) = (2*3*5*7*11*13*17*19*23)/(1*2*4*6*10*12*16*18*22) = 676039/110592. - Amiram Eldar, Sep 22 2020

A080681 17-smooth numbers: numbers whose prime divisors are all <= 17.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 39, 40, 42, 44, 45, 48, 49, 50, 51, 52, 54, 55, 56, 60, 63, 64, 65, 66, 68, 70, 72, 75, 77, 78, 80, 81, 84, 85, 88, 90, 91, 96, 98, 99, 100, 102, 104, 105
Offset: 1

Views

Author

Cino Hilliard, Mar 02 2003

Keywords

Crossrefs

For p-smooth numbers with other values of p, see A003586, A051037, A002473, A051038, A080197, A080682, A080683.

Programs

  • Magma
    [n: n in [1..150] | PrimeDivisors(n) subset PrimesUpTo(17)]; // Bruno Berselli, Sep 24 2012
    
  • Mathematica
    mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n*17^o, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}, {o, 0, Log[17, mx/(2^i*3^j*5^k*7^l*11^m*13^n)]}] (* Robert G. Wilson v, Aug 17 2012 *)
  • PARI
    test(n)= {m=n; forprime(p=2,17, while(m%p==0,m=m/p)); return(m==1)}
    for(n=1,200,if(test(n),print1(n",")))
    
  • PARI
    list(lim,p=17)=if(p==2, return(powers(2, logint(lim\1,2)))); my(v=[],q=precprime(p-1),t=1); for(e=0,logint(lim\=1,p), v=concat(v, list(lim\t,q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020
    
  • Python
    import heapq
    from itertools import islice
    from sympy import primerange
    def agen(p=17): # generate all p-smooth terms
        v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
        while True:
            v = heapq.heappop(h)
            if v != oldv:
                yield v
                oldv = v
                for p in psmooth_primes:
                    heapq.heappush(h, v*p)
    print(list(islice(agen(), 70))) # Michael S. Branicky, Nov 20 2022
    
  • Python
    from sympy import integer_log, prevprime
    def A080681(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 g(x,m): return sum((x//3**i).bit_length() for i in range(integer_log(x,3)[0]+1)) if m==3 else sum(g(x//(m**i),prevprime(m))for i in range(integer_log(x,m)[0]+1))
        def f(x): return n+x-g(x,17)
        return bisection(f,n,n) # Chai Wah Wu, Sep 16 2024

Formula

Sum_{n>=1} 1/a(n) = Product_{primes p <= 17} p/(p-1) = (2*3*5*7*11*13*17)/(1*2*4*6*10*12*16) = 17017/3072. - Amiram Eldar, Sep 22 2020
Showing 1-10 of 38 results. Next