login
Determinant of the matrix d*e/gcd(d, e)^2, where d, e run through the unitary divisors of n.
2

%I #73 Nov 08 2023 05:14:57

%S 1,-3,-8,-15,-24,576,-48,-63,-80,5184,-120,14400,-168,20736,36864,

%T -255,-288,57600,-360,129600,147456,129600,-528,254016,-624,254016,

%U -728,518400,-840,110075314176,-960,-1023,921600,746496,1327104,1440000,-1368,1166400,1806336,2286144,-1680,1761205026816,-1848

%N Determinant of the matrix d*e/gcd(d, e)^2, where d, e run through the unitary divisors of n.

%C The matrix M(n) = [d*e / gcd(d, e)^2], where d, e run through the unitary divisors of n, is a Dedekind group matrix, hence its eigenvalues can be computed using characters of finite abelian groups. It can be proven that Spec(M(n)) = Product_{p|n, p prime} {1 + p^v_p(n), 1 - p^v_p(n)}. From this it follows that the determinant is never equal to 0, hence a(n) != 0 for all n. It also follows that the eigenvalues of this matrix appear with multiplicity one and each eigenvalue is an integer.

%H David A. Corneth, <a href="/A367064/b367064.txt">Table of n, a(n) for n = 1..10000</a>

%H David A. Corneth, <a href="/A367064/a367064_1.png">Formula adapted from Orges Leka</a>

%H Orges Leka, <a href="/A367064/a367064.py.txt">SageMath script for empirical verification of the formula</a>.

%H MathOverflow, <a href="https://mathoverflow.net/questions/370207">The action of the unitary divisors group on the set of divisors</a>, 2020.

%H MathOverflow, <a href="https://mathoverflow.net/questions/369941">Sum of divisors and unitary divisors as the eigenvalue and the spectral norm of some addition matrix?</a>

%H MathOverflow, <a href="https://mathoverflow.net/questions/369751">Boolean ring of unitary divisors / Structure of unitary divisors?</a>

%H Math StackExchange, <a href="https://math.stackexchange.com/questions/3799607"> Does this characteristic polynomial factor into linear factors over the integers?</a>

%F a(n) = Product_{I subset {1,...,omega(n)}} Product_{i in I} (1 + p_i^v_{p_i}(n))* Product_{j in I^c} (1 - p_j^v_{p_j}(n)) (proven), where omega(n) denotes the number of distinct prime divisors of n, I^c denotes the complement of the subset I in {1, 2, ..., omega(n)} and v_p(n) denotes the valuation of n at the prime divisors p of n.

%F If the prime factorization of n is n = Product_i p_i^e_i, then a(n) = (Product_i 1 - p_i^(2*e_i))^(2^(omega(n) - 1)). In particular, a(n) < 0 if and only if n is a prime power (A246655). - _Chai Wah Wu_, Nov 06 2023

%F By the formula of _Chai Wah Wu_ it follows that a(n) = (uphi(n)*usigma(n))^ (2^(omega(n)-1))*(-1)^(omega(n)) where uphi(n) = A047994(n) is the unitary totient function and usigma(n) = A034448(n) is the sum of unitary divisors of n. - _Orges Leka_, Nov 07 2023

%e For n = 6 the a(6) = 576, since the eigenvalues of the matrix are 12, 2, -4, -6 and so 576 = 12*2*(-4)*(-6).

%o (SageMath)

%o def unitary_divisors(n):

%o return [d for d in divisors(n) if gcd(n//d, d) == 1]

%o def Un(n):

%o un = unitary_divisors(n)

%o M = matrix([[a*b//gcd(a,b)**2 for a in un] for b in un])

%o return M.det()

%o print([Un(n) for n in range(1, 44)])

%o # Alternative:

%o def detN(n):

%o PD = set(prime_divisors(n))

%o Sn = Subsets(PD)

%o dn = 1

%o for si in Sn:

%o sj = PD.difference(si)

%o Pi = prod(1 + p**valuation(n, p) for p in si)

%o Pj = prod(1 - p**valuation(n, p) for p in sj)

%o dn = dn*Pi*Pj

%o return dn

%o print([detN(n) for n in range(1, 44)])

%o # Based on the extensions found by David A. Corneth and Chai Wah Wu:

%o def uphi(n):

%o return prod(p**(valuation(n,p))-1 for p in prime_divisors(n))

%o def usigma(n):

%o return sum(d for d in divisors(n) if gcd(n//d, d) == 1)

%o def omega(n):

%o return len(prime_divisors(n))

%o def ddet(n):

%o return (uphi(n)*usigma(n))**(2**(omega(n)-1))*(-1)**(omega(n))

%o print([ddet(n) for n in range(1, 44)])

%o (PARI)

%o a(n) = {

%o my(d = divisors(n), m);

%o d = select(x->gcd(x, n/x)==1, d);

%o m = matrix(#d, #d, i, j, d[i]*d[j]/(gcd(d[i], d[j])^2));

%o matdet(m)

%o } \\ _David A. Corneth_, Nov 04 2023

%o (PARI)

%o a(n) = {

%o my(f = factor(n), P = vector(#f~, i, f[i,1]^(2*f[i,2])), res = 1);

%o forsubset(#f~, s,

%o s = Set(s);

%o vs = vector(#s, i, 1 - P[s[i]]);

%o res*=vecprod(vs);

%o );

%o return(res);

%o } \\ _David A. Corneth_, Nov 04 2023

%o (Python)

%o from math import prod

%o from sympy import factorint

%o def A367064(n):

%o f = factorint(n)

%o return prod(1-d**(e<<1) for d,e in f.items())**(1<<len(f)-1) if n>1 else 1

%o # _Chai Wah Wu_, Nov 06 2023

%Y Cf. A001221, A034444, A034448, A047994, A246655.

%K sign

%O 1,2

%A _Orges Leka_, Nov 04 2023

%E More terms from _David A. Corneth_, Nov 04 2023