Number of affine points on an elliptic curve I'm given an elliptic curve $y^2=x^3+ax+b \in \mathbb{Z}_p[x]$ (with numbers $a,b,p$ not greater than $10^6$).

I would like to find, using the naive approach, the number of affine points on the curve without the point at infinity . 

My approach is to go through all $x\in\mathbb{Z}_p$ and check which of them have residue (stored in counter. Since each $x$ has either two or zero residues, I multiply the number of $x$ with residues times 2.

To return the number of affine points, the following method is used:

if isSingular(a,b,p):
    return counter*2+2
return counter*2+1


It seems to work for some cases, but not for others. Is there something special we have to consider for singular curves, and why?
The final code that worked for is shown below (published here)

# Returns number of affine points on the elliptic curve (wihout point at inf).
def countResidues(a, b, p):

    counter  = 0
    for x in range(0,p):

        # the curve in 
        y2 = (x**3+a*x+b) % p

        # if x^3+ax+b = 0 (mod p), then there exists one point only.
        # i.e. the point is still on the curve even if there is no residue.
        if y2 == 0:
            counter += 1
            continue

        # if quadratic residues exists, two points are on the curve
        # i.e. the solutions to y^2 = x^3+ax+b (mod p)
        if isQuadRes(y2, p):
            counter +=2

    return counter


# Check if a curve is singular by checking that 4a^3+27b^2 = 0
# See p. 255 in course lit.
def isSingular(a, b, p):
    return ( -(4*a*a*a % p)-(27*b*b % p)) % p == 0


# Determines if there exists an x such that x^2=a (mod p), in other words,
# it will return true if the quadratic residue exists.
# Based on Euler's criterion, see course lit p. 180.
def isQuadRes(a,p):
    return squareAndMultiply(a,((p-1)//2), p) == 1

# Convert integer to binary representation
def intToBin(number):
    return '{:b}'.format(number)

# Computes number mod base
def mod(number, base):
    return number - base * (number//base)

# Calculate 'x^c mod n' using square-and-multiply algorithm.
def squareAndMultiply(x,c,n):

    # get the binary representation of c, the exponent.
    b = intToBin(c) 

    z = 1

    # the main algorithm, adjusted from p. 177 from course lit.
    for i in range(len(b)):
        z = mod(z*z, n)

        if b[i] == '1':
            z = mod(z*x, n)
    return z



Course lit refers to https://www.amazon.com/Cryptography-Practice-Discrete-Mathematics-Applications/dp/1584885084