(PARI) omnipnotsq(n, m)= local(a, x, j, fl=0); for(x=1, n, a=factor(x); for(j=1, omega(x), if(a[j, 2]>= m, fl=1, fl=0; break); ); if(fl&issquare(x)==0, print1(x", ")) )
(Python)
from math import isqrt
from sympy import integer_nthroot, mobius
def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+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
def f(x):
j = isqrt(x)
c, l = n+x+j, 0
while j>1:
k2 = integer_nthroot(x//j**2, 3)[0]+1
w = squarefreepi(k2-1)
c -= j*(w-l)
l, j = w, isqrt(x//k2**3)
c -= squarefreepi(integer_nthroot(x, 3)[0])-l
return c
