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.

Showing 1-10 of 11 results. Next

A001597 Perfect powers: m^k where m > 0 and k >= 2.

Original entry on oeis.org

1, 4, 8, 9, 16, 25, 27, 32, 36, 49, 64, 81, 100, 121, 125, 128, 144, 169, 196, 216, 225, 243, 256, 289, 324, 343, 361, 400, 441, 484, 512, 529, 576, 625, 676, 729, 784, 841, 900, 961, 1000, 1024, 1089, 1156, 1225, 1296, 1331, 1369, 1444, 1521, 1600, 1681, 1728, 1764
Offset: 1

Views

Author

Keywords

Comments

Might also be called the nontrivial powers. - N. J. A. Sloane, Mar 24 2018
See A175064 for number of ways to write a(n) as m^k (m >= 1, k >= 1). - Jaroslav Krizek, Jan 23 2010
a(1) = 1, for n >= 2: a(n) = numbers m such that sum of perfect divisors of x = m has no solution. Perfect divisor of n is divisor d such that d^k = n for some k >= 1. a(n) for n >= 2 is complement of A175082. - Jaroslav Krizek, Jan 24 2010
A075802(a(n)) = 1. - Reinhard Zumkeller, Jun 20 2011
Catalan's conjecture (now a theorem) is that 1 occurs just once as a difference, between 8 and 9.
For a proof of Catalan's conjecture, see the paper by Metsänkylä. - L. Edson Jeffery, Nov 29 2013
m^k is the largest number n such that (n^k-m)/(n-m) is an integer (for k > 1 and m > 1). - Derek Orr, May 22 2014
From Daniel Forgues, Jul 22 2014: (Start)
a(n) is asymptotic to n^2, since the density of cubes and higher powers among the squares and higher powers is 0. E.g.,
a(10^1) = 49 (49% of 10^2),
a(10^2) = 6400 (64% of 10^4),
a(10^3) = 804357 (80.4% of 10^6),
a(10^4) = 90706576 (90.7% of 10^8),
a(10^n) ~ 10^(2n) - o(10^(2n)). (End)
A proper subset of A001694. - Robert G. Wilson v, Aug 11 2014
a(10^n): 1, 49, 6400, 804357, 90706576, 9565035601, 979846576384, 99066667994176, 9956760243243489, ... . - Robert G. Wilson v, Aug 15 2014

References

  • R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics. Addison-Wesley, Reading, MA, 1990, p. 66.
  • R. K. Guy, Unsolved Problems in Number Theory, Springer, 1st edition, 1981. See section D9.
  • René Schoof, Catalan's Conjecture, Springer-Verlag, 2008, p. 1.
  • N. J. A. Sloane, A Handbook of Integer Sequences, Academic Press, 1973 (includes this sequence).
  • N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).

Crossrefs

Complement of A007916.
Subsequence of A072103; A072777 is a subsequence.
Union of A075109 and A075090.
There are four different sequences which may legitimately be called "prime powers": A000961 (p^k, k >= 0), A246655 (p^k, k >= 1), A246547 (p^k, k >= 2), A025475 (p^k, k=0 and k >= 2), and which are sometimes confused with the present sequence.
First differences give A053289.

Programs

  • Haskell
    import Data.Map (singleton, findMin, deleteMin, insert)
    a001597 n = a001597_list !! (n-1)
    (a001597_list, a025478_list, a025479_list) =
       unzip3 $ (1, 1, 2) : f 9 (3, 2) (singleton 4 (2, 2)) where
       f zz (bz, ez) m
        | xx < zz = (xx, bx, ex) :
                    f zz (bz, ez+1) (insert (bx*xx) (bx, ex+1) $ deleteMin m)
        | xx > zz = (zz, bz, 2) :
                    f (zz+2*bz+1) (bz+1, 2) (insert (bz*zz) (bz, 3) m)
        | otherwise = f (zz+2*bz+1) (bz+1, 2) m
        where (xx, (bx, ex)) = findMin m  --  bx ^ ex == xx
    -- Reinhard Zumkeller, Mar 28 2014, Oct 04 2012, Apr 13 2012
    
  • Magma
    [1] cat [n : n in [2..1000] | IsPower(n) ];
    
  • Maple
    isA001597 := proc(n)
        local e ;
        e := seq(op(2,p),p=ifactors(n)[2]) ;
        return ( igcd(e) >=2 or n =1 ) ;
    end proc:
    A001597 := proc(n)
        option remember;
        local a;
        if n = 1 then
            1;
        else
            for a from procname(n-1)+1 do
                if isA001597(a) then
                    return a ;
                end if;
             end do;
        end if;
    end proc:
    seq(A001597(n),n=1..70) ; # R. J. Mathar, Jun 07 2011
    N:= 10000: # to get all entries <= N
    sort({1,seq(seq(a^b, b = 2 .. floor(log[a](N))), a = 2 .. floor(sqrt(N)))}); # Robert FERREOL, Jul 18 2023
  • Mathematica
    min = 0; max = 10^4;  Union@ Flatten@ Table[ n^expo, {expo, Prime@ Range@ PrimePi@ Log2@ max}, {n, Floor[1 + min^(1/expo)], max^(1/expo)}] (* T. D. Noe, Apr 18 2011; slightly modified by Robert G. Wilson v, Aug 11 2014 *)
    perfectPowerQ[n_] := n == 1 || GCD @@ FactorInteger[n][[All, 2]] > 1; Select[Range@ 1765, perfectPowerQ] (* Ant King, Jun 29 2013; slightly modified by Robert G. Wilson v, Aug 11 2014 *)
    nextPerfectPower[n_] := If[n == 1, 4, Min@ Table[ (Floor[n^(1/k)] + 1)^k, {k, 2, 1 + Floor@ Log2@ n}]]; NestList[ nextPerfectPower, 1, 55] (* Robert G. Wilson v, Aug 11 2014 *)
    Join[{1},Select[Range[2000],GCD@@FactorInteger[#][[All,2]]>1&]] (* Harvey P. Dale, Apr 30 2018 *)
  • PARI
    {a(n) = local(m, c); if( n<2, n==1, c=1; m=1; while( cMichael Somos, Aug 05 2009 */
    
  • PARI
    is(n)=ispower(n) || n==1 \\ Charles R Greathouse IV, Sep 16 2015
    
  • PARI
    list(lim)=my(v=List(vector(sqrtint(lim\=1),n,n^2))); for(e=3,logint(lim,2), for(n=2,sqrtnint(lim,e), listput(v,n^e))); Set(v) \\ Charles R Greathouse IV, Dec 10 2019
    
  • Python
    from sympy import perfect_power
    def ok(n): return n==1 or perfect_power(n)
    print([m for m in range(1, 1765) if ok(m)]) # Michael S. Branicky, Jan 04 2021
    
  • Python
    import sympy
    class A001597() :
        def _init_(self) :
            self.a = [1]
        def at(self, n):
            if n <= len(self.a):
                return self.a[n-1]
            else:
                cand = self.at(n-1)+1
                while sympy.perfect_power(cand) == False:
                    cand += 1
                self.a.append(cand)
                return cand
    a001597 = A001597()
    for n in range(1,20):
        print(a001597.at(n)) # R. J. Mathar, Mar 28 2023
    
  • Python
    from sympy import mobius, integer_nthroot
    def A001597(n):
        def f(x): return int(n-2+x+sum(mobius(k)*(integer_nthroot(x,k)[0]-1) for k in range(2,x.bit_length())))
        kmin, kmax = 1,2
        while f(kmax) >= kmax:
            kmax <<= 1
        while True:
            kmid = kmax+kmin>>1
            if f(kmid) < kmid:
                kmax = kmid
            else:
                kmin = kmid
            if kmax-kmin <= 1:
                break
        return kmax # Chai Wah Wu, Aug 13 2024
  • Sage
    def A001597_list(n) :
        return [k for k in (1..n) if k.is_perfect_power()]
    A001597_list(1764) # Peter Luschny, Feb 03 2012
    

Formula

Goldbach showed that Sum_{n >= 2} 1/(a(n)-1) = 1.
Formulas from postings to the Number Theory List by various authors, 2002:
Sum_{i >= 2} Sum_{j >= 2} 1/i^j = 1;
Sum_{k >= 2} 1/(a(k)+1) = Pi^2 / 3 - 5/2;
Sum_{k >= 2} 1/a(k) = Sum_{n >= 2} mu(n)(1- zeta(n)) approx = 0.87446436840494... See A072102.
For asymptotics see Newman.
For n > 1: gcd(exponents in prime factorization of a(n)) > 1, cf. A124010. - Reinhard Zumkeller, Apr 13 2012
a(n) ~ n^2. - Thomas Ordowski, Nov 04 2012
a(n) = n^2 - 2*n^(5/3) - 2*n^(7/5) + (13/3)*n^(4/3) - 2*n^(9/7) + 2*n^(6/5) - 2*n^(13/11) + o(n^(13/11)) (Jakimczuk, 2012). - Amiram Eldar, Jun 30 2023

Extensions

Minor corrections from N. J. A. Sloane, Jun 27 2010

A136141 Decimal expansion of Sum_{p prime} 1/(p*(p-1)).

Original entry on oeis.org

7, 7, 3, 1, 5, 6, 6, 6, 9, 0, 4, 9, 7, 9, 5, 1, 2, 7, 8, 6, 4, 3, 6, 7, 4, 5, 9, 8, 5, 5, 9, 4, 2, 3, 9, 5, 6, 1, 8, 7, 4, 1, 3, 3, 6, 0, 8, 3, 1, 8, 6, 0, 4, 8, 3, 1, 1, 0, 0, 6, 0, 6, 7, 3, 5, 6, 7, 0, 9, 0, 2, 8, 4, 8, 9, 2, 3, 3, 3, 9, 7, 8, 3, 3, 7, 9, 8, 7, 5, 8, 8, 2, 3, 3, 2, 0, 8, 1, 8, 3, 2, 8, 9
Offset: 0

Views

Author

R. J. Mathar, Mar 09 2008

Keywords

Comments

Excess of prime factors with multiplicity over distinct prime factors for random (large) integers. - Charles R Greathouse IV, Sep 06 2011
Sum of reciprocals of (proper) prime powers. The sum of reciprocals of all proper powers is A072102. - Charles R Greathouse IV, Apr 24 2012
Decimal expansion of the infinite sum of the reciprocals of the prime powers which are not prime (A246547). - Robert G. Wilson v, May 13 2019
See the second 'Applications' example under the Mathematica help file for the function PrimePowerQ. - Robert G. Wilson v, May 13 2019
It easy to prove that this constant < 1 (Sum_{p prime} 1/(p*(p-1)) < Sum_{k>=2} 1/(k*(k-1)) = 1). Luthar (1969) asks for a better upper bound. The solution shows that this constant is < 3/2 - log(2) = 0.80685... . - Amiram Eldar, Feb 14 2025

Examples

			Equals 1/2 + 1/(3*2) + 1/(5*4) + 1/(7*6) + ...
= 0.7731566690497951278643674598559423956187413360831860483110060673567...
		

References

  • Henri Cohen, Number Theory, Volume II: Analytic and Modern Tools, GTM Vol. 240, Springer, 2007; see pp. 208-209.
  • Steven R. Finch, Mathematical Constants, Cambridge Univ. Press, 2003, Meissel-Mertens constants, p. 94.

Crossrefs

Cf. A152447 (over the semiprimes), A000040, A000720, A001248, A046660 (excess, see first comment), A072102, A077761, A083342, A179119, A246547.
Decimal expansion of the prime zeta function: A085548 (at 2), A085541 (at 3), A085964 (at 4) to A085969 (at 9).

Programs

  • Magma
    R := RealField(105);
    c := &+[R|(EulerPhi(n)-MoebiusMu(n))/n*Log(ZetaFunction(R,n)):n in[2..360]];
    Reverse(IntegerToSequence(Floor(c*10^103))); // Jason Kimberley, Jan 12 2017
  • Mathematica
    digits = 103; sp = NSum[PrimeZetaP[n], {n, 2, Infinity}, WorkingPrecision -> digits + 10, NSumTerms -> 2*digits]; RealDigits[sp, 10, digits] // First (* Jean-François Alcover, Sep 02 2015 *)
  • PARI
    W(x)=solve(y=log(x)/2,max(1,log(x)),y*exp(y)-x)
    eps()=2. >> (32*ceil(default(realprecision)/9.63))
    primezeta(s)=my(t=s*log(2),iter=W(t/eps())\t);sum(k=1,iter, moebius(k)/k*log(abs(zeta(k*s))))
    a(lim,e)={ \\ choose parameters to maximize speed and precision
        my(x,y=exp(W(lim)-.5));
        x=lim^e*(e*log(y))^e*(y*log(y))^-e*incgam(-e,e*log(y));
        forprime(p=2,lim,x+=1/((p*1.)^e*(p-1)));
        x+sum(n=2,e,primezeta(n))
    }; \\ Charles R Greathouse IV, Sep 07 2011
    
  • PARI
    sumeulerrat(1/(p*(p-1))) \\ Amiram Eldar, Mar 18 2021
    

Formula

Equals Sum_{n>=1} 1/A036689(n).
Equals Sum_{s>=2} P(s), where P is the prime zeta function. - Charles R Greathouse IV, Sep 06 2011
Equals A083342 - A077761, that is, Sum_{n>=2} ((EulerPhi(n) - MoebiusMu(n))/n) * log(zeta(n)). - Jean-François Alcover, Sep 02 2015
Equals 2 * Sum_{k>=2} pi(k)/(k^3-k), where pi(k) = A000720(k) (Shamos, 2011, p. 8). - Amiram Eldar, Mar 12 2024

Extensions

More terms from D. S. McNeil, Sep 06 2011
More digits from Jean-François Alcover, Sep 02 2015

A131605 Perfect powers of nonprimes (m^k where m is a nonprime positive integer and k >= 2).

Original entry on oeis.org

1, 36, 100, 144, 196, 216, 225, 324, 400, 441, 484, 576, 676, 784, 900, 1000, 1089, 1156, 1225, 1296, 1444, 1521, 1600, 1728, 1764, 1936, 2025, 2116, 2304, 2500, 2601, 2704, 2744, 2916, 3025, 3136, 3249, 3364, 3375, 3600, 3844, 3969, 4225, 4356, 4624
Offset: 1

Views

Author

Daniel Forgues, May 27 2008

Keywords

Comments

Although 1 is a square, is a cube, and so on..., 1 is sometimes excluded from perfect powers since it is not a well-defined power of 1 (1 = 1^k for any k in [2, 3, 4, 5, ...])
From Michael De Vlieger, Aug 11 2025: (Start)
This sequence is A001597 \ A246547, i.e., perfect powers without proper prime powers.
Union of {1} with the intersection of A001597 and A126706, where A126706 is the sequence of numbers that are neither prime powers nor squarefree.
Union of {1} and A286708 \ A052486, i.e., powerful numbers that are not prime powers, without Achilles numbers, but including the empty product. (End)

Crossrefs

Programs

  • Mathematica
    With[{nn = 2^20}, {1}~Join~Select[Union@ Flatten@ Table[a^2*b^3, {b, Surd[nn, 3]}, {a, Sqrt[nn/b^3]}], And[Length[#2] > 1, GCD @@ #2 > 1] & @@ {#, FactorInteger[#][[;; , -1]]} &] ] (* Michael De Vlieger, Aug 11 2025 *)
  • PARI
    isok(n) = if (n == 1, return (1), return (ispower(n, ,&np) && (! isprime(np)))); \\ Michel Marcus, Jun 12 2013
    
  • Python
    from sympy import mobius, integer_nthroot, primepi
    def A131605(n):
        def f(x): return int(n-2+x+sum(mobius(k)*((a:=integer_nthroot(x,k)[0])-1)+primepi(a) for k in range(2,x.bit_length())))
        kmin, kmax = 1,2
        while f(kmax) >= kmax:
            kmax <<= 1
        while True:
            kmid = kmax+kmin>>1
            if f(kmid) < kmid:
                kmax = kmid
            else:
                kmin = kmid
            if kmax-kmin <= 1:
                break
        return kmax # Chai Wah Wu, Aug 14 2024

Formula

Sum_{n>=1} 1/a(n) = 1 + A072102 - A136141 = 1.10130769935514973882... . - Amiram Eldar, Aug 15 2025

A052486 Achilles numbers - powerful but imperfect: if n = Product(p_i^e_i) then all e_i > 1 (i.e., powerful), but the highest common factor of the e_i is 1, i.e., not a perfect power.

Original entry on oeis.org

72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, 1125, 1152, 1323, 1352, 1372, 1568, 1800, 1944, 2000, 2312, 2592, 2700, 2888, 3087, 3200, 3267, 3456, 3528, 3872, 3888, 4000, 4232, 4500, 4563, 4608, 5000, 5292, 5324, 5400, 5408, 5488, 6075
Offset: 1

Views

Author

Henry Bottomley, Mar 16 2000

Keywords

Comments

Number of terms < 10^n: 0, 1, 13, 60, 252, 916, 3158, 10553, 34561, 111891, 359340, 1148195, 3656246, 11616582, 36851965, ..., A118896(n) - A070428(n). - Robert G. Wilson v, Aug 11 2014
a(n) = (s(n))^2 * f(n), s(n) > 1, f(n) > 1, where s(n) is not a power of f(n), and f(n) is squarefree and gcd(s(n), f(n)) = f(n). - Daniel Forgues, Aug 11 2015

Examples

			a(3)=200 because 200=2^3*5^2, both 3 and 2 are greater than 1, and the highest common factor of 3 and 2 is 1.
Factorizations of a(1) to a(20):
    72 = 2^3  3^2,  108 = 2^2 3^3,  200 = 2^3 5^2,  288 = 2^5  3^2,
   392 = 2^3  7^2,  432 = 2^4 3^3,  500 = 2^2 5^3,  648 = 2^3  3^4,
   675 = 3^3  5^2,  800 = 2^5 5^2,  864 = 2^5 3^3,  968 = 2^3 11^2,
   972 = 2^2  3^5, 1125 = 3^2 5^3, 1152 = 2^7 3^2, 1323 = 3^3  7^2,
  1352 = 2^3 13^2, 1372 = 2^2 7^3, 1568 = 2^5 7^2, 1800 = 2^3  3^2 5^2.
Examples for a(n) = (s(n))^2 * f(n): (see above comment)
s(n) = 6,  6, 10, 12, 14, 12, 10, 18, 15, 20, 12, 22, 18, 15, 24, 21,
f(n) = 2,  3,  2,  2,  2,  3,  5,  2,  3,  2,  6,  2,  3,  5,  2,  3,
		

Crossrefs

Programs

  • Maple
    filter:= proc(n) local E; E:= map(t->t[2], ifactors(n)[2]); min(E)>1 and igcd(op(E))=1 end proc:
    select(filter,[$1..10000]); # Robert Israel, Aug 11 2014
  • Mathematica
    achillesQ[n_] := Block[{ls = Last /@ FactorInteger@n}, Min@ ls > 1 == GCD @@ ls]; Select[ Range@ 5500, achillesQ@# &] (* Robert G. Wilson v, Jun 10 2010 *)
  • PARI
    is(n)=my(f=factor(n)[,2]); n>9 && vecmin(f)>1 && gcd(f)==1 \\ Charles R Greathouse IV, Sep 18 2015, replacing code by M. F. Hasler, Sep 23 2010
    
  • Python
    from math import gcd
    from itertools import count, islice
    from sympy import factorint
    def A052486_gen(startvalue=1): # generator of terms >= startvalue
        return (n for n in count(max(startvalue,1)) if (lambda x: all(e > 1 for e in x) and gcd(*x) == 1)(factorint(n).values()))
    A052486_list = list(islice(A052486_gen(),20)) # Chai Wah Wu, Feb 19 2022
    
  • Python
    from math import isqrt
    from sympy import mobius, integer_nthroot
    def A052486(n):
        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):
            c, l = n+x+1, 0
            j = isqrt(x)
            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+sum(mobius(k)*(integer_nthroot(x, k)[0]-1) for k in range(2, x.bit_length()))
            return c
        return bisection(f,n,n) # Chai Wah Wu, Sep 10 2024

Formula

a(n) = O(n^2). - Daniel Forgues, Aug 11 2015
a(n) = O(n^2 / log log n). - Daniel Forgues, Aug 12 2015
Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - Sum_{k>=2} mu(k)*(1-zeta(k)) - 1 = A082695 - A072102 - 1 = 0.06913206841581433836... - Amiram Eldar, Oct 14 2020

Extensions

Example edited by Mac Coombe (mac.coombe(AT)gmail.com), Sep 18 2010
Name edited by M. F. Hasler, Jul 17 2019

A091050 Number of divisors of n that are perfect powers.

Original entry on oeis.org

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

Views

Author

Reinhard Zumkeller, Dec 15 2003

Keywords

Comments

Not the same as A005361: a(72)=5 <> A005361(72)=6.

Examples

			Divisors of n=108: {1,2,3,4,6,9,12,18,27,36,54,108},
a(108) = #{1^2, 2^2, 3^2, 3^3, 6^2} = 5.
		

Crossrefs

Programs

  • Haskell
    a091050 = sum . map a075802 . a027750_row
    -- Reinhard Zumkeller, Dec 13 2012
    
  • Mathematica
    ppQ[n_] := GCD @@ Last /@ FactorInteger@ n > 1; ppQ[1] = True; f[n_] := Length@ Select[ Divisors@ n, ppQ]; Array[f, 105] (* Robert G. Wilson v, Dec 12 2012 *)
  • PARI
    a(n) = 1+ sumdiv(n, d, ispower(d)>1); \\ Michel Marcus, Sep 21 2014
    
  • PARI
    a(n)={my(f=factor(n)[,2]); 1 + if(#f, sum(k=2, vecmax(f), moebius(k)*(1 - prod(i=1, #f, 1 + f[i]\k))))} \\ Andrew Howroyd, Aug 30 2020

Formula

a(n) = 1 iff n is squarefree: a(A005117(n)) = 1, a(A013929(n)) > 1.
a(p^k) = k for p prime, k>0: a(A000961(n)) = A025474(n).
a(n) = Sum_{k=1..A000005(n)} A075802(A027750(n,k)). - Reinhard Zumkeller, Dec 13 2012
G.f.: Sum_{k=i^j, i>=1, j>=2, excluding duplicates} x^k/(1 - x^k). - Ilya Gutkovskiy, Mar 20 2017
Asymptotic mean: Limit_{m->oo} (1/m) * Sum_{k=1..m} a(k) = 1 + A072102 = 1.874464... . - Amiram Eldar, Dec 31 2023

Extensions

Wrong formula deleted by Amiram Eldar, Apr 29 2020

A183096 a(n) = number of divisors of n that are not perfect powers.

Original entry on oeis.org

0, 1, 1, 1, 1, 3, 1, 1, 1, 3, 1, 4, 1, 3, 3, 1, 1, 4, 1, 4, 3, 3, 1, 5, 1, 3, 1, 4, 1, 7, 1, 1, 3, 3, 3, 5, 1, 3, 3, 5, 1, 7, 1, 4, 4, 3, 1, 6, 1, 4, 3, 4, 1, 5, 3, 5, 3, 3, 1, 10, 1, 3, 4, 1, 3, 7, 1, 4, 3, 7, 1, 7, 1, 3, 4, 4, 3, 7, 1, 6, 1, 3, 1, 10, 3, 3, 3, 5, 1, 10, 3, 4, 3, 3, 3, 7, 1, 4, 4, 5
Offset: 1

Views

Author

Jaroslav Krizek, Dec 25 2010

Keywords

Comments

Sequence is not the same as A183093: a(72) = 7, A183093(72) = 6.

Examples

			For n = 12, set of such divisors is {2, 3, 6, 12}; a(12) = 4.
		

Crossrefs

Programs

  • Mathematica
    ppQ[n_] := GCD @@ FactorInteger[n][[;;, 2]] > 1; ppQ[1] = True; a[n_] := DivisorSum[n, 1 &, !ppQ[#] &]; Array[a, 100] (* Amiram Eldar, Jan 30 2025 *)
  • PARI
    A091050(n) = (1+ sumdiv(n, d, ispower(d)>1)); \\ This function from Michel Marcus, Sep 21 2014
    A183096(n) = (numdiv(n) - A091050(n)); \\ Antti Karttunen, Nov 23 2017

Formula

a(n) = A000005(n) - A091050(n).
a(1) = 0, a(p) = 1, a(pq) = 3, a(pq...z) = 2^k - 1, a(p^k) = 1, for p, q = primes, k = natural numbers, pq...z = product of k (k > 2) distinct primes p, q, ..., z.
Sum_{k=1..n} a(k) ~ n*(log(n) + 2*gamma - A072102 - 2), where gamma is Euler's constant (A001620). - Amiram Eldar, Jan 29 2025

A340585 Noncube perfect powers.

Original entry on oeis.org

4, 9, 16, 25, 32, 36, 49, 81, 100, 121, 128, 144, 169, 196, 225, 243, 256, 289, 324, 361, 400, 441, 484, 529, 576, 625, 676, 784, 841, 900, 961, 1024, 1089, 1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936, 2025, 2048, 2116, 2187, 2209, 2304, 2401, 2500
Offset: 1

Views

Author

Hugo Pfoertner, Jan 12 2021

Keywords

Comments

This was the original definition of A239870. However, the true definition of that sequence seems to be slightly different.

Crossrefs

Programs

  • Maple
    filter:= proc(n) local g;
      g:= igcd(op(ifactors(n)[2][..,2]));
      g > 1 and (g mod 3 <> 0)
    end proc:
    select(filter, [$1..10000]); # Robert Israel, Jan 12 2021
  • Mathematica
    Select[Range[2, 2500], (g = GCD @@ FactorInteger[#][[;; , 2]]) > 1 && !Divisible[g, 3] &] (* Amiram Eldar, Jan 12 2021 *)
  • PARI
    for(n=2,2500,if( ispower(n) % 3, print1(n,", ")))
    
  • Python
    from math import isqrt
    from sympy import mobius, integer_nthroot
    def A340585(n):
        def f(x): return int(n+x-isqrt(x)+sum(mobius(k)*(integer_nthroot(x,k)[0]-1) for k in range(5,x.bit_length())))
        kmin, kmax = 1,2
        while f(kmax) >= kmax:
            kmax <<= 1
        while True:
            kmid = kmax+kmin>>1
            if f(kmid) < kmid:
                kmax = kmid
            else:
                kmin = kmid
            if kmax-kmin <= 1:
                break
        return kmax # Chai Wah Wu, Aug 14 2024

Formula

Sum_{n>=1} 1/a(n) = 1 - zeta(3) + Sum_{k>=2} mu(k)*(1-zeta(k)) = 1 - A002117 + A072102 = 0.6724074652... - Amiram Eldar, Jan 12 2021

A340588 Squares of perfect powers.

Original entry on oeis.org

1, 16, 64, 81, 256, 625, 729, 1024, 1296, 2401, 4096, 6561, 10000, 14641, 15625, 16384, 20736, 28561, 38416, 46656, 50625, 59049, 65536, 83521, 104976, 117649, 130321, 160000, 194481, 234256, 262144, 279841, 331776, 390625, 456976, 531441, 614656, 707281, 810000, 923521, 1000000
Offset: 1

Views

Author

Terry D. Grant, Sep 21 2020

Keywords

Crossrefs

Cf. A153158 (complement within positive squares).

Programs

  • Maple
    q:= n-> is(igcd(seq(i[2], i=ifactors(n)[2]))<>2):
    select(q, [i^2$i=1..1000])[];  # Alois P. Heinz, Nov 26 2024
  • Mathematica
    Join[{1}, (Select[Range[2000], GCD @@ FactorInteger[#][[All, 2]] > 1 &])^2]
  • Python
    from sympy import mobius, integer_nthroot
    def A340588(n):
        def f(x): return int(n-2+x+sum(mobius(k)*(integer_nthroot(x,k)[0]-1) for k in range(2,x.bit_length())))
        kmin, kmax = 1,2
        while f(kmax) >= kmax:
            kmax <<= 1
        while True:
            kmid = kmax+kmin>>1
            if f(kmid) < kmid:
                kmax = kmid
            else:
                kmin = kmid
            if kmax-kmin <= 1:
                break
        return kmax**2 # Chai Wah Wu, Aug 14 2024

Formula

a(n) = A001597(n)^2.
a(n+1) = A062965(n) + 1. - Hugo Pfoertner, Sep 29 2020
Sum_{k>1} 1/(a(k) - 1) = 7/4 - Pi^2/6 = 7/4 - zeta(2).
Sum_{k>1} 1/a(k) = Sum_{k>=2} mu(k)*(1-zeta(2*k)).

A242977 Decimal expansion of Sum_{k>1} 1/(k*(k-1)*zeta(k)), a constant related to Niven's constant.

Original entry on oeis.org

7, 6, 6, 9, 4, 4, 4, 9, 0, 5, 2, 1, 0, 8, 8, 2, 4, 1, 6, 5, 2, 4, 1, 7, 9, 2, 3, 0, 0, 3, 1, 7, 6, 9, 3, 0, 9, 7, 4, 7, 5, 7, 8, 8, 9, 9, 3, 1, 9, 0, 5, 1, 6, 9, 6, 5, 4, 1, 2, 2, 0, 8, 1, 6, 0, 7, 8, 9, 6, 8, 4, 2, 3, 7, 5, 6, 7, 9, 5, 7, 7, 5, 7, 8, 9, 3, 7, 4, 6, 2, 9, 8, 4, 0, 9, 9, 4, 3
Offset: 0

Views

Author

Jean-François Alcover, May 28 2014

Keywords

Comments

The asymptotic mean of the reciprocals of the maximal exponent in prime factorization of the positive integers. - Amiram Eldar, Dec 15 2022

Examples

			0.766944490521088241652417923...
		

References

  • Steven R. Finch, Mathematical Constants, Cambridge University Press, 2003, Section 2.6 Niven's constant, p. 113.

Crossrefs

Programs

  • Mathematica
    digits = 98; m0 = 100; dm = 100; Clear[f]; f[m_] := f[m] = NSum[1/(k*(k - 1)*Zeta[k]), {k, 2, m}, WorkingPrecision -> digits + 10, NSumTerms -> m] + 1/m; f[m0]; f[m = m0 + dm]; While[RealDigits[f[m], 10, digits] != RealDigits[f[m - dm], 10, digits], Print["m = ", m ]; m = m + dm]; RealDigits[f[m], 10, digits] // First
  • PARI
    sumpos(k = 2, 1/(k*(k-1)*zeta(k))) \\ Amiram Eldar, Dec 15 2022

Formula

Equals lim_{n->oo} (1/n) * Sum_{k=2..n} 1/A051903(k). - Amiram Eldar, Oct 16 2020
Equals 1 + Sum_{k>=2} (1/zeta(k)-1)/(k*(k-1)). - Amiram Eldar, Dec 15 2022

A380399 The number of nonunitary divisors of n that are perfect powers (A001597).

Original entry on oeis.org

0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0
Offset: 1

Views

Author

Amiram Eldar, Jan 23 2025

Keywords

Examples

			a(16) = 2 since 16 have 2 nonunitary divisors that are perfect powers, 4 = 2^2 and 8 = 2^3.
a(32) = 3 since 32 have 3 nonunitary divisors that are perfect powers, 4 = 2^3, 8 = 2^3, and 16 = 2^4.
		

Crossrefs

Programs

  • Mathematica
    ppQ[n_] := n == 1 || GCD @@ FactorInteger[n][[;; , 2]] > 1; a[n_] := DivisorSum[n, 1 &, !CoprimeQ[#, n/#] && ppQ[#] &]; Array[a, 100]
  • PARI
    a(n) = sumdiv(n, d, gcd(d, n/d) > 1 && (d == 1 || ispower(d)));

Formula

a(n) = Sum_{d|n, gcd(d, n/d) > 1} [d in A001597], where [] is the Iverson bracket.
a(n) = A091050(n) - A380398(n).
Asymptotic mean: Limit_{m->oo} (1/m) * Sum_{k=1..m} a(k) = A072102 + Sum_{k>=2} mu(k)*(zeta(k)/zeta(k+1) - 1) = Sum_{k>=2} mu(k)*zeta(k)*(1/zeta(k+1)-1) = 0.38105110303589889319..., where mu is the Moebius function (A008683).
Showing 1-10 of 11 results. Next