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.

Previous Showing 31-40 of 112 results. Next

A053143 Smallest square divisible by n.

Original entry on oeis.org

1, 4, 9, 4, 25, 36, 49, 16, 9, 100, 121, 36, 169, 196, 225, 16, 289, 36, 361, 100, 441, 484, 529, 144, 25, 676, 81, 196, 841, 900, 961, 64, 1089, 1156, 1225, 36, 1369, 1444, 1521, 400, 1681, 1764, 1849, 484, 225, 2116, 2209, 144, 49, 100, 2601, 676, 2809
Offset: 1

Views

Author

Henry Bottomley, Feb 28 2000

Keywords

Crossrefs

Programs

  • Mathematica
    Table[i = 1; While[Mod[i^2, n] > 0, i++]; i^2, {n, 100}] (* T. D. Noe, Oct 30 2011 *)
    f[p_, e_] := p^(e + Mod[e, 2]); a[n_] := Times @@ (f @@@ FactorInteger[n]); Array[a, 100] (* Amiram Eldar, Aug 29 2019 *)
  • PARI
    a(n) = n*core(n); /* Joerg Arndt, Aug 02 2012 */

Formula

a(n) = n*A007913(n) = A019554(n)^2 = n^2/A008833(n) = (n/A000188(n))^2.
Multiplicative with p^e -> p^(e + e mod 2), p prime. - Reinhard Zumkeller, Feb 09 2003
Dirichlet g.f.: zeta(2s-2)*zeta(s-2)/zeta(2s-4). - R. J. Mathar, Oct 31 2011
Sum_{k=1..n} a(k) ~ Pi^2 * n^3 / 45. - Vaclav Kotesovec, Feb 08 2019
Sum_{n>=1} 1/a(n) = 5/2. - Amiram Eldar, Jul 29 2022

A326046 a(n) = gcd(n-A326039(n), A326040(n)-n).

Original entry on oeis.org

1, 1, 1, 1, 4, 2, 3, 1, 1, 1, 1, 4, 12, 2, 1, 1, 8, 1, 3, 1, 5, 2, 1, 4, 1, 5, 1, 24, 28, 6, 15, 1, 1, 1, 1, 1, 36, 2, 1, 1, 40, 2, 3, 4, 4, 10, 1, 4, 1, 7, 15, 3, 4, 2, 19, 4, 1, 1, 1, 8, 60, 2, 1, 1, 1, 6, 3, 1, 1, 2, 35, 1, 72, 1, 1, 12, 1, 2, 3, 1, 1, 1, 1, 4, 1, 2, 1, 4, 8, 27, 5, 8, 29, 2, 7, 60, 48, 1, 1, 1, 100, 6, 3, 1, 1
Offset: 1

Views

Author

Antti Karttunen, Jun 06 2019

Keywords

Crossrefs

Programs

Formula

a(n) = gcd(A326044(n), A326045(n)) = gcd(n-A326039(n), A326040(n)-n).

A331591 a(n) is the number of distinct prime factors of A225546(n), or equally, number of distinct prime factors of A293442(n).

Original entry on oeis.org

0, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 2, 2, 1, 1, 1, 1, 2, 1
Offset: 1

Views

Author

Antti Karttunen and Peter Munn, Jan 21 2020

Keywords

Comments

a(n) is the number of terms in the unique factorization of n into powers of squarefree numbers with distinct exponents that are powers of 2. See A329332 for a description of the relationship between this factorization, canonical (prime power) factorization and A225546.
The result depends only on the prime signature of n.
a(n) is the number of distinct bit-positions where there is a 1-bit in the binary representation of an exponent in the prime factorization of n. - Antti Karttunen, Feb 05 2020
The first 3 is a(128) = a(2^1 * 2^2 * 2^4) = 3 and in general each m occurs first at position 2^(2^m-1) = A058891(m+1). - Peter Munn, Mar 07 2022

Examples

			From _Peter Munn_, Jan 28 2020: (Start)
The factorization of 6 into powers of squarefree numbers with distinct exponents that are powers of 2 is 6 = 6^(2^0) = 6^1, which has 1 term. So a(6) = 1.
Similarly, 40 = 10^(2^0) * 2^(2^1) = 10^1 * 2^2 = 10 * 4, which has 2 terms. So a(40) = 2.
Similarly, 320 = 5^(2^0) * 2^(2^1) * 2^(2^2) = 5^1 * 2^2 * 2^4 = 5 * 4 * 16, which has 3 terms. So a(320) = 3.
10^100 (a googol) factorizes in this way as 10^4 * 10^32 * 10^64. So a(10^100) = 3.
(End)
		

Crossrefs

Sequences with related definitions: A001221, A331309, A331592, A331593, A331740.
Positions of records: A058891.
Positions of 1's: A340682.
Sequences used to express relationships between the terms: A007913, A008833, A059796, A331590.

Programs

  • Mathematica
    Array[PrimeNu@ If[# == 1, 1, Times @@ Flatten@ Map[Function[{p, e}, Map[Prime[Log2@ # + 1]^(2^(PrimePi@ p - 1)) &, DeleteCases[NumberExpand[e, 2], 0]]] @@ # &, FactorInteger[#]]] &, 105] (* Michael De Vlieger, Jan 24 2020 *)
    f[e_] := Position[Reverse[IntegerDigits[e, 2]], 1] // Flatten; a[n_] := CountDistinct[Flatten[f /@ FactorInteger[n][[;; , 2]]]]; a[1] = 0; Array[a, 100] (* Amiram Eldar, Dec 23 2023 *)
  • PARI
    A331591(n) = if(1==n,0,my(f=factor(n),u=#binary(vecmax(f[, 2])),xs=vector(u),m=1,e); for(i=1,u,for(k=1,#f~, if(bitand(f[k,2],m),xs[i]++)); m<<=1); #select(x -> (x>0),xs));
    
  • PARI
    A331591(n) = if(1==n, 0, hammingweight(fold(bitor, factor(n)[, 2]))); \\ Antti Karttunen, Feb 05 2020
    
  • PARI
    A331591(n) = if(n==1, 0, (core(n)>1) + A331591(core(n,1)[2])) \\ Peter Munn, Mar 08 2022

Formula

a(n) = A001221(A293442(n)) = A001221(A225546(n)).
From Peter Munn, Jan 28 2020: (Start)
a(n) = A000120(A267116(n)).
a(n) = a(A007913(n)) + a(A008833(n)).
For m >= 2, a(A005117(m)) = 1.
a(n^2) = a(n).
(End)
a(n) <= A331740(n) <= A048675(n) <= A293447(n). - Antti Karttunen, Feb 05 2020
From Peter Munn, Mar 07 2022: (Start)
a(n) <= A299090(n).
a(A337533(n)) = A299090(A337533(n)).
a(A337534(n)) < A299090(A337534(n)).
max(a(n), a(k)) <= a(A059796(n, k)) = a(A331590(n, k)) <= a(n) + a(k).
(End)

A055071 Largest square dividing n!.

Original entry on oeis.org

1, 1, 1, 4, 4, 144, 144, 576, 5184, 518400, 518400, 2073600, 2073600, 101606400, 914457600, 14631321600, 14631321600, 526727577600, 526727577600, 52672757760000, 52672757760000, 6373403688960000, 6373403688960000, 917770131210240000, 22944253280256000000
Offset: 1

Views

Author

Labos Elemer, Jun 13 2000

Keywords

Crossrefs

Programs

  • Maple
    seq(expand(numtheory[nthpow](n!, 2)), n=1..26); # Peter Luschny, Apr 03 2013
  • Mathematica
    a[n_] := Select[Reverse @ Divisors[n!], IntegerQ[Sqrt[#]] &, 1] // First; a /@ Range[23] (* Jean-François Alcover, May 19 2011 *)
    f[p_, e_] := p^(2*Floor[e/2]); a[n_] := Times @@ (f @@@ FactorInteger[n!]); Array[a, 30] (* Amiram Eldar, Jul 26 2024 *)
  • PARI
    a(n)=core(n!,2)[2]^2 \\ Charles R Greathouse IV, Apr 04 2012
    
  • Python
    from math import prod
    from itertools import count, islice
    from collections import Counter
    from sympy import factorint
    def A055071_gen(): # generator of terms
        c = Counter()
        for i in count(1):
            c += Counter(factorint(i))
            yield prod(p**(e-(e&1)) for p, e in c.items())
    A055071_list = list(islice(A055071_gen(),30)) # Chai Wah Wu, Jul 27 2024

Formula

a(n) = A008833(n!).
log a(n) ~ n log n. - Charles R Greathouse IV, Apr 04 2012
a(n) = A055772(n)^2. - Amiram Eldar, Jul 26 2024

Extensions

More terms from James Sellers, Jun 20 2000

A056056 Square root of largest square dividing n-th central binomial coefficient.

Original entry on oeis.org

1, 1, 1, 1, 1, 2, 1, 1, 3, 6, 1, 2, 2, 2, 3, 3, 1, 2, 1, 2, 2, 2, 1, 2, 10, 10, 30, 30, 6, 12, 3, 3, 3, 6, 5, 10, 10, 10, 3, 6, 2, 2, 2, 2, 30, 60, 15, 30, 42, 42, 42, 42, 14, 28, 2, 2, 2, 4, 2, 4, 4, 4, 21, 21, 7, 14, 7, 14, 6, 6, 1, 2, 2, 2, 10, 10, 70, 140, 7, 14, 126, 126, 6, 6, 30, 60
Offset: 1

Views

Author

Labos Elemer, Jul 26 2000

Keywords

Crossrefs

Programs

  • Mathematica
    Table[Sqrt@ Max@ Select[Divisors@ Binomial[n, Floor[n/2]], IntegerQ@ Sqrt@ # &], {n, 0, 86}] (* Michael De Vlieger, Jul 04 2016 *)
    a[n_] := Times @@ (First[#]^Floor[Last[#]/2] & /@ FactorInteger[Binomial[n, Floor[n/2]]]); Array[a, 100] (* Amiram Eldar, Sep 06 2020 *)
  • PARI
    a(n) = b = binomial(n, n\2); sqrtint(b/core(b)); \\ Michel Marcus, Dec 10 2013

Formula

a(n) = A000188(A001405(n)).

A056059 GCD of largest square and squarefree part of central binomial coefficients.

Original entry on oeis.org

1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 3, 6, 2, 1, 1, 1, 3, 3, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 6, 3, 1, 1, 1, 2, 3, 6, 2, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 2, 1, 1, 1, 1, 2, 3, 6, 6, 3, 1, 2, 2, 1, 2, 1, 3, 6, 1, 1, 1, 2, 1, 2
Offset: 1

Views

Author

Labos Elemer, Jul 26 2000

Keywords

Examples

			n=14, binomial(14,7) = 3432 = 8*3*11*13. The largest square divisor is 4, squarefree part is 858. So a(14) = gcd(4,858) = 2.
		

Crossrefs

Programs

  • Mathematica
    Table[GCD[First@ Select[Reverse@ Divisors@ #, IntegerQ@ Sqrt@ # &], Times @@ Power @@@ Map[{#1, Mod[#2, 2]} & @@ # &, FactorInteger@ #]] &@ Binomial[n, Floor[n/2]], {n, 80}] (* Michael De Vlieger, Feb 18 2017, after Zak Seidov at A007913 *)
  • PARI
    A001405(n) = binomial(n, n\2);
    A055229(n) = { my(c=core(n)); gcd(c, n/c); } \\ Charles R Greathouse IV, Nov 20 2012
    A056059(n) = A055229(A001405(n)); \\ Antti Karttunen, Jul 20 2017
    
  • Python
    from sympy import binomial, gcd
    from sympy.ntheory.factor_ import core
    def a001405(n): return binomial(n, n//2)
    def a055229(n):
        c = core(n)
        return gcd(c, n//c)
    def a(n): return a055229(a001405(n))
    print([a(n) for n in range(1, 151)]) # Indranil Ghosh, Jul 20 2017

Formula

a(n) = A055229(A001405(n)), where A055229(n) = gcd(A008833(n), A007913(n)).

Extensions

Formula clarified by Antti Karttunen, Jul 20 2017

A064179 Infinitary version of Moebius function: infinitary MoebiusMu of n, equal to mu(n) iff mu(n) differs from zero, else 1 or -1 depending on whether the sum of the binary digits of the exponents in the prime decomposition of n is even or odd.

Original entry on oeis.org

1, -1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1
Offset: 1

Views

Author

Wouter Meeussen, Sep 20 2001

Keywords

Comments

Apparently the (ordinary) Dirichlet inverse of A050377. - R. J. Mathar, Jul 15 2010
Also analog of Liouville's function (A008836) in Fermi-Dirac arithmetic, where the terms of A050376 play the role of primes (see examples). - Vladimir Shevelev, Oct 28 2013.

Examples

			G.f. = x - x^2 - x^3 - x^4 - x^5 + x^6 - x^7 + x^8 - x^9 + x^10 - x^11 + x^12 + ...
mu[45]=0 but iMoebiusMu[45]=1 because 45 = 3^2 * 5^1 and the binary digits of 2 and 1 add up to 2, an even number.
A unique representation of 48 over distinct terms of A050376 is 3*16. Since it contains even factors, then a(48)=1; for 54 such a representation is 2*3*9, thus a(54)=-1. - _Vladimir Shevelev_, Oct 28 2013
		

References

  • Vladimir S. Shevelev, Multiplicative functions in the Fermi-Dirac arithmetic, Izvestia Vuzov of the North-Caucasus region, Nature sciences 4 (1996), 28-43 (in Russian)

Crossrefs

Sequences with related definitions: A008683, A008836, A064547, A302777.
Positions of -1: A000028.
Positions of 1: A000379.
Sequences used to express relationships between the terms: A000188, A003961, A007913, A008833, A059897, A225546.

Programs

  • Mathematica
    iMoebiusMu[n_] := Switch[MoebiusMu[n], 1, 1, -1, -1, 0, If[OddQ[Plus@@(DigitCount[Last[Transpose[FactorInteger[n]]], 2, 1])], -1, 1]];
    (* The Moebius inversion formula seems to hold for iMoebiusMu and the infinitary_divisors of n: if g[ n_ ] := Plus@@(f/@iDivisors[ n ]) for all n, then f[ n_ ]===Plus@@(iMoebiusMu[ # ]g[ n/# ]&/@iDivisors[ n ]) *)
    f[p_, e_] := (-1)^DigitCount[e, 2, 1]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* Amiram Eldar, Dec 23 2023 *)
  • PARI
    {a(n) = my(A, p, e); if( n<1, 0, A = factor(n); prod(k=1, matsize(A)[1], [p, e] = A[k, ]; (-1) ^ subst( Pol( binary(e)), x, 1)))}; /* Michael Somos, Jan 08 2008 */
    
  • PARI
    a(n) = if (n==1, 1, (-1)^omega(core(n)) * a(core(n,1)[2])) \\ Peter Munn, Mar 16 2022
    
  • PARI
    a(n) = vecprod(apply(x -> (-1)^hammingweight(x), factor(n)[, 2])); \\ Amiram Eldar, Dec 23 2023
    
  • Python
    from math import prod
    from sympy import factorint
    def A064179(n): return prod(-1 if e.bit_count()&1 else 1 for e in factorint(n).values()) # Chai Wah Wu, Oct 12 2024
  • Scheme
    (define (A064179 n) (expt -1 (A064547 n))) ;; Antti Karttunen, Nov 23 2017
    

Formula

From Vladimir Shevelev Feb 20 2011: (Start)
Sum_{d runs through i-divisors of n} a(d)=1 if n=1, or 0 if n>1; Sum_{d runs through i-divisors of n} a(d)/d = A091732(n)/n.
Infinitary Moebius inversion:
If Sum_{d runs through i-divisors of n} f(d)=F(n), then f(n) = Sum_{d runs through i-divisors of n} a(d)*F(n/d). (End)
a(n) = (-1)^A064547(n). - R. J. Mathar, Apr 19 2011
Let k=k(n) be the number of terms of A050376 that divide n with odd maximal exponent. Then a(n) = (-1)^k. For example, if n=96, then the maximal exponent of 2 that divides 96 is 5, for 3 it is 1, for 4 it is 2, for 16 it is 1. Thus k(96)=3 and a(96)=-1. - Vladimir Shevelev, Oct 28 2013
From Peter Munn, Jan 25 2020: (Start)
a(A050376(n)) = -1; a(A059897(n,k)) = a(n) * a(k).
a(n^2) = a(n).
a(A003961(n)) = a(n).
a(A225546(n)) = a(n).
a(A000028(n)) = -1; a(A000379(n)) = 1.
(End)
a(n) = a(A007913(n)) * a(A008833(n)) = (-1)^A001221(A007913(n)) * a(A000188(n)). - Peter Munn, Mar 16 2022
From Amiram Eldar, Dec 23 2023: (Start)
Multiplicative with a(p^e) = (-1)^A000120(e).
Dirichlet g.f.: 1/Product_{k>=0} zeta(2^k * s) (Steuding et al., 2011). (End)

A252895 Numbers with an odd number of square divisors.

Original entry on oeis.org

1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 16, 17, 19, 21, 22, 23, 26, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 41, 42, 43, 46, 47, 48, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73, 74, 77, 78, 79, 80, 81, 82, 83, 85, 86, 87, 89, 91, 93, 94, 95, 96, 97
Offset: 1

Views

Author

Walker Dewey Anderson, Mar 22 2015

Keywords

Comments

Open lockers in the locker problem where the student numbers are the set of perfect squares.
The locker problem is a classic mathematical problem. Imagine a row containing an infinite number of lockers numbered from one to infinity. Also imagine an infinite number of students numbered from one to infinity. All of the lockers begin closed. The first student opens every locker that is a multiple of one, which is every locker. The second student closes every locker that is a multiple of two, so all of the even-numbered lockers are closed. The third student opens or closes every locker that is a multiple of three. This process continues for all of the students. [This is sometimes called the light switch problem - see A360845.]
A variant on the locker problem is when not all student numbers are considered; in the case of this sequence, only the square-numbered students open and close lockers. The sequence here is a list of the open lockers after all of the students have gone.
n is in the sequence if and only if it is the product of a squarefree number (A005117) and a fourth power (A000583). - Robert Israel, Apr 07 2015
Let D be the multiset containing d0(k), the divisor counting function, for each divisor k of n. n is in the sequence if and only if D admits a partition into two parts A and B such that the sum of the elements of A is exactly one more or less than the sum of the elements of B. For example, if n = 80, we have D = {1, 2, 2, 3, 4, 4, 5, 6, 8, 10}, and A = {1, 2, 3, 4, 4, 8} and B = {2, 5, 6, 10}. The sum of A is 22, and the sum of B is 23. - Griffin N. Macris, Oct 10 2016
From Amiram Eldar, Jul 07 2020: (Start)
Numbers k such that the largest square dividing k (A008833) is a fourth power.
The asymptotic density of this sequence is Pi^2/15 = A182448 = 0.657973... (Cesàro, 1885). (End)
Closed under the binary operation A059897(.,.), forming a subgroup of the positive integers under A059897. - Peter Munn, Aug 01 2020

Examples

			The set of divisors of 6 is {1,2,3,6}, which contains only one perfect square: 1; therefore 6 is a term.
The set of divisors of 16 is {1,2,4,8,16}, which contains three perfect squares: 1, 4, and 16; therefore 16 is a term.
The set of divisors of 4 is {1,2,4}, which contains two perfect squares: 1 and 4; therefore 4 is not a term.
		

Crossrefs

Positions of ones in A335324.

Programs

  • Haskell
    a252895 n = a252895_list !! (n-1)
    a252895_list = filter (odd . a046951) [1..]
    -- Reinhard Zumkeller, Apr 06 2015
  • Maple
    N:= 1000: # to get all terms <= N
    S:= select(numtheory:-issqrfree, {$1..N}):
    map(s -> seq(s*i^4, i = 1 .. floor((N/s)^(1/4))), S);
    # if using Maple 11 or earlier, uncomment the next line
    # sort(convert(%,list)); # Robert Israel, Apr 07 2015
  • Mathematica
    Position[Length@ Select[Divisors@ #, IntegerQ@ Sqrt@ # &] & /@ Range@ 70, Integer?OddQ] // Flatten (* _Michael De Vlieger, Mar 23 2015 *)
    a[n_] := DivisorSigma[0, Total[EulerPhi/@Select[Sqrt[Divisors[n]], IntegerQ]]]; Flatten[Position[a/@Range@100,?OddQ]] (* _Ivan N. Ianakiev, Apr 07 2015 *)
    Select[Range@ 100, OddQ@ Length@ DeleteCases[Divisors@ #, k_ /; ! IntegerQ@ Sqrt@ k] &] (* Michael De Vlieger, Oct 10 2016 *)
  • PARI
    isok(n) = sumdiv(n, d, issquare(d)) % 2; \\ Michel Marcus, Mar 22 2015
    
  • Sage
    [n for n in [1..200] if len([x for x in divisors(n) if is_square(x)])%2==1] # Tom Edgar, Mar 22 2015
    

A056057 The largest square which divides n-th central binomial coefficient.

Original entry on oeis.org

1, 1, 1, 1, 1, 4, 1, 1, 9, 36, 1, 4, 4, 4, 9, 9, 1, 4, 1, 4, 4, 4, 1, 4, 100, 100, 900, 900, 36, 144, 9, 9, 9, 36, 25, 100, 100, 100, 9, 36, 4, 4, 4, 4, 900, 3600, 225, 900, 1764, 1764, 1764, 1764, 196, 784, 4, 4, 4, 16, 4, 16, 16, 16, 441, 441, 49, 196, 49, 196, 36, 36, 1, 4
Offset: 1

Views

Author

Labos Elemer, Jul 26 2000

Keywords

Crossrefs

Programs

  • Mathematica
    Table[First@ Select[Reverse@ Divisors@ Binomial[n, Floor[n/2]], IntegerQ@ Sqrt@ # &], {n, 72}] (* Michael De Vlieger, Feb 18 2017 *)
    a[n_] := Times @@ (First[#]^(2*Floor[Last[#]/2]) & /@ FactorInteger[Binomial[n, Floor[n/2]]]); Array[a, 100] (* Amiram Eldar, Sep 06 2020 *)

Formula

a(n) = A008833(A001405(n)).
a(A046098(n)) = 1.

A056626 Number of non-unitary square divisors of n.

Original entry on oeis.org

0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0
Offset: 1

Views

Author

Labos Elemer, Aug 08 2000

Keywords

Examples

			n = p^u prime power has u+1 square divisors of which 2 (i.e., 1 and n) are unitary but u-1 are not unitary, so a(p^u) = u - 1. E.g., n = 4^4 = 256, has 5 square divisors {1, 4, 16, 64, 256} of which {4, 16, 64} are not unitary, so a(256)=3.
		

Crossrefs

Programs

  • Mathematica
    Table[DivisorSum[n, 1 &, And[IntegerQ@ Sqrt@ #, ! CoprimeQ[#, n/#]] &], {n, 105}] (* Michael De Vlieger, Jul 28 2017 *)
    f1[p_, e_] := 1 + Floor[e/2]; f2[p_, e_] := 2^(1 - Mod[e, 2]); a[1] = 0; a[n_] := Times @@ f1 @@@ (fct = FactorInteger[n])- Times @@ f2 @@@ fct; Array[a, 100] (* Amiram Eldar, Sep 26 2022 *)
  • PARI
    a(n) = {my(f = factor(n), r=0, m = 0); prod(i=1,#f~,f[i,2]>>1 + 1) - 2^(omega(f) - omega(core(f)))} \\ David A. Corneth, Jul 28 2017
    
  • PARI
    a(n) = sumdiv(n, d, if(gcd(d, n/d)!=1, issquare(d))); \\ Michel Marcus, Jul 29 2017
    
  • Python
    from math import prod
    from sympy import factorint
    def A056626(n):
        f = factorint(n).values()
        return prod((e>>1)+1 for e in f)-(1<Chai Wah Wu, Aug 04 2024

Formula

a(n) = A046951(n) - 2^r(n), where r(n) is the number of distinct prime factors of the largest unitary square divisor of n. [Corrected by Amiram Eldar, Aug 03 2024]
a(n) = A046951(n) - 2^(A162641(n)). - David A. Corneth, Jul 28 2017
From Amiram Eldar, Sep 26 2022: (Start)
a(n) = A046951(n) - A056624(n).
Asymptotic mean: Limit_{m->oo} (1/m) * Sum_{k=1..m} a(k) = zeta(2)*(1 - 1/zeta(3)) = 0.27650128922802056073... . (End)

Extensions

a(32) and a(96) corrected by Michael De Vlieger, Jul 29 2017
Previous Showing 31-40 of 112 results. Next