%I #13 Jun 06 2020 15:36:36
%S 1,1,2,2,2,1,3,1,3,2,4,3,2,2,1,4,1,2,3,2,7,4,2,2,8,2,1,4,4,2,3,7,3,2,
%T 5,2,2,4,6,2,5,2,11,6,4,1,4,1,6,2,4,12,2,5,1,4,6,4,2,5,6,4,1,2,3,4,17,
%U 6,2,3,6,1,5,6,1,3,4,6,6,13,1,2,4,8,4,4
%N Numbers of Pythagorean quadruples contained in the divisors of A330893(n).
%e a(7) = 3 because A330893(7)=168, and the set of divisors of 168: {1, 2, 3, 4, 6, 7, 8, 12, 14, 21, 24, 28, 42, 56, 84, 168} contains three Pythagorean quadruples {2, 3, 6, 7}, {4, 6, 12, 14} and {8, 12, 24, 28}.
%p with(numtheory):
%p for n from 3 to 1700 do :
%p d:=divisors(n):n0:=nops(d):it:=0:
%p for i from 1 to n0-3 do:
%p for j from i+1 to n0-2 do :
%p for k from j+1 to n0-1 do:
%p for m from k+1 to n0 do:
%p if d[i]^2 + d[j]^2 + d[k]^2 = d[m]^2
%p then
%p it:=it+1:
%p else
%p fi:
%p od:
%p od:
%p od:
%p od:
%p if it>0 then
%p printf(`%d, `,it):
%p else fi:
%p od:
%t nq[n_] := If[Mod[n, 6] > 0, 0, Block[{t, u, v, c = 0, d = Divisors[n], m}, m = Length@ d; Do[t = d[[i]]^2 + d[[j]]^2; Do[u = t + d[[h]]^2; If[u > n^2, Break[]]; If[Mod[n^2, u] == 0 && IntegerQ[v = Sqrt@ u] && Mod[n, v] == 0, c++], {h, j+1, m-1}], {i, m-3}, {j, i+1, m - 2}]; c]]; Select[Array[nq, 1638], # > 0 &] (* _Giovanni Resta_, May 04 2020 *)
%Y Cf. A027750, A169580, A330893
%K nonn
%O 1,3
%A _Michel Lagneau_, May 01 2020