Let’s say you have an array , containing (~millions) points . Perhaps it’s the output of an n-body simulation or something more complicated. Anyway, suppose you also have several other arrays of size , each listing some quantity that is associated with each point – for example , the mass at each point.

The task: find some small contiguous region , and calculate some function at every point within it. For example, approximate the integral of the density by summing over the mass at each point:

Let us introduce an operation , where , contains only zeroes or ones, and contains iff ; also , but maintains the order of elements from . Assume that this operation can be calculated efficiently (in parallel), as opposed to sequentially traversing the elements of an array, which is much slower (this is the case in the IDL programming language, for example).

The solution to the problem at hand is to produce a mask containing elements with values the characteristic function at each point – that is, if and otherwise. Then the final calculation is

and the advantage is that we can reuse the mask for additional calculations.

Now suppose that is itself a very large set, and that we have another, smaller region on which we would also like to perform calculations. Let and . Then we can compose our masks to find the values that lie within .

What are the -indices of points that are in the region ? To find the answer we do:

and now contains the desired -indices.

In the IDL programming language, is written as

A = B[where(C)]

A more general (and also efficient) way of dealing with large, unstructured lists of coordinates is to use an oct-tree.