It is ancient knowledge that every non-negative integer can be rewritten as the sum of four squared integers. For example the number 1 can be expressed as $0^2+0^2+0^2+1^2$. Or, in general, for any non-negative integer n, there exist integers a,b,c,d such that

$$n = a^2+b^2+c^2+d^2$$

Joseph-Louis Lagrange proved this in the 1700s and so it is often called Lagrange's Theorem.

This is sometimes discussed in relation to Quaternions - a type of number discovered by William Hamilton in the 1800s, represented as $$w+xtextbfi+ytextbfj+ztextbfk$$ where w,x,y,z are real numbers, and i j and k are distinct imaginary units that don't multiply associatively. Specifically, it is discussed in relation to squaring each component of the Quaternion $$w^2+x^2+y^2+z^2$$This quantity is sometimes called the norm, or squared norm, or also Quadrance. Some modern proofs of LaGrange's Theorem use Quaternions.

Rudolf Lipschitz studied Quaternions with only integer components, called Lipschitz Quaternions. Using Quadrance, we can imagine that every Lipschitz Quaternion can be thought of having a friend in the integers. For example Quaternion $0+0i+0j+1k$ can be thought of as associated with the integer $1=0^2+0^2+0^2+1^2$. Also, if we go backwards, then every integer can be thought of as having a friend in the Lipschitz Quaternions.

But there is an interesting detail of Lagrange's theorem - the summation is not unique. Each integer may have several different sets of four squares that can be summed to create it. For example, the number 1 can be expressed in 4 ways using non-negative integers (let us only consider non-negatives for this challenge):
$$1=0^2+0^2+0^2+1^2$$
$$1=0^2+0^2+1^2+0^2$$
$$1=0^2+1^2+0^2+0^2$$
$$1=1^2+0^2+0^2+0^2$$

The summands are always squares of 0, or 1, but they can be in different positions in the expression.

For this challenge, let us also "sort" our summands lowest to highest, to eliminate duplicates, so that we could consider, for this exercise, that 1 only has 1 way of being represented as the sum of four squares:

$$1=0^2+0^2+0^2+1^2$$

Another example is the number 42, which can be expressed in four ways (again, only considering non-negative a,b,c,d, and eliminating duplicate component arrangements)

$$42=0^2+1^2+4^2+5^2$$
$$42=1^2+1^2+2^2+6^2$$
$$42=1^2+3^2+4^2+4^2$$
$$42=2^2+2^2+3^2+5^2$$

What if we imagine each of these different ways of expressing an integer as being associated to a specific Quaternion? Then we could say the number 42 is associated with these four Quaternions:

$$0+1i+4j+5k$$
$$1+1i+2j+6k$$
$$1+3i+4j+4k$$
$$2+2i+3j+5k$$

If we imagine the standard computer graphics interpretation of a Quaternion, where $i$ $j$ and $k$ are vectors in three dimensional Euclidian space, and so the $x$ $y$ and $z$ components of the Quaternion are 3 dimensional Cartesian coordinates, then we can imagine that each Integer, through this thought process, can be associated with a set of 3 dimensional coordinates in space. For example, the number 42 is associated with the following four $(x,y,z)$ coordinates: $$(1,4,5),(1,2,6),(3,4,4),(2,3,5)$$

This can be thought of as a point cloud, or a set of points in space. Now, one interesting thing about a set of finite points in space is that you can always draw a minimal bounding box around them - a box that is big enough to fit all the points, but no bigger. If you imagine the box as being an ordinary box aligned with the x,y,z axes, it is called an Axis Aligned Bounding Box. The bounding box also has a volume, calculable by determining its width, length, and height, and multiplying them together.

We can then imagine the volume of a bounding box for the points formed by our Quaternions. For the integer 1, we have, using the criteria of this exercise, one Quaternion whose Quadrance is 1, $0+0i+0j+1k$. This is a very simple point cloud, it only has one point, so it's bounding box has volume 0. For the integer 42, however, we have four Quaternions, and so four points, around which we can draw a bounding box. The minimum point of the box is (1,2,4) and the maximum is (3,4,6) resulting in a width, length, and height of 2,2, and 2, giving a volume of 8.

Let's say that for an integer n, the qvolume is the volume of the Axis Aligned Bounding Box of all the 3d points formed by Quaternions that have a Quadrance equal to n, where the components of the Quaternion $w+xi+yj+zk$ are non-negative and $w<=x<=y<=z$.

Create a program or function that, given a single non-negative integer n, will output n's qvolume.

Examples:

input -> output
0 -> 0
1 -> 0
31 -> 4
32 -> 0
42 -> 8
137 -> 96
1729 -> 10032

This is code-golf, smallest number of bytes wins.

Jelly, 17 bytes

Żœċ4²S⁼ɗƇ⁸ZḊṢ€I§P

Try it online! (pretty slow - make it fast enough for all the test cases with a leading ½)

How?

Żœċ4²S⁼ɗƇ⁸ZḊṢ€I§P - Link: non-negative integer, n e.g. 27
Ż - zero-range [0,1,2,...,27]
 4 - literal four 4
 œċ - combinations with replacement [[0,0,0,0],[0,0,0,1],...,[0,0,0,27],[0,0,1,1],[0,0,1,2],...,[27,27,27,27]]
 Ƈ - filter keep those for which: e.g.: [0,1,1,5]
 ɗ - last three links as a dyad:
 ² - square (vectorises) [0,1,1,25]
 S - sum 27
 ⁼ ⁸ - equal to? chain's left argument, n 1
 - -> [[0,1,1,5],[0,3,3,3],[1,1,3,4]]
 Z - transpose [[0,0,1],[1,3,1],[1,3,3],[5,3,4]]
 Ḋ - dequeue [[1,3,1],[1,3,3],[5,3,4]]
 Ṣ€ - sort each [[1,1,3],[1,3,3],[3,4,5]]
 I - incremental differences (vectorises) [[ 0,2 ],[ 2,0 ],[ 1,1 ]]
 § - sum each [2,2,2]
 P - product 8

Python 2, 138 bytes

q=lambda n,x=0,*t:[t]*(n==0)if t[3:]else q(n-x*x,x,x,*t)+q(n,x+1,*t+(0,)*(x>n))
p=1
for l in zip(*q(input()))[:3]:p*=max(l)-min(l)
print p

Try it online!

Recursively generates the reverse-sorted quaternions with the given norm, then takes the product between the max and min of all possible values in the first three indices.

itertools might have had a shot if it didn't use ridiculously long names like itertools.combinations_with_replacement

Python 2, 161 bytes

from itertools import*
n=input();p=1
for l in zip(*[t[1:]for t in combinations_with_replacement(range(n+1),4)if sum(x*x for x in t)==n]):p*=max(l)-min(l)
print p

Try it online!

This is why itertools is never the answer.

C# (Visual C# Interactive Compiler), 229 bytes

a=>uint b=0,c=~0U,d,e,f=0,g=0,h=0,i=c,j=c,k=c;for(;b*b<=a;b++)for(c=b;c*c<=a;c++)for(d=c;d*d<=a;d++)for(e=d;e*e<=a;e++)if(b*b+c*c+d*d+e*e==a)f=c>f?c:f;g=d>g?d:g;h=e>h?e:h;i=c<i?c:i;j=d<j?d:j;k=e<k?e:k;return(f-i)*(g-j)*(h-k);

Try it online!

Wolfram Language (Mathematica), 67 58 bytes

Volume@BoundingRegion[Rest/@PowersRepresentations[#,4,2]]&

Try it online!

 ...& Pure function:
PowersRepresentations[#,4,2] Get the sorted reprs. of # as sums of 4 2nd powers
Rest/@ Drop the first coordinate of each
BoundingRegion[...] Find the bounding region, a Cuboid[] or Point[].
 By default Mathematica finds an axis-aligned cuboid.
Volume Find volume; volume of a Point[] is 0.

JavaScript (ES6),  148  143 bytes

n=>(r=[[],[],[]]).map(a=>p*=a.length+~a.indexOf(1),(g=(s,k=0,a=[])=>a[3]?s||r.map(r=>r[a.pop()]=p=1):g(s-k*k,k,[...a,++k],k>s||g(s,k,a)))(n))|p

Try it online!

Commented

We initialize an array $r$ with 3 empty arrays.

r = [ [], [], [] ]

For each valid value of $x$, we will set to $1$ the value at $x+1$ in the first array. Ditto for $y$ and $z$ with the 2nd and 3rd arrays respectively.

The dimensions of the bounding box will be deduced from the distance between the first and the last entry set to $1$ in these arrays.

Step 1

To fill $r$, we use the recursive function $g$.

g = ( // g is a recursive function taking:
 s, // s = current sum, initially set to the input n
 k = 0, // k = next value to be squared
 a = [] // a[] = list of selected values
) => //
 a[3] ? // if we have 4 values in a[]:
 s || // if s is equal to zero (we've found a valid sum of 4 squares):
 r.map(r => // for each array r[] in r[]:
 r[a.pop()] // pop the last value from a[]
 = p = 1 // and set the corresponding value in r[] to 1
 // (also initialize p to 1 for later use in step 2)
 ) // end of map()
 : // else:
 g( // do a recursive call:
 s - k * k, // subtract k² from s
 k, // pass k unchanged
 [...a, ++k], // increment k and append it to a[]
 k > s || // if k is less than or equal to s:
 g(s, k, a) // do another recursive call with s and a[] unchanged
 ) // end of outer recursive call

Step 2

We can now compute the product $p$ of the dimensions.

r.map(a => // for each array a[] in r[]:
 p *= // multiply p by:
 a.length + // the length of a[]
 ~a.indexOf(1) // minus 1, minus the index of the first 1 in a[]
) | p // end of map(); return p

Sledgehammer, 13 bytes

⣚⡩⣏⢪⠨⣹⣓⢩⣒⡀⡸⢀⣺

Use with this version of Sledgehammer. (This was pushed after the challenge was posted, but the changes actually increase the byte count.)

Decompresses to the Mathematica expression Volume[BoundingRegion[Rest /@ PowersRepresentations[Slot[], 4, 2]]] which is an optimized version of my Mathematica answer. However, Rest[] still takes 1.75 bytes to encode.