Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Powers of squarefree numbers that are not squarefree.
14

%I #57 Aug 19 2024 11:38:51

%S 4,8,9,16,25,27,32,36,49,64,81,100,121,125,128,169,196,216,225,243,

%T 256,289,343,361,441,484,512,529,625,676,729,841,900,961,1000,1024,

%U 1089,1156,1225,1296,1331,1369,1444,1521,1681,1764,1849

%N Powers of squarefree numbers that are not squarefree.

%C For all n exists k: a(n) = A072774(k) and A072776(k) > 1.

%C Numbers k such that every prime in the prime factorization of k is raised to the same power > 1; k is a term iff k/A007947(k)^m = 1 for some m > 1. - _David James Sycamore_, Jun 12 2024

%H Stanislav Sykora and Reinhard Zumkeller, <a href="/A072777/b072777.txt">Table of n, a(n) for n = 1..20000</a> (first 10000 terms from Reinhard Zumkeller)

%F Sum_{n>=1} 1/a(n) = Sum_{n>=2} mu(n)^2/(n*(n-1)) = Sum_{n>=2} (zeta(n)/zeta(2*n) - 1) = 0.8486338679... (A368250). - _Amiram Eldar_, Jul 22 2020

%e The number 144 = 12^2 is not a member because 12 is not squarefree.

%e 64 = 2^6 and 49 = 7^2 are members because, though not squarefree, they are powers of the squarefree numbers 2 and 7, respectively. Note that 64 is included even though it is also a square of a nonsquarefree number. - _Stanislav Sykora_, Jul 11 2014

%t Select[Range[2000], Length[u = Union[FactorInteger[#][[All, 2]]]] == 1 && u[[1]] > 1 &] (* _Jean-François Alcover_, Mar 27 2013 *)

%o (Haskell)

%o import Data.Map (singleton, findMin, deleteMin, insert)

%o a072777 n = a072777_list !! (n-1)

%o a072777_list = f 9 (drop 2 a005117_list) (singleton 4 (2, 2)) where

%o f vv vs'@(v:ws@(w:_)) m

%o | xx < vv = xx : f vv vs' (insert (bx*xx) (bx, ex+1) $ deleteMin m)

%o | xx > vv = vv : f (w*w) ws (insert (v^3) (v, 3) m)

%o where (xx, (bx, ex)) = findMin m

%o -- _Reinhard Zumkeller_, Apr 06 2014

%o (PARI) BelongsToA(n) = {my(f, k, e); if(n == 1, return(0));

%o f = factor(n); e = f[1, 2]; if(e == 1, return(0));

%o for(k = 2, #f[, 2], if(f[k, 2] != e, return(0))); return(1);}

%o Ntest(nmax, test) = {my(k = 1, n = 0, v); v = vector(nmax); while(1, n++; if(test(n), v[k] = n; k++; if(k > nmax, break)); ); return(v); }

%o a = Ntest(20000, BelongsToA) \\ Note: not very efficient. - _Stanislav Sykora_, Jul 11 2014

%o (PARI) is(n)=ispower(n,,&n) && issquarefree(n) \\ _Charles R Greathouse IV_, Oct 16 2015

%o (Python)

%o from math import isqrt

%o from sympy import mobius, integer_nthroot

%o def A072777(n):

%o def g(x): return int(sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1)))-1

%o def f(x): return n-1+x-sum(g(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length()))

%o kmin, kmax = 1,2

%o while f(kmax) >= kmax:

%o kmax <<= 1

%o while True:

%o kmid = kmax+kmin>>1

%o if f(kmid) < kmid:

%o kmax = kmid

%o else:

%o kmin = kmid

%o if kmax-kmin <= 1:

%o break

%o return kmax # _Chai Wah Wu_, Aug 19 2024

%Y Cf. A013929, A368250.

%Y Cf. A005117, subsequence of A001597 and A072774.

%Y Cf. A007947.

%K nonn,nice

%O 1,1

%A _Reinhard Zumkeller_, Jul 10 2002