# Count sums of two squares

Given a non-negative number n, output the number of ways to express n as the sum of two squares of integers n == a^2 + b^2 (OEIS A004018). Note that a and b can be positive, negative, or zero, and their order matters. Fewest bytes wins.

For example, n=25 gives 12 because 25 can be expressed as

(5)^2  + (0)^2
(4)^2  + (3)^2
(3)^2  + (4)^2
(0)^2  + (5)^2
(-3)^2 + (4)^2
(-4)^2 + (3)^2
(-5)^2 + (0)^2
(-4)^2 + (-3)^2
(-3)^2 + (-4)^2
(0)^2  + (-5)^2
(3)^2  + (-4)^2
(4)^2  + (-3)^2


Here are the values up to n=25. Be careful that your code works for n=0.

0 1
1 4
2 4
3 0
4 4
5 8
6 0
7 0
8 4
9 4
10 8
11 0
12 0
13 8
14 0
15 0
16 4
17 8
18 4
19 0
20 8
21 0
22 0
23 0
24 0
25 12


Here are the values up to n=100 as a list.

[1, 4, 4, 0, 4, 8, 0, 0, 4, 4, 8, 0, 0, 8, 0, 0, 4, 8, 4, 0, 8, 0, 0, 0, 0, 12, 8, 0, 0, 8, 0, 0, 4, 0, 8, 0, 4, 8, 0, 0, 8, 8, 0, 0, 0, 8, 0, 0, 0, 4, 12, 0, 8, 8, 0, 0, 0, 0, 8, 0, 0, 8, 0, 0, 4, 16, 0, 0, 8, 0, 0, 0, 4, 8, 8, 0, 0, 0, 0, 0, 8, 4, 8, 0, 0, 16, 0, 0, 0, 8, 8, 0, 0, 0, 0, 0, 0, 8, 4, 0, 12]


Fun facts: The sequence contains terms that are arbitrarily high, and the limit of its running average is π.

    var QUESTION_ID=64812,OVERRIDE_USER=20260;function answersUrl(e){return"https://api.stackexchange.com/2.2/questions/64812/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function commentUrl(e,s){return"https://api.stackexchange.com/2.2/answers/"+s.join(";")+"/comments?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+COMMENT_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(answer_page++),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){answers.push.apply(answers,e.items),answers_hash=[],answer_ids=[],e.items.forEach(function(e){e.comments=[];var s=+e.share_link.match(/\d+/);answer_ids.push(s),answers_hash[s]=e}),e.has_more||(more_answers=!1),comment_page=1,getComments()}})}function getComments(){jQuery.ajax({url:commentUrl(comment_page++,answer_ids),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){e.items.forEach(function(e){e.owner.user_id===OVERRIDE_USER&&answers_hash[e.post_id].comments.push(e)}),e.has_more?getComments():more_answers?getAnswers():process()}})}function getAuthorName(e){return e.owner.display_name}function process(){var e=[];answers.forEach(function(s){var r=s.body;s.comments.forEach(function(e){OVERRIDE_REG.test(e.body)&&(r="<h1>"+e.body.replace(OVERRIDE_REG,"")+"</h1>")});var a=r.match(SCORE_REG);a&&e.push({user:getAuthorName(s),size:+a[2],language:a[1],link:s.share_link})}),e.sort(function(e,s){var r=e.size,a=s.size;return r-a});var s={},r=1,a=null,n=1;e.forEach(function(e){e.size!=a&&(n=r),a=e.size,++r;var t=jQuery("#answer-template").html();t=t.replace("{{PLACE}}",n+".").replace("{{NAME}}",e.user).replace("{{LANGUAGE}}",e.language).replace("{{SIZE}}",e.size).replace("{{LINK}}",e.link),t=jQuery(t),jQuery("#answers").append(t);var o=e.language;/<a/.test(o)&&(o=jQuery(o).text()),s[o]=s[o]||{lang:e.language,user:e.user,size:e.size,link:e.link}});var t=[];for(var o in s)s.hasOwnProperty(o)&&t.push(s[o]);t.sort(function(e,s){return e.lang>s.lang?1:e.lang<s.lang?-1:0});for(var c=0;c<t.length;++c){var i=jQuery("#language-template").html(),o=t[c];i=i.replace("{{LANGUAGE}}",o.lang).replace("{{NAME}}",o.user).replace("{{SIZE}}",o.size).replace("{{LINK}}",o.link),i=jQuery(i),jQuery("#languages").append(i)}}var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe",COMMENT_FILTER="!)Q2B_A2kjfAiU78X(md6BoYk",answers=[],answers_hash,answer_ids,answer_page=1,more_answers=!0,comment_page;getAnswers();var SCORE_REG=/<h\d>\s*([^\n,]*[^\s,]),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/,OVERRIDE_REG=/^Override\s*header:\s*/i;
    body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px}
    <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"> <div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr></thead> <tbody id="answers"> </tbody> </table> </div><div id="language-list"> <h2>Winners by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr></thead> <tbody id="languages"> </tbody> </table> </div><table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody> </table> <table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody> </table>

• Wait, what?? "The sequence contains terms that are arbitrarily high, and the limit of its running average is π." Nov 26, 2015 at 11:58
• @StewieGriffin The two statements are consistent. Consider the sequence 1,0,2,0,0,3,0,0,0,4,0,0,0,0,5,.... Cutting the sequence off after any nonzero number, the average so far is 1. And, the runs of 0's have less and less impact later in the sequence.
– xnor
Nov 26, 2015 at 12:04
• I know it's consistent.. =) I had checked the 10.000 first numbers when I posted the comment. What I don't get is: Why on earth does it equal Pi? Nov 26, 2015 at 12:19
• @StewieGriffin The sum of the terms up to N corresponds to the points (a,b) with a^2+b^2<=N. These are the lattice points in the circle of radius sqrt(N), whose area is πN.
– xnor
Nov 26, 2015 at 12:23
• @xnor and there goes the magic:( Nov 29, 2015 at 1:02

## Python (59 57 56 bytes)

lambda n:0**n+sum((-(n%(x-~x)<1))**x*4for x in range(n))


Online demo

As with my CJam answer, this uses Möbius inversion and runs in pseudoquasilinear time.

Thanks to Sp3000 for 2 bytes' savings, and feersum for 1.

• Those parentheses are annoying. Nov 30, 2015 at 22:15
• @ThomasKwa, tell me about it. The thing which really surprised me, in one of the versions I passed through on the way to the first one I posted, was that -1**x is always -1. I expected the -1 to be a single integer literal token rather than a low-precedence unary minus followed by a one. Nov 30, 2015 at 22:21
• Congrats on the bounty! Your code is shorter than anything I had come up with. Your solution is based on a totally new and unexpected mathematical idea, and it brings me joy to that this is possible even in such a straightforward-looking challenge .The Mobius-inverse result is quite pretty and I took the pleasure of deriving a proof for myself.
– xnor
Dec 5, 2015 at 8:45
• 1 byte can be saved by moving the multiplication by 4 to after **x. Dec 5, 2015 at 9:05
• @PeterTaylor Can you elaborate how your algorithm works/can you point me to a resource? I cannot quite see how you can apply the möbius inversion to number of sums of suqares problem. Dec 5, 2015 at 12:55

# Python 2, 44 bytes

f=lambda n,x=1:+(x>n)or(n%x<1)-f(n,x+2)/4<<2


This is almost the same as xsot's solution (which is based on Peter Taylor's solution), but saves 8 bytes by simplifying the way signs are handled.

Note that for a full program, we can save 2 bytes in the function without incurring a cost outside the function:

f=lambda n,x=1:x>n or(n%x<1)-f(n,x+2)/4<<2
print+f(input())


Two additional bytes for a full program this way:

n=input()
f=lambda x:x>n or(n%x<1)-f(x+2)/4<<2
print+f(1)


For n > 0 there is a very legible 40-byte solution:

f=lambda n,x=1:n/x and(n%x<1)*4-f(n,x+2)

• Congrats for winning the bounty! The recursive subtraction is a clean and short way to express the alternating sum for odd divisors without needing to extract a sign from the divisor itself. Also, credit to xsot for streamlining Peter Taylor's solution to a recursive one with clever handling of n=0.
– xnor
Dec 13, 2015 at 7:33

# Mathematica, 13 bytes

If built-ins are allowed, this is how to do it in Mathematica.

2~SquaresR~#&


For 0 < = n <= 100

2~SquaresR~# & /@ Range[0, 100]


{1, 4, 4, 0, 4, 8, 0, 0, 4, 4, 8, 0, 0, 8, 0, 0, 4, 8, 4, 0, 8, 0, 0, 0, 0, 12, 8, 0, 0, 8, 0, 0, 4, 0, 8, 0, 4, 8, 0, 0, 8, 8, 0, 0, 0, 8, 0, 0, 0, 4, 12, 0, 8, 8, 0, 0, 0, 0, 8, 0, 0, 8, 0, 0, 4, 16, 0, 0, 8, 0, 0, 0, 4, 8, 8, 0, 0, 0, 0, 0, 8, 4, 8, 0, 0, 16, 0, 0, 0, 8, 8, 0, 0, 0, 0, 0, 0, 8, 4, 0, 12}

• Because of course Mathematica has a built-in for this. Dec 5, 2016 at 15:01

# Pyth, 13 bytes

/sM^^R2}_QQ2Q


Test suite

/sM^^R2}_QQ2Q
Q = eval(input())
}_QQ      Inclusive range from -Q to Q (all possible a and b)
^R2          Map to their squares
^       2     Form all pairs
sM              Sum pairs
/           Q    Count occurances of Q

• Late, but I don't think you need the last Q. Dec 5, 2016 at 13:07

## J, 16 bytes

+/@,@:=+&*:/~@i:


This is a monadic verb (in other words, a unary function). Try it online or see it pass all test cases.

## Explanation

+/@,@:=+&*:/~@i:  Denote input by n
i:  The array of integers from -n to n
/~@    Take outer product wrt the following function:
+           the sum of
&*:        squares of both inputs
This results in a 2D array of a^2+b^2 for all a, b between -n and n
=           Replace by 1 those entries that are equal to n, and others by 0
,@:            Flatten the binary matrix
+/@               Take its sum


# Python 2, 695553 52 bytes

f=lambda n,x=1:+(x>n)or(2-x%4)*(n%x<1)+f(n,x+2)/4<<2


This is a recursive function based off Peter Taylor's excellent solution.

• This is a great improvement. But, there is still a way to make it shorter, and I encourage you to look for it.
– xnor
Dec 8, 2015 at 22:29
• @xnor Another byte down. I hope you don't have any more tricks up your sleeves.
– xsot
Dec 9, 2015 at 11:32
• I don't know whether I should make an answer of it, it's just your solution plus one trick: f=lambda n,x=1:+(x>n)or(n%x<1)-f(n,x+2)/4<<2. Also, I guess we don't care about exceeding the default maximum recursion depth? Dec 9, 2015 at 20:38
• @MitchSchwartz I think that's an incredible improvement worthy of the bounty and likely the final optimisation xnor had in mind.
– xsot
Dec 9, 2015 at 23:07
• @MitchSchwartz Yes, that's the optimization I was thinking of! And xsot's /4<<2 trick makes it shorter than what I had.
– xnor
Dec 9, 2015 at 23:30

# Julia, 40 bytes

n->n>0?4sum(i->(n-i^2)^.5%1==0,1:n^.5):1


Ungolfed:

function f(n)
if n==0
return 1           # Handle special case of n=0
else
m=0                # Start the counter at zero
for i=1:sqrt(n)    # Loop over the values (i) whose squares are
# less than n (except zero)
k=sqrt(n-i^2)    # Find k such that n=k^2+i^2
if k==floor(k)   # if k is an integer, we've found a pair
m+=4           # Add one for each of k^2+i^2, (-k)^2+(-i)^2, and the other two
end
end
return m           # Return the resulting count
end
end


Note that the loop doesn't include i==0, because when n is a square, it's already included by i=sqrt(n), and there are only four, not eight, for that form (0^2+k^2, 0^2+(-k)^2, k^2+0^2, (-k)^2+0^2).

# CJam, 25 23 bytes

Zri:R#Ym*{Rf-Yf#:+R=},,


This is a theoretical solution that requires O(9n) time and memory for input n.

At the cost of one extra byte – for a total of 24 bytes – we can reduce the complexity to O(n2):

ri:R)Y*Ym*{Rf-Yf#:+R=},,


Try it online!

### How it works

Either

Z                  Push 3.
ri:R              Read an integer from STDIN and save it in R.
#             Compute 3**R.


or

ri:R               Read an integer from STDIN and save it in R.
)Y*            Add 1 and multiply by 2.


Then

Ym*                Take the second Cartesian power, i.e., compute all pairs.
{          },   Filter the pairs:
Rf-              Subtract R from each.
Yf#           Square the differences.
R=       Compare with R.
If = pushed 1, keep the pair.
,  Count the kept pairs.

• And at a saving of one byte it's possible to get the complexity down to Õ(n) Nov 30, 2015 at 14:18
• Yes, I saw. That's amazing. Nov 30, 2015 at 14:20

## CJam (25 24 22 21 bytes)

{:X!X{X\2*)%!4*\-}/z}


Online demo

This runs in pseudoquasilinear time* and uses the statement from the OEIS that

Moebius transform is period 4 sequence [ 4, 0, -4, 0, ...]. - Michael Somos, Sep 17 2007

The input 0 is obviously a special case (Möbius tranforms and annihilators don't go well together), but ended up only costing one char.

* Pseudo- because it's quasilinear in the value of the input, not the size of the input; quasi because it does Theta(n) operations on integers of size on the order of n; a b-bit mod operation should take b lg b time, so overall this takes Theta(n lg n lg lg n) time.

f n|q<-[-n..n]=sum[1|a<-q,b<-q,a*a+b*b==n]


Usage exapmle:

*Main> map f [0..25]
[1,4,4,0,4,8,0,0,4,4,8,0,0,8,0,0,4,8,4,0,8,0,0,0,0,12]
*Main>

• Binding q in a guard is clever, I'll remember this trick.
– xnor
Nov 26, 2015 at 11:26

# Japt, 4237 33 bytes

Japt is a shortened version of JavaScript. Interpreter

V=Un oU+1;Vr@X+(Vf_p2 +Y*Y¥U)l ,0


### How it works

           // Implicit: U = input number
V=Un oU+1  // Set variable V to range(-U, U+1). Ends up like [-U,-U+1,...,U-1,U]
Vr@    ,0  // Reduce each item Y in V with this function, starting at 0:
X+(     l  //  Return the previous value X + the length of:
Vf_p2      //   V filtered by items Z where Z*Z
+Y*Y==U)   //    + Y*Y equals U.
// This ends up as the combined length of all fitting pairs of squares.
// Implicit: return last expression


Perhaps there's a better technique; suggestions are welcome.

## Python 3, 6861 60 bytes

lambda n:0**n+4*sum(i**.5%1+(n-i)**.5%1==0for i in range(n))


Using two nested list comprehensions is too expensive. Instead, this checks if both coordinates on the circle of radius sqrt(n) are integers.

Peter Taylor has beaten this with a Moebius-inversion based approach.

• Well done. I was tinkering with a recursive function but couldn't resolve n=0 elegantly.
– xsot
Nov 30, 2015 at 1:48

# Octave, 28 bytes

@(n)nnz((a=(-n:n).^2)'+a==n)


# Julia, 8979 63 bytes

g(x)=cos(π*x^.5)^2÷1
a(n)=(n==0)+4sum([g(i)g(n-i)for i=1:n])


This is a named function a which accepts an integer and returns a float. It calls a helper function g.

Ungolfed:

function g(x::Integer)
floor(cos(π*sqrt(x))^2)
end

function a(n::Integer)
(n == 0) + 4*sum([g(i)*g(n-i) for i=1:n])
end


The approach here uses a simplification of Wesley Ivan Hurt's formula listed on OEIS. The simplification was found by Glen O, the very same person who shaved 26 bytes off of this answer!

• Use x^.5 rather than sqrt(x) to save 3 bytes. And (n==0) saves 2 bytes over 1÷(n+1). And you can save 4 more characters by using cos(π*sqrt(x))^2÷1 rather than floor(cos(π*sqrt(x))^2). Also, use 1:n/2 rather than 1:n÷2, because there's no harm using a float in g(x) and it'll be locked to the integers for i anyway. And sum(i->g(i)g(n-i),1:n/2) will shave some more characters, too. Nov 26, 2015 at 3:12
• @GlenO Great suggestions, thanks. The sum trick fails for n=0 though, so I kept array comprehension. Nov 26, 2015 at 3:45
• So, it can be recovered - if you let the i=0 case happen in the sum, you can switch the sign on 4g(n). So (n==0)-4g(n)-4g(n/2)+8sum(i->g(i)g(n-i),0:n/2), which will not run into the error. But you can do even better, by noting the symmetries - (n==0)+4sum([g(i)g(n-i)for i=1:n]) Nov 27, 2015 at 3:11
• @GlenO That simplification is seriously genius. I recommend you submit that as an alternative formula for the sequence on OEIS! Nov 27, 2015 at 18:44

# Pyth, 16 15 bytes

lfqQs^R2T^}_QQ2


Try it online in the Pyth Compiler.

### How it works

lfqQs^R2T^}_QQ2

}_QQ   Compute the inclusive range from -Q to Q (input).
^    2  Take the second Cartesian power, i.e., compute all pairs.
f               Filter; for each T in the list of pairs:
^R2T          Compute the squares of T's elements.
qQ               Compare the sum with Q.
If q returned True, keep T.
l                Count the kept pairs.


## TI-BASIC, 23 bytes

sum(seq(Σ(X²+Y²=Ans,X,-Ans,Ans),Y,-Ans,Ans


Pretty straightforward. Σ( is summation.

Strangely, sum(seq(sum(seq( throws an ERR:ILLEGAL NEST, and so does Σ(Σ(, but sum(seq(Σ( is fine. I chose to put the Σ( on the inside to save a close-paren.

• What's the difference between sum and Σ? Nov 26, 2015 at 9:28
• @alephalpha Σ( takes a summation, adding up all of the X²+Y²=Ans from values of X between -Ans and Ans. sum( is sum of a list, so we need to create the list first using seq(...,Y,-Ans,Ans Nov 26, 2015 at 18:11

# JavaScript (ES6), 66 60 bytes

n=>eval("for(r=0,a=~n;a++<n;)for(b=~n;b++<n;)r+=a*a+b*b==n")


6 bytes saved thanks to @edc65!

## Explanation

n=>eval(              // use eval to allow for loops in an unparenthesised arrow function
for(r=0,             // r = number of pairs
a=~n;a++<n;        // a = first number to square
)
for(b=~n;b++<n;) // b = second number to square
r+=a*a+b*b==n  // add one to the result if a^2 + b^2 == n
// implicit: return r
)


## Test

n = <input type="number" oninput='result.innerHTML=(

n=>eval("for(r=0,a=~n;a++<n;)for(b=~n;b++<n;)r+=a*a+b*b==n")

)(+this.value)' /><pre id="result"></pre>

• 60: n=>eval('for(r=0,a=~n;a++<n;)for(b=~n;b++<n;)r+=a*a+b*b==n') Nov 26, 2015 at 12:52
• @edc65 Nice! I didn't think of using eval to put the for loops into an arrow function without parentheses. I also forgot about the ~ operator haha. Nov 26, 2015 at 14:20

# Python 3, 9362 69 bytes

Itertools wasn't working so I used two ranges again, but moved the range out to save bytes.

Edit: Previous code didn't actually work, since I defined range over n before I defined n.

lambda n:sum(i*i+j*j==n for i in range(-n,n+1)for j in range(-n,n+1))


## APL, 2320 19 bytes

{+/⍵=∊∘.+⍨×⍨0,,⍨⍳⍵}


Explanation:

{+/⍵=∊∘.+⍨×⍨0,,⍨⍳⍵}        Monadic function:
⍳⍵          1 2 3 ... ⍵
,⍨            Duplicate
0,              Concatenate to 0
×⍨                 Square everything
∊                       Flatten
⍵=                        1s where equal to the input
+/                          Sum up the 1s


Other than the fact that APL doesn't have J's i: (numbers from -n to n) function, this works pretty much like the J answer.

We can't use a train because getting the -\⍳2×⍵ to not parse as (-\) ⍳ (2×⍵) would cost three bytes; similarly with other pair of atops. All those parentheses make the regular function shorter.

Try it here. The output of 1 means all values match.

# Matlab 41 bytes

Even smaller as the previous answers

@(n)nnz(~mod(sqrt(n-(1:n^.5).^2),1))*4+~n


Essentially Agawa001's answer with power and sqrt replaced

# Candy, 17 14 bytes

Input pushed onto stack initially

~TbAT1C(sWs+Aeh)Z

~T0C(sWs+Aeh)Z

peekA    # copy arg from stack to register A
range2   # create double sided range on stack, -A, 1-A, ... A-1, A
digit0   # prefix argument to 'cart',
cart     # cartesian product of current stack(0), and other stack(0)
while    # while stack not empty
sqr    # pop and square and push
swap   # swap two stack elements
sqr    # pop and square and push
pushA  # push original argument
equal  # equality test 0/1
popAddZ  # Z := Z + pop
endwhile
pushZ    # push Z onto stack, will be output to stdout on termination


# CJam, 28

qi_mF{3a>},{~)\4%2-&}%4+:*1?


Not really short, but efficient. E.g. the result for 15625 is instantly 28. Uses the factorization-based formula from OEIS.
Try it online

Explanation:

qi       read input and convert to integer
_        make a copy (will be used to handle the 0 case at the end)
mF       factorize into [prime exponent] pairs
{…},     filter the array of pairs
3a>    with the condition that the pair is greater than [3]
which means the prime factor must be ⩾3
{…}%     transform each pair as follows:
~      dump the prime factor and exponent onto the stack
)      increment the exponent
\      swap with the prime
4%     get the remainder mod 4 (it will be 1 or 3)
2-     subtract 2 (resulting in -1 or 1)
&      bitwise AND with the incremented exponent (see below)
4+       append a 4 to the array
:*       multiply all
1?       if the input was 0, use 1, else use the above result


• if the prime is 1 mod 4, the code calculates (exponent + 1) & -1, which is exponent + 1
• if the prime is 3 mod 4, the code calculates (exponent + 1) & 1, which is 0 if the exponent is odd, and 1 if even

All these values multiplied together and multiplied by 4 are exactly the OEIS formula.

# Python 2, 68 bytes

def x(n):r=range(-n,n+1);print sum(a*a+b*b==n for a in r for b in r)


Defines a function called x() that takes a number n.

Try it online. http://ideone.com/aRoxGF

• You are missing a print or return statement. Nov 28, 2015 at 23:47
• Thanks, I completely forgot. The link has the print statement though. I edited my code while I was making the code. Nov 28, 2015 at 23:51
• Ok, no worries. But this also seems to give incorrect results for n=0 and n=1 (0 and 2 instead of 1 and 4). Maybe the range limits need adjusting? Nov 28, 2015 at 23:56
• @Zgarb Yeah, they should end at n+1. Nov 28, 2015 at 23:56
• I'll look for it. Nov 29, 2015 at 1:23

# Pyth, 41353330 27 bytes

Try it online.

Edit: Thanks to isaacg, I got m and *F to work! YES!

?Q*F+4m.&tt%ed4hhdr-PQ2 8 1
(implicit) Q = input()
?Q                              If Q != 0
m                           Map to d (exponent, prime) from ...
r-PQ2 8         run-length-encoded(PQ with 2's removed)
.&                           Bitwise and
%ed4                       d[-1] % 4
tt                           -2
hd                  with d[0]
h                      +1
+4                            Append 4 to the resulting array
*F                              Then multiply it all together
1     Else 1


Edit: Put the bitwise and back in for more byte savings! Also I removed all the "Formerly" stuff. It was starting to get cluttered.

Thanks to aditsu and his CJam solution, and to maltysen and his tips (One day I will get m*Fd to work. One day...)

J4Vr-PQ2 8=J*J.&tt%eN4hhN;?QJ1
(implicit) Q=input()
J4                              J=4
-PQ2                        Remove all 2's from the prime factorization of Q
r     8                      run-length encode (exponent, prime factor)
V                      ;      For N in range( the above ):
=J*J                      J = J * ...
tt%eN4                N[-1] mod 4 -2
hhN             (exponent + 1)
.&                    bitwise and
?QJ1  if Q != 0 print(J) else print(1)


Note that,

• if the prime is 1 mod 4, we get -1 & (exponent + 1), which is exponent + 1

• but if the prime is 3 mod 4, we get 1 & (exponent + 1), which is 0 if the exponent is odd, and 1 if even

Multiply it all together (times 4 at the beginning) and we get the number of sums of two squares that add up to our input.

# Python, 57 bytes

Nice challenge. Unfortunately I'm not getting it shorter than this at the moment.

lambda n:0**n+sum(2-d%4for d in range(1,n+1)if d%2>n%d)*4


# PARI/GP, 34 28 bytes

Using generating functions:

Saved 6 bytes thanks to Mitch Schwartz.

n->sum(i=-n,n,x^i^2)^2\x^n%x


Using built-ins, 33 bytes (saved 1 byte thanks to Mitch Schwartz.):

n->if(n,2*qfrep(matid(2),n)[n],1)


qfrep(q,B,{flag=0}): vector of (half) the number of vectors of norms from 1 to B for the integral and definite quadratic form q. If flag is 1, count vectors of even norm from 1 to 2B.

• matid(2) saves a byte. Dec 9, 2015 at 22:39
• And down to 28 for the generating function approach: n->sum(i=-n,n,x^i^2)^2\x^n%x Dec 9, 2015 at 22:48

## Matlab, 72 bytes

n=input('');m=fix(sqrt(n));m=(-m:m).^2;disp(nnz(bsxfun(@plus,m,m')==n))

• You don't need disp in this challenge Luis! =) Nov 26, 2015 at 12:21
• @StewieGriffin Thanks! But in this case it's a program, not a function. So according to the accepted answer in your link it's needed, isn't it? Nov 26, 2015 at 13:06

# Matlab, 63 50 bytes

@(y)nnz(~mod(sqrt(y-power((1:sqrt(y)),2)),1))*4+~y


• This does beat the other same-entitled code, thus I put it :D.

• The program finds the positive integer solutions then multiply by 4 to encompass negative ones.

• It can perform all 25 first test cases

for i=1:25 ans(i)
end

1

4

4

0

4

8

0

0

4

4

8

0

0

8

0

0

4

8

4

0

8

0

0

0

0

12

• Nov 26, 2015 at 12:22
• thanks @StewieGriffin i did include it just as a fair-game relating to luis' one Nov 26, 2015 at 12:40
• Tips: When you're planning on posting the results from MATLAB, use format compact =) Nov 27, 2015 at 6:57

## Python, 133 129 bytes

from itertools import*
n=input()
l=[];s=0
for i in range(-n,n+1):l.append(i**2)
for j in product(l,l):
if sum(j)==n:s+=1
print s


This demonstrates the use of itertools cartesian product. Try it online

• @Mego I tried, I don't think there is anything to further golf. I will take that part out if you insist. Nov 28, 2015 at 21:01

# JavaScript, 89 bytes

n=prompt()
p=Math.pow
for (x=c=(+n?0:1);x<=n;x++)if(x&&p(n-p(x,2),.5)%1===0)c+=4