is_square?

from  itertools import combinations

def distsq(p1, p2):
    x1, y1 = p1
    x2, y2 = p2
    return (x2 - x1) ** 2 + (y2 - y1) ** 2

def is_square(p1, p2, p3, p4):
    # check for double points first
    if any(a == b for a, b in combinations((p1, p2, p3, p4), 2)):
        return False
    dis = [distsq(p1, p) for p in (p2, p3, p4)]
    pmin = min(dis)
    pmax = max(dis)
    if not (pmax == 2 * pmin and dis.count(pmin) == 2):
        return False
    # find point furthest away and assert that distances of that point to the 2
    #  remaining points are the same as pmin
    far_pt = (p2, p3, p4)[dis.index(pmax)]
    return all(distsq(far_pt, p) == pmin for p in (p2, p3, p4) if p != far_pt)

if __name__ == "__main__":
    """random test against brute force function (is_square0)"""
    import random
    ran = random.randint
    X = 5

    def sub(a, b):
        return (a[0] -  b[0], a[1] - b[1])

    def add(a, b):
        return (a[0] +  b[0], a[1] + b[1])

    def is_square0(p1, p2, p3, p4):
        """for testing - check everything to death"""
        if any(a == b for a, b in combinations((p1,p2,p3,p4), 2)):
            return False
        dis = [distsq(p1, p) for p in (p2,p3,p4)]
        pmin = min(dis)
        pmax = max(dis)
        if not (pmax == (2 * pmin) and dis.count(pmin) == 2):
            return False
        dis = [distsq(p2, p) for p in (p1, p3,p4)]
        if not (pmin == min(dis) and pmax == max(dis)):
            return False
        dis = [distsq(p3, p) for p in (p1, p2,p4)]
        if not (pmin == min(dis) and pmax == max(dis)):
            return False
        dis = [distsq(p4, p) for p in (p1, p3,p2)]
        if not (pmin == min(dis) and pmax == max(dis)):
            return False
        return True

    def rect():
        b = a = (ran(-X, X), ran(-X, X))
        while a == b:
            b = (ran(-X, X), ran(-X, X))
        dx,dy = sub(b, a)
        c = add(a, (dy, -dx))
        d = add(b, (dy, -dx))
        return a, b, c, d

    def norect():
        """random - still could be rect"""
        return [(ran(-X, X), ran(-X, X)) for i in range(4)]

    for i in range(20):
        q= rect()
        assert is_square(*q), "{:d} {:s}".format(i, str(q))

    for i in range(20000):
        q= norect()
        s1 = is_square(*q)
        s2 = is_square0(*q)
        assert s1 == s2, str(s1) + " " + str(s2) + " " + str(q)
    print "OK"