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.

Previous Showing 11-20 of 47 results. Next

A089233 Number of coprime pairs of divisors > 1 of n.

Original entry on oeis.org

0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 1, 1, 0, 0, 2, 0, 2, 1, 1, 0, 3, 0, 1, 0, 2, 0, 6, 0, 0, 1, 1, 1, 4, 0, 1, 1, 3, 0, 6, 0, 2, 2, 1, 0, 4, 0, 2, 1, 2, 0, 3, 1, 3, 1, 1, 0, 11, 0, 1, 2, 0, 1, 6, 0, 2, 1, 6, 0, 6, 0, 1, 2, 2, 1, 6, 0, 4, 0, 1, 0, 11, 1, 1, 1, 3, 0, 11, 1, 2, 1, 1, 1, 5, 0, 2, 2, 4, 0, 6, 0, 3, 6
Offset: 1

Views

Author

Reinhard Zumkeller, Dec 11 2003

Keywords

Comments

Also the number of divisors of n^2 which do not divide n and which are less than n. See link for proof. - Andrew Weimholt, Dec 06 2009
a(A000961(n)) = 0; a(A006881(n)) = 1; a(A054753(n)) = 2; a(A065036(n)) = 3. - Robert G. Wilson v, Dec 16 2009
First occurrence of k beginning with 0: 1, 6, 12, 24, 36, 96, 30, 384, 144, 216, 288, 60, 432, 24576, 1152, 864, 120, 393216, 1728, 1572864, 180, 240, 18432, 25165824, 5184, 210, 480, 13824, 10368, 360, 15552, 960, 20736, 55296, 1179648, 31104, 900, ..., . Except for 1, each is divisible by 6. Also the first occurrence of k must occur at or before 6*2^(n-1). - Robert G. Wilson v, Dec 16 2009
a(3*2^n) = n; if x = 2^n, then a(x) = a(2*x); and if x is not a power of two, then a(x) = y then a(2*x) > y. - Robert G. Wilson v, Dec 16 2009
a(n) = 0 iff n is a prime power. - Franklin T. Adams-Watters, Aug 20 2013

Programs

  • Haskell
    a089233 n = sum $ [a063524 $ gcd u v | let ds = tail $ a027750_row n,
                                           u <- ds, v <- dropWhile (<= u) ds]
    -- Reinhard Zumkeller, Sep 04 2013
    
  • Magma
    [(NumberOfDivisors(n^2)-1)/2 - NumberOfDivisors(n)+1: n in [1..100]]; // Vincenzo Librandi, Dec 23 2018
  • Mathematica
    a[n_] := (DivisorSigma[0, n^2] - 1)/2 - DivisorSigma[0, n] + 1; Array[a, 104] (* Robert G. Wilson v, Dec 16 2009 *)
  • PARI
    a(n) = (numdiv(n^2)-1)/2 - numdiv(n) + 1; \\ Michel Marcus, Feb 17 2016
    

Formula

a(n) = #{(x,y): 1 < x < y, x|n, y|n and gcd(x, y) = 1}.
a(n) = A063647(n) - A000005(n) + 1.
a(n) = A018892(n) - A000005(n). - Franklin T. Adams-Watters, Aug 20 2013

A095904 Triangular array of natural numbers (greater than 1) arranged by prime signature.

Original entry on oeis.org

2, 3, 4, 5, 9, 6, 7, 25, 10, 8, 11, 49, 14, 27, 12, 13, 121, 15, 125, 18, 16, 17, 169, 21, 343, 20, 81, 24, 19, 289, 22, 1331, 28, 625, 40, 30, 23, 361, 26, 2197, 44, 2401, 54, 42, 32, 29, 529, 33, 4913, 45, 14641, 56, 66, 243, 36, 31, 841, 34, 6859, 50, 28561, 88, 70
Offset: 0

Views

Author

Alford Arnold, Jul 10 2004

Keywords

Comments

The unit, 1, has the empty prime signature { } (thus not in triangle).
Downwards diagonals:
* Rightmost diagonal: smallest numbers of a given prime signature in increasing order (A025487). This defines the order of signatures used.
This special ordering of prime signatures (by increasing smallest numbers of a given prime signature, A181087) is unrelated to any of the 8 variants of graded lexicographic or colexicographic orderings (based on the exponents only) since it depends on the magnitudes of the prime numbers. It is not even graded by Omega(n).
* Second rightmost diagonal: second smallest numbers of a given prime signature (A077560). (They are not increasing anymore.)
Upwards diagonals:
* Leftmost diagonal: primes. {1} (A000040)
* 2nd leftmost diagonal: squares of primes. {2} (A001248)
* 3rd leftmost diagonal: squarefree biprimes. {1,1} (A006881)
* 4th leftmost diagonal: cubes of primes. {3} (A030078)
* 5th leftmost diagonal: signature (Achilles numbers) {1,2} (A054753)
* 6th leftmost diagonal: fourth powers of primes. {4} (A030514)
* 7th leftmost diagonal: signature (Achilles numbers) {1,3} (A065036)
* 8th leftmost diagonal: squarefree triprimes. {1,1,1} (A007304)
The Achilles numbers are nonsquarefree while not perfect powers.
Prime signatures are often expressed in increasing order of exponents. The decreasing order of exponents (as on the Wiki page, see links) has the advantage of listing the exponents in the same order (with the canonical factorization convention) as the smallest number of a given prime signature.

Examples

			343 is in the 4th left- and 4th rightmost diagonal, because it is the 4th value with the 4th prime signature {3}.
First 8 rows of triangular array (Cf. table link for this sequence):
                                   2
                              3         4
                         5         9         6
                    7        25        10         8
               11       49        14        27        12
          13      121        15       125        18        16
     17       169       21       343        20        81        24
19       289       22       1331       28       625        40        30
		

Crossrefs

Extensions

Extended by Ray Chandler, Jul 31 2004
Corrected (minor) by Daniel Forgues, Jan 21 2011
Example, comments by Daniel Forgues, Jan 21 2011
Edited by Alois P. Heinz, Jan 23 2011
Edited by Daniel Forgues, Jan 23 2011

A179668 Products of the 8th power of a prime and a distinct prime (p^8*q).

Original entry on oeis.org

768, 1280, 1792, 2816, 3328, 4352, 4864, 5888, 7424, 7936, 9472, 10496, 11008, 12032, 13122, 13568, 15104, 15616, 17152, 18176, 18688, 20224, 21248, 22784, 24832, 25856, 26368, 27392, 27904, 28928, 32512, 32805, 33536, 35072, 35584, 38144, 38656, 40192
Offset: 1

Views

Author

Keywords

Crossrefs

Programs

  • Mathematica
    f[n_]:=Sort[Last/@FactorInteger[n]]=={1,8}; Select[Range[40000], f]
    With[{nn=40},Take[Union[#[[1]]^8 #[[2]]&/@Flatten[Permutations/@Subsets[ Prime[Range[nn]],{2}],1]],nn]] (* Harvey P. Dale, Jan 20 2016 *)
  • PARI
    list(lim)=my(v=List(),t);forprime(p=2,(lim\2)^(1/8),t=p^8;forprime(q=2,lim\t,if(p==q,next);listput(v,t*q)));vecsort(Vec(v)) \\ Charles R Greathouse IV, Jul 20 2011
    
  • Python
    from sympy import primepi, primerange, integer_nthroot
    def A179668(n):
        def bisection(f,kmin=0,kmax=1):
            while f(kmax) > kmax: kmax <<= 1
            kmin = 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(primepi(x//p**8) for p in primerange(integer_nthroot(x,8)[0]+1))+primepi(integer_nthroot(x,9)[0])
        return bisection(f,n,n) # Chai Wah Wu, Feb 21 2025

A179692 Numbers of the form p^9*q where p and q are distinct primes.

Original entry on oeis.org

1536, 2560, 3584, 5632, 6656, 8704, 9728, 11776, 14848, 15872, 18944, 20992, 22016, 24064, 27136, 30208, 31232, 34304, 36352, 37376, 39366, 40448, 42496, 45568, 49664, 51712, 52736, 54784, 55808, 57856, 65024, 67072, 70144, 71168, 76288, 77312, 80384, 83456
Offset: 1

Views

Author

Keywords

Crossrefs

Programs

  • Mathematica
    f[n_]:=Sort[Last/@FactorInteger[n]]=={1,9}; Select[Range[90000], f]
  • PARI
    list(lim)=my(v=List(),t);forprime(p=2, (lim\2)^(1/9), t=p^9;forprime(q=2, lim\t, if(p==q, next);listput(v,t*q))); vecsort(Vec(v)) \\ Charles R Greathouse IV, Jul 24 2011
    
  • Python
    from sympy import primepi, integer_nthroot, primerange
    def A179692(n):
        def bisection(f,kmin=0,kmax=1):
            while f(kmax) > kmax: kmax <<= 1
            kmin = 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(primepi(x//p**9) for p in primerange(integer_nthroot(x,9)[0]+1))+primepi(integer_nthroot(x,10)[0])
        return bisection(f,n,n) # Chai Wah Wu, Feb 21 2025

A189990 Numbers with prime factorization p^2*q^6.

Original entry on oeis.org

576, 1600, 2916, 3136, 7744, 10816, 18225, 18496, 23104, 33856, 35721, 53824, 61504, 62500, 87616, 88209, 107584, 118336, 123201, 140625, 141376, 179776, 210681, 222784, 238144, 263169, 287296, 322624, 341056, 385641, 399424, 440896, 470596
Offset: 1

Views

Author

Keywords

Crossrefs

Subsequence of A137484.

Programs

  • Mathematica
    f[n_]:=Sort[Last/@FactorInteger[n]]=={2,6}; Select[Range[800000],f]
  • PARI
    list(lim)=my(v=List(),t);forprime(p=2, (lim\4)^(1/6), t=p^6;forprime(q=2, sqrt(lim\t), if(p==q, next);listput(v,t*q^2))); vecsort(Vec(v)) \\ Charles R Greathouse IV, Jul 24 2011
    
  • Python
    from math import isqrt
    from sympy import primepi, integer_nthroot, primerange
    def A189990(n):
        def bisection(f,kmin=0,kmax=1):
            while f(kmax) > kmax: kmax <<= 1
            kmin = 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(primepi(isqrt(x//p**6)) for p in primerange(integer_nthroot(x,6)[0]+1))+primepi(integer_nthroot(x,8)[0])
        return bisection(f,n,n) # Chai Wah Wu, Feb 21 2025

Formula

Sum_{n>=1} 1/a(n) = P(2)*P(6) - P(8) = A085548 * A085966 - A085968 = 0.003658..., where P is the prime zeta function. - Amiram Eldar, Jul 06 2020
a(n) = A065036(n)^2. - Chai Wah Wu, Mar 27 2025

A255231 The number of factorizations n = Product_i b_i^e_i, where all bases b_i are distinct, and all exponents e_i are distinct >=1.

Original entry on oeis.org

1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 1, 4, 1, 2, 1, 2, 1, 1, 1, 3, 2, 1, 2, 2, 1, 1, 1, 5, 1, 1, 1, 4, 1, 1, 1, 3, 1, 1, 1, 2, 2, 1, 1, 5, 2, 2, 1, 2, 1, 3, 1, 3, 1, 1, 1, 2, 1, 1, 2, 7, 1, 1, 1, 2, 1, 1, 1, 6, 1, 1, 2, 2, 1, 1, 1, 5, 4, 1, 1, 2, 1, 1, 1, 3, 1, 2, 1, 2, 1, 1, 1, 6, 1, 2, 2, 4, 1, 1, 1, 3, 1, 1, 1, 6, 1, 1, 1, 5, 1, 1, 1, 2, 2, 1, 1, 3
Offset: 1

Views

Author

Saverio Picozzi, Feb 18 2015

Keywords

Comments

Not multiplicative: a(48) = a(2^4*3) = 5 <> a(2^4)*a(3) = 4*1 = 4. - R. J. Mathar, Nov 05 2016

Examples

			From _R. J. Mathar_, Nov 05 2016: (Start)
a(4)=2: 4^1 = 2^2.
a(8)=2: 8^1 = 2^3.
a(9)=2: 9^1 = 3^2.
a(12)=2: 12^1 = 2^2*3^1.
a(16)=4: 16^1 = 4^2 = 2^2*4^1 = 2^4.
a(18)=2: 18^1 = 2*3^2.
a(20)=2: 20^1 = 2^2*5^1.
a(24)=3: 24^1 = 2^2*6^1 = 2^3*3^1.
a(32)=5: 32^1 = 2^1*4^2 = 2^2*8^1 = 2^3*4^1 = 2^5.
a(36)=4: 36^1 = 6^2 = 3^2*4^1 = 2^2*9^1.
a(48)=5: 48^1 = 3^1*4^2 = 2^2*12^1 = 2^3*6^1 = 2^4*3^1.
a(60)=2 : 60^1 = 2^2*15^1.
a(64)=7: 64^1 = 8^2 = 4^3 = 2^2*16^1 = 2^3*8^1 = 2^4*4^1 = 2^6.
a(72)=6 : 72^1 = 3^2*8^1 = 2^1*6^2 = 2^2*18^1 = 2^3*9^1 = 2^3*3^2.
(End)
		

Crossrefs

Cf. A000688 (b_i not necessarily distinct).

Programs

  • Maple
    # Count solutions for products if n = dvs_i^exps(i) where i=1..pividx are fixed
    Apiv := proc(n,dvs,exps,pividx)
        local dvscnt, expscopy,i,a,expsrt,e ;
        dvscnt := nops(dvs) ;
        a := 0 ;
        if pividx > dvscnt then
            # have exhausted the exponent list: leave of the recursion
            # check that dvs_i^exps(i) is a representation
            if n = mul( op(i,dvs)^op(i,exps),i=1..dvscnt) then
                # construct list of non-0 exponents
                expsrt := [];
                for i from 1 to dvscnt do
                    if op(i,exps) > 0 then
                        expsrt := [op(expsrt),op(i,exps)] ;
                    end if;
                end do;
                # check that list is duplicate-free
                if nops(expsrt) = nops( convert(expsrt,set)) then
                    return 1;
                else
                    return 0;
                end if;
            else
                return 0 ;
            end if;
        end if;
        # need a local copy of the list to modify it
        expscopy := [] ;
        for i from 1 to nops(exps) do
            expscopy := [op(expscopy),op(i,exps)] ;
        end do:
        # loop over all exponents assigned to the next base in the list.
        for e from 0 do
            candf := op(pividx,dvs)^e ;
            if modp(n,candf) <> 0 then
                break;
            end if;
            # assign e to the local copy of exponents
            expscopy := subsop(pividx=e,expscopy) ;
            a := a+procname(n,dvs,expscopy,pividx+1) ;
        end do:
        return a;
    end proc:
    A255231 := proc(n)
        local dvs,dvscnt,exps ;
        if n = 1 then
            return 1;
        end if;
        # candidates for the bases are all divisors except 1
        dvs := convert(numtheory[divisors](n) minus {1},list) ;
        dvscnt := nops(dvs) ;
        # list of exponents starts at all-0 and is
        # increased recursively
        exps := [seq(0,e=1..dvscnt)] ;
        # take any subset of dvs for the bases, i.e. exponents 0 upwards
        Apiv(n,dvs,exps,1) ;
    end proc:
    seq(A255231(n),n=1..120) ; # R. J. Mathar, Nov 05 2016

Formula

a(n)=1 for all n in A005117. a(n)=2 for all n in A001248 and for all n in A054753 and for all n in A085987 and for all n in A030078. a(n)=3 for all n in A065036. a(n)=4 for all n in A085986 and for all n in A030514. a(n)=5 for all n in A178739, all n in A179644 and for all n in A050997. a(n)=6 for all n in A143610, all n in A162142 and all n in A178740. a(n)=7 for all n in A030516. a(n)=9 for all n in A189988 and all n in A189987. a(n)=10 for all n in A092759. a(n) = 11 for all n in A179664. a(n)=12 for all n in A179646. - R. J. Mathar, Nov 05 2016, May 20 2017

Extensions

Values corrected. Incorrect comments removed. - R. J. Mathar, Nov 05 2016

A343443 If n = Product (p_j^k_j) then a(n) = Product (k_j + 2), with a(1) = 1.

Original entry on oeis.org

1, 3, 3, 4, 3, 9, 3, 5, 4, 9, 3, 12, 3, 9, 9, 6, 3, 12, 3, 12, 9, 9, 3, 15, 4, 9, 5, 12, 3, 27, 3, 7, 9, 9, 9, 16, 3, 9, 9, 15, 3, 27, 3, 12, 12, 9, 3, 18, 4, 12, 9, 12, 3, 15, 9, 15, 9, 9, 3, 36, 3, 9, 12, 8, 9, 27, 3, 12, 9, 27, 3, 20, 3, 9, 12, 12, 9, 27, 3, 18
Offset: 1

Views

Author

Ilya Gutkovskiy, Apr 15 2021

Keywords

Comments

Inverse Moebius transform of A056671.
a(n) depends only on the prime signature of n (see formulas). - Bernard Schott, May 03 2021

Crossrefs

Programs

  • Mathematica
    a[1] = 1; a[n_] := Times @@ ((#[[2]] + 2) & /@ FactorInteger[n]); Table[a[n], {n, 80}]
    a[n_] := Sum[If[GCD[d, n/d] == 1, DivisorSigma[0, d], 0], {d, Divisors[n]}]; Table[a[n], {n, 80}]
  • PARI
    a(n) = sumdiv(n, d, if(gcd(d, n/d)==1, numdiv(d))) \\ Andrew Howroyd, Apr 15 2021
    
  • PARI
    for(n=1, 100, print1(direuler(p=2, n, (1 + X - X^2)/(1-X)^2)[n], ", ")) \\ Vaclav Kotesovec, Feb 11 2023
    
  • Python
    from math import prod
    from sympy import factorint
    def A343443(n): return prod(e+2 for e in factorint(n).values()) # Chai Wah Wu, Feb 21 2025

Formula

a(n) = 2^omega(n) * tau_3(n) / tau(n), where omega = A001221, tau = A000005 and tau_3 = A007425.
a(n) = Sum_{d|n, gcd(d, n/d) = 1} tau(d).
From Bernard Schott, May 03 2021: (Start)
a(p^k) = k+2 for p prime, or signature [k].
a(A006881(n)) = 9 for signature [1, 1].
a(A054753(n)) = 12 for signature [2, 1].
a(A065036(n)) = 15 for signature [3, 1].
a(A085986(n)) = 16 for signature [2, 2].
a(A178739(n)) = 18 for signature [4, 1].
a(A143610(n)) = 20 for signature [3, 2].
a(A007304(n)) = 27 for signature [1, 1, 1]. (End)
Dirichlet g.f.: zeta(s)^2 * Product_{primes p} (1 + 1/p^s - 1/p^(2*s)). - Vaclav Kotesovec, Feb 11 2023
From Amiram Eldar, Sep 01 2023: (Start)
a(n) = A000005(A064549(n)).
a(n) = A363194(A348018(n)). (End)

A374459 Nonsquarefree exponentially odd numbers.

Original entry on oeis.org

8, 24, 27, 32, 40, 54, 56, 88, 96, 104, 120, 125, 128, 135, 136, 152, 160, 168, 184, 189, 216, 224, 232, 243, 248, 250, 264, 270, 280, 296, 297, 312, 328, 343, 344, 351, 352, 375, 376, 378, 384, 408, 416, 424, 440, 456, 459, 472, 480, 486, 488, 512, 513, 520, 536
Offset: 1

Views

Author

Amiram Eldar, Jul 09 2024

Keywords

Comments

First differs from A301517 at n = 1213. A301517(1213) = 12500 = 2^2 * 5^5 is not an exponentially odd number.
Numbers whose exponents in their prime factorization are all odd and at least one of them is larger than 1.
All the squarefree numbers (A005117) are exponentially odd. Therefore, the sequence of exponentially odd numbers (A268335) is a disjoint union of the squarefree numbers and this sequence.
The asymptotic density of this sequence is A065463 - A059956 = 0.096515099145... .

Examples

			8 = 2^3 is a term since 3 is odd and larger than 1.
		

Crossrefs

Intersection of A013929 (or A046099) and A268335.
Subsequence of A301517.
Subsequences: A062838 \ {1}, A065036, A102838, A113850, A113852, A179671, A190011, A335988 \ {1}.

Programs

  • Mathematica
    q[n_] := Module[{e = FactorInteger[n][[;; , 2]]}, AllTrue[e, OddQ] && ! AllTrue[e, # == 1 &]]; Select[Range[1000], q]
  • PARI
    is(k) = {my(e = factor(k)[, 2]); for(i = 1, #e, if(!(e[i] %2), return(0))); for(i = 1, #e, if(e[i] >1, return(1))); 0;}

Formula

a(n) = A268335(A374460(n)).
Sum_{n>=1} 1/a(n)^s = zeta(2*s) * (Product_{p prime} (1 + 1/p^s - 1/p^(2*s))) - zeta(s)/zeta(2*s) for s > 1.

A179689 Numbers with prime signature {7,2}, i.e., of form p^7*q^2 with p and q distinct primes.

Original entry on oeis.org

1152, 3200, 6272, 8748, 15488, 21632, 36992, 46208, 54675, 67712, 107163, 107648, 123008, 175232, 215168, 236672, 264627, 282752, 312500, 359552, 369603, 445568, 476288, 574592, 632043, 645248, 682112, 703125, 789507, 798848, 881792, 1013888
Offset: 1

Views

Author

Keywords

Crossrefs

Programs

  • Maple
    a:= proc(n) option remember; local k;
          for k from 1+ `if` (n=1, 1, a(n-1))
            while sort (map (x-> x[2], ifactors(k)[2]), `>`)<>[7, 2]
          do od; k
        end:
    seq (a(n), n=1..32);  # Alois P. Heinz, Jan 23 2011
  • Mathematica
    f[n_]:=Sort[Last/@FactorInteger[n]]=={2,7}; Select[Range[10^6], f]
  • PARI
    list(lim)=my(v=List(),t);forprime(p=2, (lim\4)^(1/7), t=p^7;forprime(q=2, sqrt(lim\t), if(p==q, next);listput(v,t*q^2))); vecsort(Vec(v)) \\ Charles R Greathouse IV, Jul 20 2011
    
  • Python
    from math import isqrt
    from sympy import primepi, integer_nthroot, primerange
    def A179689(n):
        def bisection(f,kmin=0,kmax=1):
            while f(kmax) > kmax: kmax <<= 1
            kmin = 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(primepi(isqrt(x//p**7)) for p in primerange(integer_nthroot(x,7)[0]+1))+primepi(integer_nthroot(x,9)[0])
        return bisection(f,n,n) # Chai Wah Wu, Feb 21 2025

Formula

Sum_{n>=1} 1/a(n) = P(2)*P(7) - P(9) = A085548 * A085967 - A085969 = 0.001741..., where P is the prime zeta function. - Amiram Eldar, Jul 06 2020

Extensions

Title edited by Daniel Forgues, Jan 22 2011

A179696 Numbers with prime signature {7,1,1}, i.e., of form p^7*q*r with p, q and r distinct primes.

Original entry on oeis.org

1920, 2688, 4224, 4480, 4992, 6528, 7040, 7296, 8320, 8832, 9856, 10880, 11136, 11648, 11904, 12160, 14208, 14720, 15232, 15744, 16512, 17024, 18048, 18304, 18560, 19840, 20352, 20608, 21870, 22656, 23424, 23680, 23936, 25728, 25984, 26240, 26752, 27264
Offset: 1

Views

Author

Keywords

Crossrefs

Programs

  • Maple
    a:= proc(n) option remember; local k;
          for k from 1+ `if` (n=1, 1, a(n-1))
            while sort (map (x-> x[2], ifactors(k)[2]), `>`)<>[7, 1, 1]
          do od; k
        end:
    seq (a(n), n=1..40); # Alois P. Heinz, Jan 23 2011
  • Mathematica
    f[n_]:=Sort[Last/@FactorInteger[n]]=={1,1,7}; Select[Range[30000], f]
  • PARI
    list(lim)=my(v=List(),t1,t2);forprime(p=2, (lim\6)^(1/7), t1=p^7;forprime(q=2, lim\t1, if(p==q, next);t2=t1*q;forprime(r=q+1, lim\t2, if(p==r,next);listput(v,t2*r)))); vecsort(Vec(v)) \\ Charles R Greathouse IV, Jul 20 2011
    
  • Python
    from math import isqrt
    from sympy import primerange, primepi, integer_nthroot
    def A179696(n):
        def bisection(f,kmin=0,kmax=1):
            while f(kmax) > kmax: kmax <<= 1
            kmin = 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((t:=primepi(s:=isqrt(y:=x//r**7)))+(t*(t-1)>>1)-sum(primepi(y//k) for k in primerange(1, s+1)) for r in primerange(integer_nthroot(x,7)[0]+1))+sum(primepi(x//p**8) for p in primerange(integer_nthroot(x,8)[0]+1))-primepi(integer_nthroot(x,9)[0])
        return bisection(f,n,n) # Chai Wah Wu, Mar 27 2025

Extensions

Title edited by Daniel Forgues, Jan 22 2011
Previous Showing 11-20 of 47 results. Next