## Python (59 bytes)

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

[Online demo](http://ideone.com/dXMOrd)

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