# Generating all integer lattice points (Python)

### From LiteratePrograms

We will describe how to generate all integer lattice points; that is, all tuples (*x _{1}*,

*x*, ...,

_{2}*x*) where the variables run through all integer values {..., -3, -2, -1, 0, 1, 2, 3, ...}. There is obviously more than one possible way to order the sequence of such tuples, but the most sensible choice is to choose an origin such as (0, 0, ..., 0) and then proceed outwards, generating points with monotonically increasing distance from the origin. Even with this restriction, there are multiple approaches depending on the metric used to calculate distance.

_{n}## Contents |

## Generating Z^{n} from N^{n}

It is sufficient to initially generate all points with nonnegative coordinates (), since the rest () can be obtained by reflecting each point across all combinations of the coordinate axes. In the case of *n* = 2, there are three such reflections for each point and the routine looks as follows:

<<generate Z2>>=defgenerate_Z2(limit=None, origin=(0, 0)):forpingenerate_N2(limit): p = (p[0]+origin[0], p[1]+origin[1]) yield pifp[0] != 0: yield -p[0], p[1]ifp[1] != 0: yield p[0], -p[1]ifp[0]andp[1]: yield -p[0], -p[1]

The function `generate_N2` called here is assumed to generate points up to a distance of at most `limit` units from the origin (continuing forever if set to `None`). A custom origin other than (0, 0) may be provided. Note that points lying on the *x* and *y* axes should not be reflected since they would be reflected onto themselves and thus yielded twice; hence the checks for nonzero coordinates.

## In order of Manhattan distance

### In two dimensions

The Manhattan distance is the simplest metric along a perpendicular grid, for two points (*a*, *b*) and (*c*, *d*) defined as |*a*-*c*| + |*b*-*d*|. In other words, the Manhattan distance is the minimum distance required to walk from *a* to *b* if one can only move strictly horizontally or vertically.

For points (*x*, *y*) located in the first quadrant (), the distance *x*+*y* is simply the sum *x*+*y*. (The geometric interpretation of this formula is that points with the same distance lie on a 45° diagonal whose side faces the origin.) Thus, to generate all points in the first quadrant, we may loop through all distances *d* = 0, 1, 2, ... and for each *d* find all pairs *x* and *y* that sum to *d*. An implementation in Python might look like this:

<<generate N2, Manhattan distance>>=defgenerate_N2(limit=None): d = 0whilelimitisNoneord <= limit:forxinrange(d+1): yield x, d-x d += 1

We create a module for this generator:

<<manhattan2d.py>>=generate N2, Manhattan distance generate Z2

Here is an animation of the output from this generator:

### In *n* dimensions

To be done.

## In order of Euclidean distance

### In two dimensions

The Euclidean distance to the origin of a point (*x*, *y*) is defined as . With this metric, the points that share a common distance from the origin are those that lie on the circle with radius *d*. If we generate points by increasing Euclidean distance, they will form the interior of an expanding circle.

We define a helper function that computes the square of the Euclidean distance between a point and the origin. (Taking the square root of this expression is unnecessary, and if we skip it, we only have to work with integer arithmetic.)

<<distance calculation>>=defdistance(p):returnp[0]**2 + p[1]**2

The main generator will yield all points in the first quadrant (). We will treat this quadrant as a list of vertical columns, one for each *x* coordinate. We now create a list whose *n*-th entry holds the maximum *y*-coordinate of any point with that *x*-coordinate yielded so far; that is, each entry will represent how far up the corresponding column has been "filled". For each distance *d* = 0, 1, 3, ..., we proceed by finding all points whose distance from the origin is in the interval [*d*-1, *d*] and store them in a list `yieldable`. Then we sort the list of acquired points by distance and yield them in order.

<<generate N2, Euclidean distance>>=defgenerate_N2(limit=None): ymax = [0] d = 0whilelimitisNoneord <= limit: yieldable = [] find all remaining points closer than d sort list by distanceforpinyieldable: yield p d += 1 ymax.append(0) # Extend to make room for column[d]

To find the intended points, we create an infinite loop that repeatedly fetches points whose distance is in the sought interval. The loop ends when there are no more points to yield.

<<find all remaining points closer than d>>=while1: batch = []forxinrange(d+1): y = ymax[x]ifdistance((x, y)) <= d**2: # Note: distance squared batch.append((x, y)) ymax[x] += 1ifnotbatch:breakyieldable += batch

Sorting by distance can be achieved by providing a custom `key` function:

<<sort list by distance>>=yieldable.sort(key=distance)

To create a complete package, we re-use the routine that generates from implemented earlier.

<<euclidean2d.py>>=distance calculation generate N2, Euclidean distance generate Z2

The functions can be tested by counting the number of points they produce for given distances:

[len(list(generate_N2(n))) for n in range(20)] [len(list(generate_Z2(n))) for n in range(20)]

The respective outputs should be 1, 3, 6, 11, 17, 26, 35, 45, 58, ..., and 1, 5, 13, 29, 49, 81, 113, 149, 197, ...; corresponding to A000603 ("Number of nonnegative solutions to *x*² + *y*² ≤ *n*") and A000328 ("Number of points of norm ≤ *n*² in square lattice.") in the *On-Line Encyclopedia of Integer Sequences*.

The animated output is the following:

### In *n* dimensions

The general case can be handled similarly to the two-dimensional case. In *n* dimensions, the squared Euclidean distance becomes , in Python:

<<squared Euclidean distance, n dimensions>>=defdistance(p):returnsum(x**2forxinp)

Instead of slicing a two-dimensional plane into a one-dimensional list of "heights", we will now split the *n*-dimensional space into an (*n*-1)-dimensional array ("*x* coordinates") of heights ("*y* coordinates"). Python does not have an efficient multidimensional array type built-in (see SciPy), but we can easily simulate a multidimensional array by indexing a `dict` with tuples. Values that have not been set should default to zero, for which purpose we can use Python 2.5's `defaultdict` class:

<<custom dictionary class>>=fromcollectionsimportdefaultdict

We also need a way to generate all array indices up to a bound *m*. This can be implemented by counting to *m*^{n} and writing the number as a base-*m* integer with *n* digits:

<<multirange>>=defmultirange(dims, m): c = 0foriinxrange(m**dims): x = i r = []forninrange(dims): x, rem = divmod(x, m) r.append(rem) yield tuple(r)[::-1]

The main function becomes:

<<generate Nn, Euclidean distance>>=defgenerate_Nn(n=2, limit=None): custom dictionary class ymax = defaultdict(int) # unknown keys are mapped to 0 d = 0whilelimitisNoneord <= limit: yieldable = []while1: batch = []forxinmultirange((n-1), d+1): y = ymax[x] pt = x + (y,)ifdistance(pt) <= d**2: batch.append(pt) ymax[x] += 1ifnotbatch:breakyieldable += batch yieldable.sort(key=distance)forpinyieldable: yield p d += 1

The only essential differences from the two-dimensional case are the iteration over an (*n*-1)-dimensional `multirange`, and the use of tuple concatenation (`x + (y,)`).

To test, we can verify that all generated points are correctly ordered, and check the number of generated points against known data. The sequences for *n* = 3, 4, 5, 6, and 7 are A000604, A055403, A055404, A055405, and A055406.

<<test generate Nn, Euclidean distance>>=testdata ={\ 0 : [1, 1, 1, 1, 1, 1, 1, 1], 1 : [1, 2, 3, 4, 5, 6, 7, 8], 2 : [1, 3, 6, 11, 17, 26, 35, 45], 3 : [1, 4, 11, 29, 54, 99, 163, 239], 4 : [1, 5, 20, 70, 165, 357, 688, 1154], 5 : [1, 6, 36, 157, 482, 1203, 2673, 5139], 6 : [1, 7, 63, 337, 1319, 3819, 9763, 21374], 7 : [1, 8, 106, 702, 3390, 11496, 33792, 83877]}if__name__ == "__main__":fordiminrange(8):forliminrange(8): a = list(generate_Nn(n=dim, limit=lim))foriinrange(1, len(a)):assertdistance(a[i-1]) <= distance(a[i])asserttestdata[dim][lim] == len(a)

<<euclidean.py>>=squared Euclidean distance, n dimensions multirange generate Nn, Euclidean distance test generate Nn, Euclidean distance

Although the results are difficult to visualize in higher dimensions, the three-dimensional case is simple. The points generated by `generate_Nn` with *n*=3 make up one-eighth of a growing ball:

Download code |