def add(a,b): return a*b/gcd(a,b)**2 def Un(n): return matrix([[a*b/gcd(a,b)**2 for a in divisors(n) if gcd(a,n/a)==1] for b in divisors(n) if gcd(b,n/b)==1]) print("Sequence entries:", [Un(n).det() for n in range(1,19)]) def KK(d,e): return len(prime_divisors(gcd(d,e))) def chi(d,e): return (-1)**KK(d,e) def usigma(n): return sum([d for d in divisors(n) if gcd(n/d,d)==1]) def minkowskiMul(A,B): return set([a*b for a in A for b in B]) def unitary_divisors(n): return [d for d in divisors(n) if gcd(n/d,d)==1] def detN(n): un = unitary_divisors(n) PD = prime_divisors(n) Sn = list(Subsets(PD)) dn = 1 for si in Sn: sj = set(PD).difference(set(si)) Pi = prod([1+p**valuation(n,p) for p in si]) Pj = prod([1-p**valuation(n,p) for p in sj]) dn = dn*Pi*Pj return dn print("\nEmpirical check of the formula:") for n in range(1,100): print(n,detN(n),detN(n)==Un(n).det())