(PARI)
A054747(n)={local(p=vector(n)); local(q=matrix(2, 2)); q[1, 1] = 2; q[1, 2] = 0; q[2, 1]=0; q[2, 2]=1; my(S=0, A() = sum(j=1, 2, prod(r=1, n, prod(s=1, 2, (2*sumdiv(lcm(r, s), d, if(d < n+1, d*p[d], 0)))^(p[r]*q[j, s]*gcd(r, s)))))/2,
inc()=!forstep(i=n, 1, -1, p[i]<n\i && p[i]++ && return; p[i]=0), t); until(inc(), t=0; for( i=1, n, if( n < t+=i*p[i], until(i++>n, p[i]=n); next(2))); t==n && S+ = A()/prod(i=1, n, i^p[i]*p[i]!)); S} \\ This is a modification of
M. F. Hasler's PARI program from
A002854. -
Petros Hadjicostas, Mar 08 2021
