cp's OEIS Frontend

A059404 Numbers with different exponents in their prime factorizations.

%I A059404 #56 Apr 16 2025 04:33:16
%S A059404 12,18,20,24,28,40,44,45,48,50,52,54,56,60,63,68,72,75,76,80,84,88,90,
%T A059404 92,96,98,99,104,108,112,116,117,120,124,126,132,135,136,140,144,147,
%U A059404 148,150,152,153,156,160,162,164,168,171,172,175,176,180,184,188,189,192,198,200
%N A059404 Numbers with different exponents in their prime factorizations.
%C A059404 Former name: Numbers k such that k/(largest power of squarefree kernel of k) is larger than 1.
%C A059404 Also numbers k = p(1)^alpha(1)* ... * p(r)^alpha(r) such that RootMeanSquare(alpha(1), ..., alpha(r)) is not an integer. - _Ctibor O. Zizka_, Sep 19 2008
%H A059404 Donald Alan Morrison, <a href="/A059404/b059404.txt">Table of n, a(n) for n = 1..10000</a>
%H A059404 Donald Alan Morrison, <a href="/A059404/a059404.txt">Sage program</a>
%F A059404 A062760(a(n)) > 1, i.e., a(n)/(A007947(a(n))^A051904(a(n))) = a(n)/A062759(a(n)) > 1.
%F A059404 A071625(a(n)) > 1. - _Michael S. Branicky_, Sep 01 2022
%F A059404 Sum_{n>=1} 1/a(n)^s = zeta(s) - 1 - Sum_{k>=1} (zeta(k*s)/zeta(2*k*s)-1) for s > 1. - _Amiram Eldar_, Mar 20 2025
%e A059404 440 is in the sequence because 440 = 2^3*5*11 and it has two distinct exponents, 3 and 1.
%p A059404 isA := n -> 1 < nops({seq(padic:-ordp(n, p), p in NumberTheory:-PrimeFactors(n))}): select(isA, [seq(1..190)]);  # _Peter Luschny_, Apr 14 2025
%t A059404 A059404Q[n_] := Length[Union[FactorInteger[n][[All, 2]]]] > 1;
%t A059404 Select[Range[200], A059404Q] (* _Paolo Xausa_, Jan 07 2025 *)
%o A059404 (PARI) is(n)=#Set(factor(n)[,2])>1 \\ _Charles R Greathouse IV_, Sep 18 2015
%o A059404 (Python)
%o A059404 from sympy import factorint
%o A059404 def ok(n): return len(set(factorint(n).values())) > 1
%o A059404 print([k for k in range(190) if ok(k)]) # _Michael S. Branicky_, Sep 01 2022
%o A059404 (Python)
%o A059404 from math import isqrt
%o A059404 from sympy import mobius, integer_nthroot
%o A059404 def A059404(n):
%o A059404     def g(x): return int(sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1)))
%o A059404     def f(x): return n+1-(y:=x.bit_length())+sum(g(integer_nthroot(x,k)[0]) for k in range(1,y))
%o A059404     kmin, kmax = 1,2
%o A059404     while f(kmax) >= kmax:
%o A059404         kmax <<= 1
%o A059404     while True:
%o A059404         kmid = kmax+kmin>>1
%o A059404         if f(kmid) < kmid:
%o A059404             kmax = kmid
%o A059404         else:
%o A059404             kmin = kmid
%o A059404         if kmax-kmin <= 1:
%o A059404             break
%o A059404     return kmax # _Chai Wah Wu_, Aug 19 2024
%o A059404 (SageMath)
%o A059404 def isA(n): return 1 < len(set(valuation(n, p) for p in prime_divisors(n)))
%o A059404 print([n for n in range(1, 190) if isA(n)])  # _Peter Luschny_, Apr 14 2025
%Y A059404 Complement of A072774 (powers of squarefree numbers).
%Y A059404 Cf. A003557, A007947, A051904, A062759, A062760, A071625.
%K A059404 nonn,easy
%O A059404 1,1
%A A059404 _Labos Elemer_, Jul 18 2001