A340313 The n-th squarefree number is the a(n)-th squarefree number having its number of primes.

1, 1, 2, 3, 1, 4, 2, 5, 6, 3, 4, 7, 8, 5, 6, 9, 7, 10, 1, 11, 8, 9, 10, 12, 11, 12, 13, 2, 14, 13, 15, 14, 16, 15, 16, 17, 17, 18, 18, 19, 3, 19, 20, 4, 20, 21, 21, 22, 5, 22, 23, 23, 24, 25, 26, 24, 27, 28, 29, 30, 25, 26, 6, 27, 7, 31, 28, 29, 8, 32, 30, 9

Offset: 1

Author

Peter Dolland, Jan 04 2021

Keywords

• nonn

Comments

The sequence gives the column index of A005117(n) in the array A340316 and may be understood as a complementary addition to A072047 giving the row index.

Examples

			{x|x <= 6, A072047(x) = A072047(6) = 1} = {2,3,4,6}, therefore a(6) = 4.
{x|x <= 28, A072047(x) = A072047(28) = 3} = {19,28}, therefore a(28) = 2.
		

Links

• David A. Corneth, Table of n, a(n) for n = 1..10000
• Alois P. Heinz, Plot of n, a(n) for n = 1..1000000

Crossrefs

Cf. A001221, A001222, A005117 (squarefree numbers), A058933, A067003, A072047 (number of prime factors), A340316 (squarefree numbers array).

Programs

Haskell

a340313 n = a340313_list !! (n-1)
a340313_list = repetitions a072047_list
    where
    repetitions [] = []
    repetitions (a:as) = 1 : h a as (repetitions as)
    h  []  = []
    h b (c:cs) (r:rs) = (if c == b then succ else id) r : h b cs rs

Maple

with(numtheory):
b:= proc(n) option remember; local k; if n=1 then 1 else
      for k from 1+b(n-1) while not issqrfree(k) do od; k fi
    end:
p:= proc() 0 end:
a:= proc(n) option remember; local h; a(n-1);
      h:= bigomega(b(n)); p(h):= p(h)+1;
    end: a(0):=0:
seq(a(n), n=1..100);  # Alois P. Heinz, Jan 06 2021

Mathematica

b[n_] := b[n] = Module[{k}, If[n == 1, 1,
     For[k = 1 + b[n - 1], !SquareFreeQ[k], k++]; k]];
p[_] = 0;
a[n_] := a[n] = Module[{h}, a[n - 1];
     h = PrimeOmega[b[n]]; p[h] = p[h]+1];
a[0] = 0;
Table[a[n], {n, 1, 100}] (* Jean-François Alcover, Mar 28 2022, after Alois P. Heinz *)

PARI

first(n) = {v = vector(5); n--; res = vector(n); t = 0; for(i = 2, oo, f = factor(i)[,2]; if(vecmax(f) == 1, if(#f > #v, v = concat(v, vector(#f - #v)) ); t++; v[#f]++; res[t] = v[#f]; if(t >= n, return(concat(1, res)) ) ) ) } \\ David A. Corneth, Jan 07 2021

Python

from math import isqrt, prod
from sympy import primerange, integer_nthroot, mobius, primenu, primepi
def A340313(n):
    if n == 1: return 1
    def g(x,a,b,c,m): yield from (((d,) for d in enumerate(primerange(b+1,isqrt(x//c)+1),a+1)) if m==2 else (((a2,b2),)+d for a2,b2 in enumerate(primerange(b+1,integer_nthroot(x//c,m)[0]+1),a+1) for d in g(x,a2,b2,c*b2,m-1)))
    def f(x): return n+x-sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))
    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
    kmax = bisection(f)
    return int(sum(primepi(kmax//prod(c[1] for c in a))-a[-1][0] for a in g(kmax,0,1,1,m)) if (m:=primenu(kmax)) > 1 else primepi(kmax)) # Chai Wah Wu, Aug 31 2024

Formula

a(n) = #{x|x <= n, A072047(x) = A072047(n)}.