Count sums of two squares

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

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.

Answer 2

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)

Answer 3

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 4

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

Answer 5

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

Answer 6

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.

Answer 7

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 8

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.

Answer 9

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.

Answer 10

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 11

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.

Answer 12

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.

Answer 13

Octave, 28 bytes

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

Answer 14

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 15

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.

Answer 16

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 17

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>

Answer 18

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 19

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.

Answer 20

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 21

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 22

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.

Answer 23

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 24

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.

Answer 25

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 26

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.

---

Answer 27

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.

Answer 28

Matlab, 72 bytes

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

Answer 29

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 30

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