IntegerPoints.py

"""IntegerPoints.py

Space-efficient streaming algorithms to generate all integer points (x,y)
of the Euclidean plane, in order by distance from the origin,
as described at http://11011110.livejournal.com/332116.html

D. Eppstein, July 2016.
"""

import unittest
from BucketQueue import BucketQueue
from Sequence import Sequence

def IntegerPointsByDistance():
    """All integer points in order by distance, regardless of sign.
    The space needed to generate the first N points is O(N^{1/3}).
    Each point is generated in constant amortized time."""

    # Start out by generating the two initial points of the hull
    yield (0,0)
    yield (0,1)

    # Data structures: hull, a sequence of edges, and queue, a bucket queue
    # Each point is represented as a tuple (x,y) of integer coordinates
    # Each edge is an object with instance attributes
    #    corner -- the point at its counterclockwise end
    #    unit -- the difference between the corner and the
    #            next grid point that lies on the same edge
    #            (which is not necessarily the other endpoint)
    #    beyond -- the next point beyond the edge, in its rectangle
    class edge:
        def __init__(self,corner,unit,beyond):
            self.corner = corner
            self.unit = unit
            self.beyond = beyond

    e = edge((0,0),(0,1),(-1,0))
    f = edge((0,1),(0,-1),(1,0))
    hull = Sequence([e,f])

    def dist2(p):
        """Squared Euclidean distance of a point from the origin"""
        return p[0]**2 + p[1]**2

    def prioritize(e):
        """Include edge e in the priority queue with its correct priority"""
        queue[e] = dist2(e.beyond)

    queue = BucketQueue()
    prioritize(e)
    prioritize(f)

    def box(p,u,b):
        """Given an edge with corner p and unit u, find a boxed translate of b."""
        dot = (b[0]-p[0])*u[0]+(b[1]-p[1])*u[1]
        shift = dot//dist2(u)
        b = (b[0]-shift*u[0],b[1]-shift*u[1])
        if abs(u[0]) + abs(u[1]) == 1:  # Unit square has 2 boxed xlates, pick best
            c = (b[0]+u[0],b[1]+u[1])
            if dist2(c) < dist2(b):
                return c
        return b

    def nonconvex(e):
        """Are e and its successor not strictly convex?"""
        u = e.unit
        v = hull.successor(e).unit
        return u[1]*v[0]-u[0]*v[1] <= 0

    def pop(e):
        """Merge e into its successor (removing the successor from the hull)"""
        f = hull.successor(e)
        if e.unit == f.unit:    # Special case flat vertex
            if dist2(f.beyond) < dist2(e.beyond):
                e.beyond = f.beyond
        else:                   # Concave vertex, use opp vertex of parallelogram
            p = e.corner
            q = f.corner
            r = hull.successor(f).corner
            e.unit = (r[0]-p[0],r[1]-p[1])
            e.beyond = (p[0]-q[0]+r[0],p[1]-q[1]+r[1])
        hull.remove(f)
        del queue[f]
        prioritize(e)
    
    def split(e):
        """Update the hull after adding e's beyond point"""
        p = e.corner
        q = hull.successor(e).corner
        u = e.unit
        b = e.beyond
        e.unit = (b[0]-p[0],b[1]-p[1])
        e.beyond = box(p,e.unit,(b[0]-u[0],b[1]-u[1]))
        fu = (q[0]-b[0],q[1]-b[1])
        fb = box(b,fu,(q[0]+u[0],q[1]+u[1]))
        f = edge(b,fu,fb)
        hull.insertAfter(e,f)
        prioritize(e)
        prioritize(f)
        while nonconvex(hull.predecessor(e)):
            pop(hull.predecessor(e))
            e = hull.predecessor(f)
        while nonconvex(f):
            pop(f)

    for e in queue:
        yield e.beyond
        split(e)

# ============================================================
#     If run from command line, perform unit tests
# ============================================================

class IntegerPointsByDistanceTest(unittest.TestCase):
    radius = 100
    threshold = radius**2
    trange = range(-radius,radius+1)

    def testOrdered(self):
        """Test whether point ordering is by distance."""
        oldDistance = 0
        for x,y in IntegerPointsByDistance():
            newDistance = x**2 + y**2
            self.assertTrue(newDistance >= oldDistance)
            oldDistance = newDistance
            if newDistance > IntegerPointsByDistanceTest.threshold:
                break

    def testCircle(self):
        """Test whether each point in a disk is listed exactly once."""
        points = []
        for x,y in IntegerPointsByDistance():
            if x**2 + y**2 > IntegerPointsByDistanceTest.threshold:
                break
            points.append((x,y))
        self.assertEqual(len(points),len(set(points)))
        self.assertEqual(len(points),len(
            [None for x in IntegerPointsByDistanceTest.trange
             for y in IntegerPointsByDistanceTest.trange
             if x**2 + y**2 <= IntegerPointsByDistanceTest.threshold]))

if __name__ == "__main__":
    unittest.main()