Perl 5, 45 + 1 (-a) = 46 bytes

$_=!$_+4*grep{$_=sqrt"@F"-$_*$_;!/\./}1..sqrt

Try it online!

Japt -x, 11 bytes

õUn)ï £¶Xx²

Try it

õUn)ï £¶Xx²     :Implicit input of integer U
õ               :Inclusive range from U to
 Un             :  -U
   )            :End range
    ï           :Cartesian product with itself
      £         :Map each X
       ¶        :  Test U for equality with
        Xx      :    X reduced by addition after
          ²     :      Squaring each
                :Implicit output of sum of resulting array

R, 39 bytes

\(n)sum(rowSums(merge(-n:n,-n:n)^2)==n)

Attempt This Online!

I think this is a really compact solution for an R answer. Here, merge function creates the cartesian product of the range from the -n to n with itself (quite a generous estimation to avoid creating of new variables). Then each value is squared, values in each row are summed and the sum is compared to n. Thus the number of matches is produced.

Thunno 2 S, 7 bytes

ṄIDṖ²ʂ=

Try it online!

Explanation

ṄIDṖ²ʂ=  # Implicit input
ṄI       # Push [-input..input]
  DṖ     # Cartesian product
    ²ʂ   # Square and sum
      =  # Equals input?
         # Sum the list
         # Implicit output

Husk, 8 bytes

#¹´×+m□ṡ

Try it online! or Verify first 25 values

Husk, 8 bytes

#¹mṁ□π2ṡ

Try it online!

Jelly, 7 bytes

rNp`²§ċ

Try it online!

How it works

rNp`²§ċ - Main link. Takes n on the left
 N      - Yield -n
r       - Take the range [-n, -n+1, ..., -1, 0, 1, ..., n-1, n]
   `    - Use this list for both arguments for:
  p     -   Cartesian product
    ²   - Square each number
     §  - Take the sums of each pair
      ċ - Count the number of times n appears

05AB1E, 7 bytes

(ŸãnOQO

Try it online or verify the test cases in the range \$[0,100]\$.

Explanation:

(        # Get the negative of the (implicit) input-integer
 Ÿ       # Push a list in the range [(implicit) input-integer, -input]
  ã      # Get the cartesian product of this list, creating all possible pairs
   n     # Square each value in each pair
    O    # Sum each inner pair
     Q   # Check for each sum whether it's equal to the (implicit) input-integer
         # (1 if truthy; 0 if falsey)
      O  # And sum those
         # (after which the result is output implicitly)

MathGolf, 9 bytes

╤■mæ²Σk=Σ

Try it online.

Explanation:

╤         # Take the (implicit) input-integer, and push a list in the range [-input, input]
 ■        # Take the cartesian product of this, creating a list of all possible pairs
  mæ      # Map these pairs to, using the following four commands:
    ²     #  Take the square of both values in the pair
     Σ    #  Sum those
      k=  #  And check whether it's equal to the input-integer (1 if truthy; 0 if falsey)
  Σ       # After the map, sum the list
          # (after which the entire stack joined together is output implicitly)

Add++, 31 bytes

D,f,@,.5^1+iR2€Ω^d0BFB]d‽+A€=¦+

Try it online!

How it works

This defines a function, \$f\$, that takes the input, \$x\$, as an argument and returns the correct output. We start by yielding \$y := \lfloor\sqrt{x}+1\rfloor\$. We then push the range \$a := [1, 4, ..., y^2]\$, the list of square numbers up to the smallest square number larger than \$x\$. We then duplicate this array and push \$0\$ to the stack. At this point, the stack looks like this, for an input of \$25\$:

[[1 4 9 16 25 36] [1 4 9 16 25 36] 0]

We then collect these values into a single list, which yields the list of \$n^2\$ for each \$n \in [-y, y]\$. We then duplicate this list and operate the table operator over the addition command. The table operator takes a dyad, \$g(p, q)\$, and two arrays, \$A\$ and \$B\$. It then takes the Cartesian Product of \$A\$ and \$B\$ and iterates \$g(a, b)\$ over each pair \$(a, b)\$ where \$a \in A\$ and \$b \in B\$.

In this code, this yields the array \$\big[a^2+b^2 \: | \: a, b \in [-y, y]\big]\$. We then compare each element of this list with the input, yielding a boolean array. Finally, we count the number of \$1\$s in this array by taking its sum, and returning that total.

Most of the code should be understandable when paired with the explanation. A few of the overlooked commands:

• Ω : The reverse operator. Takes a dyad and reverses the order of the arguments
• ‽ : The table operator, as described above.
• ¦ : The fold operator. Takes an array and a dyad and reduces by the dyad. ¦+ is an alias for sum.

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.

ASP, 53 + 4 = 57 bytes

#show N:N=#count{o(A,B):k=A**2+B**2,A=-k..k,B=-k..k}.

Answer Set Programming is a logical language, similar to prolog. I use here the Potassco implementation, clingo.

Input is taken from parameters (-ck= is 4 bytes long). Call example:

clingo -ck=25

Output sample:

12

You can try it in your browser ; unfortunately, this method doesn't handle call flags, so you need to add the line #const k=25 in order to make it work.

Axiom 115 bytes

g(n)==(v:INT:=truncate(sqrt(n)::Float);c:=0;for i in -v..v repeat for j in -v..v repeat if i^2+j^2=n then c:=c+1;c)

ungolfed

gg(n:NNI):NNI== 
   v:NNI:=truncate(sqrt(n)::Float)
   c:NNI:=0
   for i in -v..v repeat
      for j in -v..v repeat
          if i^2+j^2=n then c:=c+1
   c

results

(5) -> [i,g(i)]  for i in 0..25
   Compiling function g with type NonNegativeInteger -> NonNegativeInteger
   (5)
   [[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]]
                                      Type: Tuple List NonNegativeInteger

dobious results

(13) -> g(10018)
   (13)  8
                                                    Type: PositiveInteger
(14) -> g(10019)
   (14)  0
                                                 Type: NonNegativeInteger
(15) -> g(10020)
   (15)  0

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

C++, 81 bytes

#define F(x) for(int x=~n;x++<n;)
void f(int n,int&r){r=0;F(a)F(b)r+=a*a+b*b==n;}

Returns via reference parameter. Simply two ranges over [-n,n], so it works for n=0.

PHP, 80 Bytes

for($m=-$a=1+$argv[1];++$m<$a;)for($n=-$a;$n++<$a;)$c+=$a-1==$m**2+$n**2;echo$c;

PHP, 70 bytes, not competing

for($x=-1;$x++<=$n=$argv[1];)$s+=(-($n%($x-~$x)<1))**$x*4;echo$n?$s:1;

algorithm stolen from one of the Python answers ... I forgot which one; wanted to at least partially understand what´s happening before I posted.

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)

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.

Swift 2.0, 110 108 107 bytes

var n=Int(readLine()!)!,c=0;for(var i = -n;i<=n;i++){for(var j = -n;j<=n;j++){if i*i+j*j==n{c++}}};print(c)

First attempt

Ruby, 66 58 56 54 bytes

->n{r=-n..n;r.map{|x|r.count{|y|n==x*x+y*y}}.reduce:+}

Thanks to sherlock9.

56 bytes

->n{r=-n..n;r.map{|x|r.count{|y|n==x**2+y**2}}.reduce:+}

58 bytes

->n{c=0;r=(-n..n);r.map{|x|c+=r.count{|y|n==x**2+y**2}};c}

66 bytes

->n{c=0;(-n..n).each{|x|c+=n.downto(-n).count{|y|n==x**2+y**2}};c}

Ungolfed:

-> n {
  r = -n..n
  r.map { |x|
    r.count { |y|
      n == x*x + y*y
    }
  }.reduce:+
}

Usage:

->n{r=-n..n;r.map{|x|r.count{|y|n==x*x+y*y}}.reduce:+}[25]
=> 12

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

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.

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.

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.

C++, 175 bytes

#include<iostream> using namespace std; int main() {int i,j,k=0,n; cin>>n; for(i=1;i<n;++i) for(j=1;j<n;++j) if(i*i+j*j==n) ++k; for(i=1;i<=n;++i) if(i*i==n) ++k; cout<<4*k;}

Ungolfed

#include<iostream>
using namespace std;
int main()
{
int i,j,k=0,n;
cin>>n;
for(i=1;i<n;++i)
 for(j=1;j<n;++j)
  if(i*i+j*j==n)
   ++k;
for(i=1;i<=n;++i)
 if(i*i==n)
  ++k;
cout<<4*k;
}

Pyth, 41 35 33 30 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.

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).

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

C, 224 219 182 175 158 bytes

Non-golfed version (with descriptive variables):

#include <math.h>

main (argc, argv) 
    char **argv;
{

    int n, count, max, i, j;
    n = atoi (argv[1]);
    count = 0;
    max = sqrt ((float) n);

    for (i = -max; i <= max; i++)
        for (j = -max; j <= max; j++)
            if (i * i + j * j == n)
                count++;

    printf ("%d", count);
}

Golfed version (with short variables):

#include<math.h>
main(c,v)char**v;{int n,o,x,i,j;n=atoi(v[1]);o=0;x=sqrt((float)n);for(i=-x;i<=x;i++)for(j=-x;j<=x;j++)if(i*i+j*j==n)o++;printf("%d",o);}

Update Sorry about the multiple edits, I'm learning how to abuse K&R-style C ;)

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

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))

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

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.

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

Japt, 42 37 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.

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!

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

Python, 76 74 bytes

I'm sure this can be golfed more, but I need to get back to work.

lambda n:sum([1for a in range(-n,n+1)for b in range(-n,n+1)if a*a+b*b==n])

Try it online

Thanks to @mathmandan for taking off 2 bytes.

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
  

Perl 5, 52 bytes

56 bytes:

$n=pop;for$a(@a=-$n..$n){map$i+=$_*$_+$a*$a==$n,@a}say$i

If output can be in base 1, then 52 bytes:

$n=pop;for$a(@a=-$n..$n){print$_*$_+$a*$a==$n for@a}

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}

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>

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> 

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. The output of 1 means all values match.

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.

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

Octave, 28 bytes

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

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(9ⁿ) 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(n²):

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.
          :+         Add the squares.
            R=       Compare with R.
                   If = pushed 1, keep the pair.
                ,  Count the kept pairs.

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.
    s              Add the squares.
  qQ               Compare the sum with Q.
                 If q returned True, keep T.
l                Count the kept pairs.

Matlab, 72 bytes

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