A059379 - OEIS

%I #31 Dec 16 2023 13:11:04

%S 1,1,1,2,3,1,2,8,7,1,4,12,26,15,1,2,24,56,80,31,1,6,24,124,240,242,63,

%T 1,4,48,182,624,992,728,127,1,6,48,342,1200,3124,4032,2186,255,1,4,72,

%U 448,2400,7502,15624,16256,6560,511,1,10,72,702,3840

%N Array of values of Jordan function J_k(n) read by antidiagonals (version 1).

%D L. Comtet, Advanced Combinatorics, Reidel, 1974, p. 199, #3.

%D R. Sivaramakrishnan, "The many facets of Euler's totient. II. Generalizations and analogues", Nieuw Arch. Wisk. (4) 8 (1990), no. 2, 169-187.

%H Enrique Pérez Herrero, <a href="/A059379/b059379.txt">Table of n, a(n) for n = 1..10000</a>

%F J_k(n) = sum( d divides n, d^k*mu(n/d)) - _Benoit Cloitre_ and Michael Orrison (orrison(AT)math.hmc.edu), Jun 07 2002

%e Array begins:

%e   1,  1,  2,   2,   4,    2,    6,    4,   6,  4, 10, 4, ...

%e   1,  3,  8,  12,  24,   24,   48,   48,  72, 72, ...

%e   1,  7, 26,  56, 124,  182,  342,  448, 702, ...

%e   1, 15, 80, 240, 624, 1200, 2400, 3840, ...

%p J := proc(n,k) local i,p,t1,t2; t1 := n^k; for p from 1 to n do if isprime(p) and n mod p = 0 then t1 := t1*(1-p^(-k)); fi; od; t1; end;

%p #alternative

%p A059379 := proc(n,k)

%p     add(d^k*numtheory[mobius](n/d),d=numtheory[divisors](n)) ;

%p end proc:

%p seq(seq(A059379(d-k,k),k=1..d-1),d=2..12) ; # _R. J. Mathar_, Nov 23 2018

%t JordanTotient[n_,k_:1]:=DivisorSum[n,#^k*MoebiusMu[n/#]&]/;(n>0)&&IntegerQ[n];

%t A004736[n_]:=Binomial[Floor[3/2+Sqrt[2*n]],2]-n+1;

%t A002260[n_]:=n-Binomial[Floor[1/2+Sqrt[2*n]],2];

%t A059379[n_]:=JordanTotient[A004736[n],A002260[n]]; (* _Enrique Pérez Herrero_, Dec 19 2010 *)

%o (PARI)

%o jordantot(n,k)=sumdiv(n,d,d^k*moebius(n/d));

%o A002260(n)=n-binomial(floor(1/2+sqrt(2*n)),2);

%o A004736(n)=binomial(floor(3/2+sqrt(2*n)),2)-n+1;

%o A059379(n)=jordantot(A004736(n),A002260(n)); \\ _Enrique Pérez Herrero_, Jan 08 2011

%o (Python)

%o from functools import cache

%o def MoebiusTrans(a, i):

%o     @cache

%o     def mb(n, d = 1):

%o           return d % n and -mb(d, n % d < 1) + mb(n, d + 1) or 1 // n

%o     def mob(m, n): return mb(m // n) if m % n == 0 else 0

%o     return sum(mob(i, d) * a(d) for d in range(1, i + 1))

%o def Jrow(n, size):

%o     return [MoebiusTrans(lambda m: m ** n, k) for k in range(1, size)]

%o for n in range(1, 8): print(Jrow(n, 13))

%o # Alternatively:

%o from sympy import primefactors as prime_divisors

%o from fractions import Fraction as QQ

%o from math import prod as product

%o def J(n: int, k: int) -> int:

%o     t = QQ(pow(k, n), 1)

%o     s = product(1 - QQ(1, pow(p, n)) for p in prime_divisors(k))

%o     return (t * s).numerator  # the denominator is always 1

%o for n in range(1, 8): print([J(n, k) for k in range(1, 13)])

%o # _Peter Luschny_, Dec 16 2023

%Y See A059379 and A059380 (triangle of values of J_k(n)), A000010 (J_1), A059376 (J_3), A059377 (J_4), A059378 (J_5). Columns give A000225, A024023, A020522, A024049, A059387, etc.

%Y Main diagonal gives A067858.

%K nonn,tabl

%O 1,4

%A _N. J. A. Sloane_, Jan 28 2001