The qvolume of an integer

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+x\textbf{i}+y\textbf{j}+z\textbf{k}$$ where $$\w,x,y,z\$$ are real numbers, and $$\\textbf{i}, \textbf{j}\$$ and $$\\textbf{k}\$$ are distinct imaginary units that don't multiply commutatively. 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+0\textbf{i}+0\textbf{j}+1\textbf{k}\$$ 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 one 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+1\textbf{i}+4\textbf{j}+5\textbf{k}$$ $$1+1\textbf{i}+2\textbf{j}+6\textbf{k}$$ $$1+3\textbf{i}+4\textbf{j}+4\textbf{k}$$ $$2+2\textbf{i}+3\textbf{j}+5\textbf{k}$$

If we imagine the standard computer graphics interpretation of a quaternion, where $$\\textbf{i}\$$, $$\\textbf{j}\$$ and $$\\textbf{k}\$$ are vectors in three dimensional Euclidean 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+0\textbf{i}+0\textbf{j}+1\textbf{k}\$$. 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+x\textbf{i}+y\textbf{j}+z\textbf{k}\$$ 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.

• what do i need to add? i meant to indicate that smallest number of bytes would win – don bright May 27 at 18:52
• You forgot the code-golf tag, I helped you add it – Embodiment of Ignorance May 27 at 18:52
• This is a nice challenge but it would be even better IMHO if it was a bit less verbose. Also, beware of irrelevant links (I'm not saying that all your links are irrelevant, but only a few of them really bring meaningful information for the challenge, while the other ones are just distracting). – Arnauld May 28 at 0:45
• Yes, but why take only i, j, k as 3D space but not 4D space? – tsh May 28 at 3:29
• @tsh because Quaternions don't necessarily represent a 4 dimensional Euclidean space. Hamilton discovered them while searching for a way to work with 3 dimensional space. It would be possible to do a 4d version but i was pondering their use in 3d space when i made the question – don bright May 28 at 4:18

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.

• wow, i had no idea that something like PowersRepresentations would be a built in in a language. i actually thought about making a challenge to show the different ways to sum an integer as four squares but im glad i did not. – don bright May 28 at 2:18
• Lol, Mathematica even has a builtin for determining goats in an image, so having a builtin for this really doesn't surprise me. xD – Kevin Cruijssen May 28 at 6:33

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]
²             -     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


z=zipWith
(!)=foldr1.z
h n=[0..n]
f n|p<-[[c,b,a]|a<-h n,b<-h a,c<-h b,d<-h c,a^2+b^2+c^2+d^2==n]=product$z(-)(max!p)$min!p


Try it online!

Pretty straightforward solution. Brute force all the possible solutions by iterating over all the values from 0 to n (way overkill but shorter bytecount). I output the point as a list so we can use @Lynn's magic (!) operator. That operator collapses each dimension with the function on the left side so max!p returns a list of size 3 which consists of the maximums along each dimension and min!p does the same for minimum. Then we just find the minimum size in each dimension (by subtracting the min value from the max with z(-)) and multiply them together.

Thanks a lot to @Lynn for taking off 9 bytes with some folding zip magic!

• I shaved off a few bytes by forgoing the transposition in favor of some zipWith logic. 123 bytes – Lynn May 29 at 15:26

Sledgehammer 0.2, 12 bytes

⡌⢎⣟⡊⢘⣚⡏⡨⠍⠁⡇⠨


Use with Mathematica 11.2 and this version of Sledgehammer, which predates the challenge. See edit history for a version that works in version 0.3, which has a GUI and generates a Mathematica expression.

This pushes the input to the stack and calls the sequence of commands

{intLiteral[4], intLiteral[2], call["PowersRepresentations", 3], call["Thread", 1], call["Rest", 1], call["Thread", 1], call["BoundingRegion", 1], call["Volume", 1]}


which is equivalent to evaluating the following Wolfram code derived from my Wolfram Language answer:

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


Try it online!

• does this require mathematica to test it? – don bright May 28 at 4:20
• @don bright Yes, the repository has instructions. It's a work in progress so not very user-friendly yet. After running setup.wls you can test either with wolframscript or interactive_app.wls. – lirtosiast May 28 at 4:38
• @Downgoat Yes. I plan to implement a golfing library, but currently it decompresses to plain Mathematica. – lirtosiast May 28 at 7:23
• @pipe The older version should work (now that I think of it, the code is exactly the same on one older version), but I'd have to download it and run setup again. (The changes since then have mostly been writing the GUI and refactoring code with no major change in functionality.) Since this answer is shortest it seems important to prove the eligibility, so I'll do that tomorrow morning. – lirtosiast May 28 at 9:24
• can anyone else run this? i kind of would like to verify it works before giving it the ol' checkmark. – don bright May 29 at 23:16

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.

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


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!

05AB1E, 18 bytes

Ý4ãʒnOQ}€{ζ¦ε{¥O}P


Try it online!

Port of Jonathan Allan's Jelly answer.

n%i=sum[maximum[t!!i*b|t<-mapM([0..n]<$f)[0..3],sum(map(^2)t)==n,scanr1 max t==t]|b<-[-1,1]] f n=n%0*n%1*n%2  Try it online! (times out on the larger test cases) There's some strange optimizations here. To compute maximum l-minimum l for the list l of elements at a given position, it turns out shorter in context to convert them both to maxima by negating the second term: maximum l+maximum(map((-1)*))l, or equivalently sum[maximum$map(b*)l||b<-[-1,1]].

To multiply the three dimensions, it turns out shorter to just write the product f n=n%0*n%1*n%2 than to use any sort of loop. Here, n%i is the difference between the min and max of the i'th coordinate values, which are extracted with indexing !!i.

To generate the valid four-tuples, we take lists of four numbers from [0..n] whose squares sum to n and are in decreasing order. We check the reverse-sortedness of t with scanr1 max t==t, which sees if the running maximum of the reverse is itself, as Haskell doesn't have a built-in sort without a costly import. I tried various ways to recursively generate the four-tuples like in my Python answers, but they were all longer than this brute-force generate-and-filter way.