Uiua, 13 bytes

/×≠7◿8[⍥⊸÷₄∞]

Test cases

Explanation

Uses the same approach as the Python answer.

/×≠7◿8[⍥⊸÷₄∞]

       ⍥   ∞  # Iterate until fixed point (0)
        ⊸÷₄   # Divide by four
      [     ] # Collect intermediate results into array
    ◿8        # Take all mod 8
/×≠7          # All are not equal to seven

SAKO, 79 bytes

PODPROGRAM:F(N)
**)LINIIENT(4*A×(8×B+7)-N)
POWTORZ:A=0(1)N
POWTORZ:B=0(1)N
WROC

Outputs 0 carriage returns for truthy and a positive number of carriage returns for falsy. An easy way to check how many were printed is by using this command: ./filename | grep $'\r' -c.

Full programme version, 81 bytes

CZYTAJ:N
**1)LINIIENT(4*A×(8×B+7)-N)
POWTORZ:A=0(1)N
POWTORZ:B=0(1)N
STOP1
KONIEC

Thumb-2 assembly; 18 bytes, 8 insns

f:
    rbit r1, r0       // 4 (0xfa90f1a0)
    clz  r1, r1       // 4 (0xfab1f181)
    movs r2, #7       // 2 (0x2207)
    lsrs r0, r0, r1   // 2 (0x40c8)
    ands r0, r0, r2   // 2 (0x4010)
    subs r0, r0, r2   // 2 (0x1a80)
    bx lr             // 2 (0x4770)

This is a fairly straightforward translation of another answer I just posted but using the Thumb-2 32-bit instruction set. This is a mixed-width instruction set, with instructions either 2 or 4 bytes wide. As it's cheaper to do operations with a register (2B) rather than an immediate (4B), I sacrifice r2 and use it to store 0x7 instead of just using it as an immediate in the ands and subs, saving two bytes. Even though the -s suffixed instructions wastefully set flags, they're all 2 bytes instead of 4, so the tradeoff is worth it.

It should handle numbers up to 2^32 - 1.

64-bit Aarch64 assembly, 24 bytes / 6 insns

f:
    rbit  x9, x0
    clz   x9, x9
    lsr   x9, x0, x9
    and   x9, x9, #7
    sub   x0, x9, #7
    ret

It follows a normal calling convention and can be called from regular C code as if it were extern int f(int);

All Aarch64 instructions are 4 bytes/32 bits long, so byte count = instruction count * 4 for all instructions. Like most of the rest of the answers I used Legendre's theorem. The code returns a non-zero value if the number can be represented as a sum of three squares, and 0 if it cannot. Adding 4 more bytes and replacing sub x0, x1, #7 with cmp x1, #7; cset x0, ne can guarantee that it returns 1 instead of a number of values, but those negatives are considered truthy as far as C is concerned, so I decided to cut 4 bytes there.

As Aarch64 doesn't include a "count trailing zeros" instruction, the input number is reversed and leading zeros are counted, then the input is shifted right by that amount, a modulo-8 is performed (or, equivalently, & 0b111), and 7 is subtracted from the remainder to either give a falsy zero, or a truthy negative integer.

It handles positive integers up to 2^64 - 1

Vyxal M, 7 bytes

3ÞẊ²Ṡ=a

Try it Online!

Here is a pretty compact solution (inspired by an answer to a similar challenge)

Explanation:

3ÞẊ     # make 3d cartesian product of 0..input
   ²    # square
    Ṡ   # sum
     =a # check whether any value of the list is equal to the input

C (gcc), 32 bytes

f(a){return!a|a%4?a%8-7:f(a/4);}

Try it online!

It should be ok for 0 too. Thank you to @ceilingcat for suggestion from || to |.

Python 3, 31 bytes

lambda n:(c:=n&-n)%3==1>~n//c%8

Try it online!

This uses Legendre's three-square theorem, that \$n\$ is a sum of 3 squares unless it has the form \$n=4^a(8b+7)\$. We output True for such n, and False otherwise.

We first extract c, the largest power-of-2 divisor in n, using the classic bit trick n&-n. Check that this c is \$1 \bmod 3\$, which indicates that c is a power of 4, the \$4^a\$ above, and not double one. This also short-circuit eliminates n==0 which would divide-by-0 error in the next part.

Next, we check that n//c%8==7. This is equivalent to ~(n//c)%8<1, and it turns out ~n//c%8<1 works too.

29 bytes

lambda n:(c:=n&-n)%3%2>~n&7*c

Try it online!

Python 2, 27 bytes

lambda n:(n&-n&n/2&n/4)%3%2

Try it online!

Thanks to @Albert.Lang for this shorter version. This bitwise approach is similar one found on the anarchy golf challenge by tails and ported to Python by mitchs.

As before, n&-n extracts the largest power-of-2 factor of n. Then, we & it with n/2 and n/4, which zeroes it unless the two bits to its left are also set. This checks that it's the rightmost bit of a 111, which indicates that the value is \$7 \bmod 8\$ once powers of 2 are removed.

It now remains to check that the number of zero bits to the right is even, that is, the 1 we identified is in an even bit position. We could & it with a bit pattern of ...01010101, which we can make as 4**n/3 (the anarchy solution uses 85 which works suffices for n<=1000). But as Albert.Lang points out, it's shorter to use %3 as in the Python 3 code to tell if the power of 2 we extracted is a power of 4, via %3%2. This also works to output 0 (False) if the &n/2&n/4 zeroed our bit, and handles n=0.

A Python 3 version comes out to 29 bytes, with two /'s becoming //.

Bespoke, 251 bytes

total I do with squares is an easy task,when applying adrian legendre^s theorem
using math,when dividing repeatedly by four
then,with quotient,adding together a numeral one
imagine the quantity matches;when matching,minusing eights produces a single O

Output is 0 for true, and 1 for false. This program uses Legendre's three-square theorem, and is descriptive of its own algorithm.

C (GCC) (Compile time), 36 bytes

n[-(((f&-f)%3&1)&&!(f/(f&-f)%8-7))];

Try it online! Based on @xnor's Python answer. If the number is representable as the sum of three squares, the program compiles, otherwise, the program will not compile. Takes input via a preprocessor macro f.

Maple, 57 bytes

n->orseq(iquo(n,4^i,'r')mod 8=7 and r=0,i=0..ilog[4](n));

Testing for 4^a*(8*b+7) as @xnor proposed. Can be shorter (48 bytes) if replace ilog[4](n) by n, but then won't handle largest n example. Outputs true if n cannot be represented as sum of three squares; false otherwise.

R, 34 bytes

f=\(n)!n||`if`(n%%4,n%%8<7,f(n/4))

Attempt This Online!

Port of l4m2's JavaScript answer.

R, 34 bytes

\(n,k=2^(0:n)^2)2^n%in%(k%o%k%o%k)

Attempt This Online!

Adapting Robin Ryder's answer to similar challenge.

R, 40 38 bytes

\(n)n%in%combn(rep(0:n,3),3,\(r)r%*%r)

Attempt This Online!

Straightforward brute force.

x86 32Bits, 44 bytes

00000838  56                push esi
00000839  8B442408          mov eax,[esp+0x8]
0000083D  31C9              xor ecx,ecx
0000083F  B104              mov cl,0x4
00000841  89C6              mov esi,eax
00000843  31D2              xor edx,edx
00000845  09C0              or eax,eax
00000847  7406              jz 0x84f
00000849  F7F1              div ecx
0000084B  09D2              or edx,edx
0000084D  74F2              jz 0x841
0000084F  89F0              mov eax,esi
00000851  31D2              xor edx,edx
00000853  31C9              xor ecx,ecx
00000855  B108              mov cl,0x8
00000857  F7F1              div ecx
00000859  31C0              xor eax,eax
0000085B  40                inc eax
0000085C  49                dec ecx
0000085D  39CA              cmp edx,ecx
0000085F  7502              jnz 0x863
00000861  31C0              xor eax,eax
00000863  5E                pop esi
00000864  C3                ret

Optimized for number of bytes, not for number of instructions... The stack at the begin of the function has to have the return address in esp+4, and than the value of input in esp+8. The caller has to do something as "push value; call abovef; add esp, 4". Possible errors in something.

APL(NARS), 18 chars

{0=4∣⍵:∇⍵÷4⋄7≠8∣⍵}

The input of above should be >0. Test:

  (⍳100)/⍨∼f ¨ ⍳100
7 15 23 28 31 39 47 55 60 63 71 79 87 92 95
  f 78454578451231202548745194641511x
0

that above is the recursive version of iterative and better one below...

APL(NARS), 48 chars

r←f w
→2×⍳∼(w≠0)∧0=4∣w⋄w÷←4⋄→1
r←1⋄→0×⍳7≠8∣w⋄r←0

//6+25+17=48 It use the theorem

"exist x,y,zeN0 : n=x^2+y^2+z^2 <=> exist a,beN0 : n=(4^a)*(8b+7)"

Because gcd(8b+7,4)=1 because 8b+7 is odd, than I divide n for 4 until find that is not divisible for 4. That number if is of type 8b+7 or mod 8 = 7, return 0 else return 1.

x86 32-bit machine code, 13 bytes

0F BC C8 D3 E8 24 07 D1 E9 1C 07 D6 C3

Try it online!

Following the regparm(1) calling convention, this function takes a 32-bit integer in EAX and returns an 8-bit integer in AL, which is -1 if the input number is a sum of three squares and 0 if it is not.

This uses Legendre's three-square theorem.

In assembly:

f:  bsf ecx, eax    # Set ECX to the lowest position of a 1 bit in EAX.
    shr eax, cl     # Shift EAX right by that much, removing the ending 0s.
    and al, 7       # Bitwise-AND AL (low byte of EAX) with 7.
    shr ecx, 1      # Shift ECX right by 1 bit. The last bit goes into CF.
    sbb al, 7       # Subtract 7+CF from AL, giving a negative value and CF=1
                    #  unless AL=7 and CF=0, in which case CF=0 afterwards.
    .byte 0xD6      # Undocumented SALC instruction -- sets AL to -CF.
    ret             # Return.

05AB1E, 9 7 bytes

Ý3ãnOIå

-2 bytes using my answer for the related Sum of two squares challenge.

Try it online or verify all smaller test cases.

Here is an additional version that's much faster (but also longer): 12 bytes

(&©÷>8Ö®3%Θ*

Also uses Legendre's three-square theorem taken from @xnor's Python answer.

Outputs an inverted 05AB1E boolean (1 if falsey; 0 or a positive integer if truthy).

Try it online or verify all test cases.

Explanation:

Ý        # Push a list in the range [0, (implicit) input]
 3ã      # Get all possible triplets using these values
   n     # Square each inner value
    O    # Sum each inner list
     Iå  # Check whether this list contains the input-integer
         # (which is output implicitly as result)

(         # Negate the (implicit) input-integer
 &        # Bitwise-AND it with the (implicit) input
  ©       # Store this value in variable `®` (without popping)
   ÷      # Integer-divide the (implicit) input by this n&-n
    >     # Increase it by 1
     8Ö   # Check whether it's now divisible by 8
  ®       # Push n&-n from variable `®` again
   3%     # Modulo-3
       *  # Multiply it to the earlier check
          # (only 1 is truthy in 05AB1E)
          # (after which it is output implicitly as result)

Ruby, 23 bytes

->n{~n&7*(n&=-n)<n%3%2}

Try it online!

Ported from xnor's python solution

Ruby, 48 bytes

f=->n,r=3{r<1?n==0:(0..n).any?{|x|f[n-x*x,r-1]}}

Try it online!

Naive brute-force solution.

Wolfram Language (Mathematica), 61 15 bytes

3~SquaresR~#>0&

Try it online!

Tips from @Neil and @att would have taken my original solution to 35 bytes with Solve[a^2+b^2+c^2==#,Integers]!={}&

but all credit to @alephalpha for pointing out this built in.

For your convenience, the regex test harnesses print numbers aren't a sum of 3 squares, since most numbers are.

Regex (Perl / PCRE2 _{^{v10.34 or later}} / .NET), 25 22 bytes

((\1|^)(\3xx|^x)*){3}$

Takes its input in unary, as a string of x characters whose length represents the number.

Try it online! / Attempt This Online! - Perl
Attempt This Online! - PCRE2
Try it online! - .NET

This is a modification of my Sum of two squares answer. Checking for 3 squares instead of 2 required using (\1|^) instead of \1?, dropping Java compatibility (apparently due to bugs in Java's regex engine). This is only a +2 byte change thanks to it allowing the initial ^ to be removed.

This can be generalized to sums of any number of squares by changing the 3 in {3} to the desired count. (With any count ≥4, it will always return true.)

Commented and indented:

                   # tail = N = input number
(                  # The following will be captured in \1
    (\1|^)         # If not on the first iteration, subtract \1, which is the
                   # previous perfect square.
    (\3xx|^x)*     # If on the first iteration, subtract a perfect square from
                   # tail; if on a later iteration and "(\1|^)" already
                   # subtracted \1, combine with it to subtract a perfect square
                   # from tail that is greater than or equal to the previous one.
){3}               # Iterate this three times
$                  # Assert that tail == 0

Retina 0.8.2, 31 28 bytes

.+
$*
((\1|^)(11\3|^1)*){3}$

Try it online!

This adapts the Perl / PCRE2 / .NET version from above.

Regex (ECMAScript or better), 45 31 bytes

^((x+)\2\2(?=\2$))*x{7}(x{8})*$

Try it online! - ECMAScript
Try it online! - Perl
Try it online! - Java
Try it online! - Python
Try it online! - Ruby
Try it online! - PCRE
Try it online! - .NET

Uses Legendre's three-square theorem. Outputs truthy if the number can be represented as \${4^a}(8b+7)\$, meaning it outputs falsey if the input number can be represented as the sum of 3 squares.

The algorithm used to repeatedly divide by 4 is related to the third form of Powers of 2, which repeatedly divides by 2.

^                 # tail = N = input number
(
    (x+)\2\2(?=\2$)   # Assert tail is divisible by 4;
                      # tail /= 4
)*                # Iterate the above any number of times, min zero
x{7}              # tail -= 7
(x{8})*$          # Assert tail is divisible by 8

Regex (Perl / PCRE / Java / .NET), 31 bytes

^((\2{4}|^xxx)*x)\1{6}(\1{8})*$

Try it online! / Attempt This Online - Perl
Try it online! / Attempt This Online - PCRE
Try it online! / Attempt This Online - Java
Try it online! - .NET

Uses Legendre's three-square theorem. Outputs truthy if the number can be represented as \${4^a}(8b+7)\$, meaning it outputs falsey if the input number can be represented as the sum of 3 squares.

The algorithm used to match powers of 4 is the same as the first form of Powers of 2, which for powers of 2 is ^(\1\1|^x)*x$.

^                  # tail = N = input number
((\2{4}|^xxx)*x)   # \1 = any arbitrary power of 4; tail -= \1
\1{6}              # tail -= \1 * 6
(\1{8})*           # tail -= \1 * 8 * {any arbitrary positive number}
$                  # Assert tail == 0

Regex (ECMAScript or better), 37 bytes

^((x(x*))(?=.*(?=\2*$)(\3+$))\4|){3}$

Try it online! - ECMAScript
Try it online! - Perl
Try it online! - Java
Try it online! - Python
Try it online! - Ruby
Try it online! - PCRE
Try it online! - .NET

Uses a variant of the multiplication/squaring algorithm explained in this post.

This is a trivial modification of my Sum of two squares answer, changing the {2} to {3}. It can be generalized to sums of any number of squares by changing that to the desired count. (With any count ≥4, it will always return true.)

^                  # tail = N = input number
(                  # subtract a perfect square from tail
    # subtract a nonzero perfect square from tail
    (x(x*))        # \2 = any positive integer; \3 = \2-1; tail -= \2
    (?=
        .*         # find the smallest value of tail that satisfies the following
        (?=\2*$)   # assert \2 divides tail
        (\3+$)     # assert \3 divides tail at least once; \4 = \2*\2 - \2
    )
    \4             # along with "tail -= \2" above, acts as "tail -= \2*\2"
|
    # or just subtract 0 from tail
){3}               # do this three times
$                  # assert that tail == 0

K (ngn/k), 12 bytes

|/7=8!%[;4]\

Try it online!

Repeatedly divides by 4 until the number becomes 0, collects all the intermediate values, and tests if any of them is 7 modulo 8.

This does not lead to infinite loop because there exists the smallest nonzero positive number that can be represented in float. This is theoretically correct if we assume arbitrary (but not infinite) precision floats, because the only moment the number can be 7 modulo 8 is right before it has nonzero fractional part.

K (ngn/k), 19 bytes

3>#&~>/*:\-/&+2\|4\

Try it online!

It's much more mind-boggling than I thought to test if a number ends with three ones followed by an even number of zeros in binary.

|4\    convert to base 4, and then reverse to put least significant digit first
+2\    convert each digit to base 2, making sure that each digit becomes two bits
       (actually it isn't correct for 0 or numbers that don't have 2 or higher in
       base 4, but they are treated as original number times 2, which doesn't change
       the answer (truthy))
-/&    take the 2d positions of one bits and calculate (row index - column index)
       high bit becomes (row index) and low bit becomes (row index - 1)
~>/*:\ mark each number that is less than or equal to the first element as 1
       if the first 1 is a high bit, there are at most 3 places that this can happen:
       itself, its low bit counterpart, and the next higher digit's low bit;
       if all three are present, it is precisely 4^k*(8m+7)
       otherwise, only itself can meet the condition
3>#&   test if there are less than three ones

Japt -x¡, 13 11 bytes

ô²ïÕà3 £¶Xx

Try it

ô²ïÕà3 £¶Xx     :Implicit input of integer U
ô               :Range [0,U]
 ²              :Square each
  ï             :Cartesian product
   Õ            :Reduce each pair by GCD
    à3          :Combinations of length 3
       £        :Map each X
        ¶       :  Is U equal to
         Xx     :    Sum of X
                :Implicit output of sum of resulting array as a Boolean

Charcoal, 18 bytes

≔⮌⍘Ｎ⁴θ⬤θ⊖↔⁻³²Ｉ…θ⊕κ

Try it online! Link is to verbose version of code. Explanation: Uses the theorem from @xnor's Python answer.

≔⮌⍘Ｎ⁴θ⬤θ

Convert the input to base 4 and reverse it so that impossible values start with ...00031 or ...00033.

⊖↔⁻³²Ｉ…θ⊕κ

Check that no prefix of the base 4 string has a decimal value that differs by 1 from 32.

Husk, 8 bytes

€mṁ□π3ṡ¹

Try it online! (outputs first 20 truthy integers)

      ṡ¹ # symmetric range of input:
         #  -input ... +input
    π3   # cartesian power of 3
 m       # map over this:
  ṁ□     #  sums of squares 
€        # index of input, or zero if absent

Retina 0.8.2, 39 bytes

.+
$*
^((^1|11\2)*\1?){2}(1(?(3)1\3))*$

Try it online! Link includes simple test cases. Explanation: My original Retina 0.8.2 answer to Sum of two squares was readily extended to a third term but I then applied @Deadcode's golf for the first two terms.

.+
$*

Convert to unary. (Very inefficient for large numbers of course.)

^((^1|11\2)*\1?){2}

Match a square number, but optionally also match it again or a previous square number. (Unfortunately this trick doesn't work for a third square number.)

(1(?(3)1\3))*$

Match a third square number.

A port of @xnor's Python answer is only 37 bytes (positive result) or 35 bytes (inverted result).

.+
$*
+`^(.)\1{3}$
$1
^(.{8})*.{0,6}$

Try it online! Invert the result by replacing 0,6 with 7. Explanation:

.+
$*

Convert to unary.

+`^(.)\1{3}$
$1

Divide out any powers of 4.

^(.{8})*.{0,6}$

Check the remainder after dividing by 8.

K (ngn/k), 26 24 bytes

{:[4!x|~x;7>8!x;o -4!x]}

Try it online!

-2 bytes thanks to Bubbler

Return 1 for true and 0 for false. Direct implementation of l4m2's JS solution. Definitely could be golfed further, considering the approach.

Explanation:

{:[4!x|~x;7>8!x;o -4!x]}    Main function. x is input.
 :[                   ]     If
       ~x                   Not x (Check if x is 0)
      |                     Or
   4!x                      x is divisible by 4 (x mod 4)
         ;  8!x             Then, x mod 8...
          7>                Smaller than 7?
               ;o -4!x      Else, recursively call function with argument x div 4

Uiua, 20 chars

/↥=0⍥₃(♭⬚η∵(-ⁿ₂⇡⌊.))

Pad Link

Nekomata + -e, 5 bytes

Ṗ√4~L

Attempt This Online!

Times out for large test cases.

Ṗ√4~L
Ṗ       Find an integer partition of the input
 √      Square root, fail if not perfect square
    L   Check if the length is
  4~    0, 1, 2, or 3

Nekomata + -e, 9 bytes

Zʷ{4¦}→8¦

Attempt This Online!

Using Legendre's three-square theorem like @xnor's answer.

This swaps truthy and falsy.

Zʷ{4¦}→8¦
Z           Check that the input is nonzero
 ʷ{  }      Loop until failure
   4¦           Divide by 4, fail if not divisible
      →     Increment
       8¦   Check if divisible by 8

Pyth, 10 bytes

/sM^^R2hQ3

Attempt This Online!

/sM^^R2hQ3Q­⁡​‎‎⁡⁠⁢⁡‏⁠‎⁡⁠⁢⁢‏⁠‎⁡⁠⁢⁣‏⁠‎⁡⁠⁢⁤‏⁠‎⁡⁠⁣⁡‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁤‏⁠‎⁡⁠⁣⁢‏‏​⁡⁠⁡‌⁣​‎‎⁡⁠⁢‏⁠‎⁡⁠⁣‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁡‏⁠‎⁡⁠⁣⁣‏‏​⁡⁠⁡‌­
    ^R2hQ    # ‎⁡Generate the first input square numbers
   ^     3   # ‎⁢3d cartesian product with generated squares
 sM          # ‎⁣Reduce each element on sum
/         Q  # ‎⁤Count occurrences of input 
💎

Created with the help of Luminespire.

JavaScript (Node.js), 24 bytes

f=x=>!x|x%4?x%8<7:f(x/4)

Try it online!

Theorem xnor used