cp's OEIS Frontend

A375004 Number of ordered primitive solutions (x,y,z,w) to x*y + x*z + x*w + y*z + y*w + z*w = n with x,y,z,w >= 1.

%I A375004 #22 Jul 29 2024 06:20:53
%S A375004 0,0,0,0,0,1,0,0,4,0,0,4,6,0,4,0,12,8,0,0,16,6,12,4,12,0,16,12,24,8,0,
%T A375004 12,34,0,24,8,30,12,16,0,36,32,24,12,32,6,36,16,36,12,40,12,72,8,0,24,
%U A375004 64,24,48,32,30,24,56,12,72,8,48,24,70,24,60,32,54,24,40,12,120,62,24,24,76,24,96,32
%N A375004 Number of ordered primitive solutions (x,y,z,w) to x*y + x*z + x*w + y*z + y*w + z*w = n with x,y,z,w >= 1.
%H A375004 Chai Wah Wu, <a href="/A375004/b375004.txt">Table of n, a(n) for n = 1..10000</a>
%o A375004 (PARI) a(n) = sum(x=1, n, sum(y=1, n, sum(z=1, n, sum(w=1, n, (gcd([x, y, z, w])==1)*(x*y+x*z+x*w+y*z+y*w+z*w==n)))));
%o A375004 (Python)
%o A375004 from math import gcd
%o A375004 from sympy import divisors, integer_nthroot
%o A375004 def A375004(n):
%o A375004     k = 0
%o A375004     for c in range(1,n-1):
%o A375004         for d in divisors(c,generator=True):
%o A375004             for x in range(1,d):
%o A375004                 y = d-x
%o A375004                 xy = x*y
%o A375004                 a = (c//d)**2
%o A375004                 b = a-(n-c-xy<<2)
%o A375004                 if b>=0:
%o A375004                     q,r = integer_nthroot(b,2)
%o A375004                     if r:
%o A375004                         w = c//d+q>>1
%o A375004                         z = c//d-w
%o A375004                         if 1<=w<c//d and gcd(x,y,z,w)==1:
%o A375004                             k += 1
%o A375004                         if q:
%o A375004                             w = c//d-q>>1
%o A375004                             z = c//d-w
%o A375004                             if 1<=w<c//d and gcd(x,y,z,w)==1:
%o A375004                                 k += 1
%o A375004     return k # _Chai Wah Wu_, Jul 27 2024
%Y A375004 Cf. A066958, A374970, A375003.
%K A375004 nonn
%O A375004 1,9
%A A375004 _Seiichi Manyama_, Jul 27 2024