[Question]
  [
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](https://oeis.org/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 π.
**Leaderboard:**
```
    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>
```
  
[Answer]
## Python (59 57 56 bytes)
```
lambda n:0**n+sum((-(n%(x-~x)<1))**x*4for x in range(n))
```
[Online demo](http://ideone.com/nwd5fE)
As with my CJam answer, this uses Möbius inversion and runs in pseudoquasilinear time.
Thanks to [Sp3000](https://codegolf.stackexchange.com/users/21487/sp3000) for 2 bytes' savings, and [feersum](https://codegolf.stackexchange.com/users/30688/feersum) for 1.
[Answer]
# 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](https://codegolf.stackexchange.com/a/65077/2537) (which is based on Peter Taylor's [solution](https://codegolf.stackexchange.com/a/65326/2537)), 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)
```
[Answer]
# 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}
> 
> 
> 
[Answer]
# Pyth, 13 bytes
```
/sM^^R2}_QQ2Q
```
[Test suite](https://pyth.herokuapp.com/?code=%2FsM%5E%5ER2%7D_QQ2Q&test_suite=1&test_suite_input=0%0A1%0A2%0A3%0A4%0A5%0A6%0A7%0A8%0A9%0A10&debug=0)
```
/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
```
[Answer]
## J, 16 bytes
```
+/@,@:=+&*:/~@i:
```
This is a monadic verb (in other words, a unary function).
[Try it online](https://tio.run/nexus/j#@5@mYGuloK3voONgZautpmWlX@eQafU/NTkjXyFNwcj0PwA "J – TIO Nexus") or [see it pass all test cases](https://tio.run/nexus/j#@5@mYGuloK3voONgZautpmWlX@eQafU/NTkjX0FDR08hTclAUyFTT8HI7D8A "J – TIO Nexus").
## 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
```
[Answer]
# Python 2, ~~69~~ ~~55~~ ~~53~~ 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](https://codegolf.stackexchange.com/a/65326/45268).
[Answer]
# 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`).
[Answer]
# 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!](http://cjam.tryitonline.net/#code=cmk6UikyKlltKntSZi1ZZiM6K1I9fSws&input=MjU)
### 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.
          :+         Add the squares.
            R=       Compare with R.
                   If = pushed 1, keep the pair.
                ,  Count the kept pairs.
```
[Answer]
## CJam (25 24 22 21 bytes)
```
{:X!X{X\2*)%!4*\-}/z}
```
[Online demo](http://cjam.aditsu.net/#code=17%2C%0A%0A%7B%3AX!X%7BX%5C2*)%25!4*%5C-%7D%2Fz%7D%0A%0A%25%60)
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.
[Answer]
## Haskell, 42 bytes
```
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> 
```
[Answer]
# [Japt](https://github.com/ETHproductions/Japt), ~~42~~ ~~37~~ 33 bytes
**Japt** is a shortened version of **Ja**vaScri**pt**. [Interpreter](https://codegolf.stackexchange.com/a/62685/42545)
```
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.
[Answer]
## Python 3, ~~68~~ ~~61~~ 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](https://codegolf.stackexchange.com/a/65326/39328).*
[Answer]
# Octave, 28 bytes
```
@(n)nnz((a=(-n:n).^2)'+a==n)
```
[Answer]
# Julia, ~~89~~ ~~79~~ 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!
[Answer]
# Pyth, ~~16~~ 15 bytes
```
lfqQs^R2T^}_QQ2
```
Try it online in the [Pyth Compiler](https://pyth.herokuapp.com/?code=lfqQs%5ER2T%5E%7D_QQ2&input=25&debug=0).
### 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.
    s              Add the squares.
  qQ               Compare the sum with Q.
                 If q returned True, keep T.
l                Count the kept pairs.
```
[Answer]
## 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.
[Answer]
# 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](https://codegolf.stackexchange.com/users/21348/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>
```
[Answer]
# Python 3, ~~93~~ ~~62~~ 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))
```
[Answer]
## APL, ~~23~~ ~~20~~ 19 bytes
```
{+/⍵=∊∘.+⍨×⍨0,,⍨⍳⍵}
```
Explanation:
```
{+/⍵=∊∘.+⍨×⍨0,,⍨⍳⍵}        Monadic function:
                 ⍳⍵          1 2 3 ... ⍵
               ,⍨            Duplicate
             0,              Concatenate to 0
          ×⍨                 Square everything
      ∘.+⍨                   Make an addition table
     ∊                       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](http://tryapl.org/?a=%28%7B+/%u2375%3D%u220A%u2218.+%u2368%D7%u23680%2C%2C%u2368%u2373%u2375%7D%A80%2C%u2373100%29%u22611%204%204%200%204%208%200%200%204%204%208%200%200%208%200%200%204%208%204%200%208%200%200%200%200%2012%208%200%200%208%200%200%204%200%208%200%204%208%200%200%208%208%200%200%200%208%200%200%200%204%2012%200%208%208%200%200%200%200%208%200%200%208%200%200%204%2016%200%200%208%200%200%200%204%208%208%200%200%200%200%200%208%204%208%200%200%2016%200%200%200%208%208%200%200%200%200%200%200%208%204%200%2012&run). The output of `1` means all values match.
[Answer]
# 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
[Answer]
# [Candy](https://github.com/dale6john/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
  add    # pop and pop and add 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
```
[Answer]
# 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](http://cjam.aditsu.net/#code=qi_mF%7B3a%3E%7D%2C%7B~)%5C4%252-%26%7D%254%2B%3A*1%3F&input=15625)
**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
```
Some details about the calculation:
* 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.
[Answer]
# 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>
[Answer]
# Pyth, ~~41~~ ~~35~~ ~~33~~ ~~30~~ 27 bytes
[Try it online.](https://pyth.herokuapp.com/?code=%3FQ*F%2B4m.%26tt%25ed4hhdr-PQ2%208%201&test_suite=1&test_suite_input=0%0A1%0A2%0A3%0A4%0A5%0A6%0A7%0A8%0A9%0A10%0A21%0A65%0A362880&debug=1)
**Edit:** Thanks to [isaacg](http://chat.stackexchange.com/users/115644/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](https://codegolf.stackexchange.com/users/7416/aditsu) and his CJam solution, and to [maltysen](http://chat.stackexchange.com/users/138131/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.
[Answer]
# 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
```
[Answer]
# PARI/GP, ~~34~~ 28 bytes
Using generating functions:
Saved 6 bytes thanks to [Mitch Schwartz](https://codegolf.stackexchange.com/users/2537/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](https://codegolf.stackexchange.com/users/2537/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.
> 
> 
> 
---
[Answer]
## Matlab, 72 bytes
```
n=input('');m=fix(sqrt(n));m=(-m:m).^2;disp(nnz(bsxfun(@plus,m,m')==n))
```
[Answer]
# 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
```
[Answer]
## 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](http://ideone.com/iq0N9h)
[Answer]
# 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
alert(c)
```
I know this isn't the shortest JavaScript answer even if I were to remove the i/o lines, but I do think it is the best performing JS answer giving me the result for a million within a few seconds (ten million took about a minute).
]