It is, in fact, a difficult problem. But not in dimension 2.
From a given basis (a set of independant vectors) of an integer lattice, you are trying to find a shortest non-zero) vector of that lattice. This is the shortest vector problem (SVP), and for larger dimensions it is in fact hard (known results include NP-hardness for randomized reductions and infinity norm).
The problem is that you could potentially get a short vector from a weird integer combination of the basis vectors, as the way they interact with each other is really hard to predict (especially in larger dimensions).
In dimensions larger than two, there are in fact no better known ways to generate a short vector than to enumerate all possible vectors in a given sphere. There are ways to avoid generating too many vectors, but the complexity is still exponential. The LLL algorithm is a polynomial approach that exploits the basic idea of trying to create shorter vectors by working on only two vectors at a time (a generalisation of the 2-dimensional case), but it only gives you a relativly short vector (best proof gives you a factor of $2^{O(n)}$ on the size :().
In dimension two however, the LLL algorithm (Called Lagrange's algorithm, sometimes erroneously attributed to Gauss, a generalisation of the euclidean algorithm.) gives you a basis of the lattice made of the shortest and the second shortest (independant) vector of the lattice in polynomial time.
Short story in dimension 2 : from the largest vector of your basis, you add/remove the shortest one as to make it shorter. Rinse and repeat. (works the same way as the euclidean algorithm to find the pgcd of two integers)
Exemple : from $\displaystyle \binom{5}{5}$ and $\displaystyle \binom{14}{13}$.
you remove two times $\displaystyle \binom{5}{5}$ from $\displaystyle \binom{14}{13}$, forming $\displaystyle \binom{4}{3}$. as an integer combination, this new vector is inside your lattice.
your new vectors are $\displaystyle \binom{5}{5}$ and $\displaystyle \binom{4}{3}$. You remove the latter from the former, and you get $\displaystyle \binom{1}{2}$.
now you have $\displaystyle \binom{1}{2}$ and $\displaystyle \binom{4}{3}$. You remove the former twice from the latter, and get $\displaystyle \binom{2}{-1}$.
now you have $\displaystyle \binom{2}{-1}$ and $\displaystyle \binom{1}{2}$. You can't make any of them shorter by adding/removing the other. The LLL algorithm in dimension 2 (called the Lagrange or Gauss algorithm) tells you that those are the shortest vectors in the lattice (non dependant).
This is all usually done by using unitary transformation (multiplication by integer lattices of determinant $\pm 1$) on the lattice basis (a $n\times n$ matrix). Here we did :
$$ \left(\begin{array}{cc} 5&14 \\5 &13\end{array}\right) \left(\begin{array}{cc}1 & -2\\ 0&1\end{array}\right)\left(\begin{array}{cc}1 &0 \\-1 &1\end{array}\right)\left(\begin{array}{cc}1 & -2\\ 0&1\end{array}\right)=\left(\begin{array}{cc} 1&2 \\ 2&-1\end{array}\right) $$
the original basis is on the left, and each step is one of the multiplying matrices (from left to right). The end result is on the right. There usually are $\left(\begin{array}{cc}0 &1 \\ 1&0\end{array}\right)$ matrices in between for swapping the two vectors (that way the shortest one is on the left), but I removed those here for clarity.
For larger dimensions, it gets ugly. Couldn't find anything nice for the 2 dimensional case, so let's go :
Inner working :
You have two vectors $u$ and $v$, sorted by euclidean length. We define the coefficient $\mu=\dfrac{u\cdot v}{\|u\|^2}$. that coefficient $\mu$ represents the relative length of the projection of $v$ over $u$ (see picture below). Then making the vector $v$ as short as possible is done by removing $\lceil \mu \rfloor \times u$ from $v$. We get a new vector $v$ such that one of the following happens :
the new $\mu$ is in $\left[-\dfrac{1}{2};\dfrac{1}{2}\right]$, meaning that the $v$ vector is inside the two vertical bars. This is called size-reduced.
On the left side the vector $v$ is inside the circle of radius $\|u\|$, meaning that $v$ is now shorter than $u$. Therefore swapping position of $u$ and $v$ will yield $|\mu|>\dfrac{1}{2}$ (you can see it by projecting $u$ on $v$)
On the right side however, $v$ is outside the circle (the Lovasz condition is said to be verified) and therefore swapping the places of $u$ and $v$ wouldn't yield anything. The basis is said LLL-reduced, and $u$ is the shortest non-zero vector of the lattice.
Proof :
let's say it isn't. There are $x$ and $y$ integers ($x\neq 0$) such that $xu+yv$ is shorter than $u$.
$$
\begin{aligned}
\|xu+yv\|^2 &= x^2\|u\|^2+2xy(u\cdot v)+y^2\|v\|^2\\
&= x^2\|u\|^2+2\mu xy \|u\|^2+y^2\|v\|^2\\
&\geq (x^2+2\mu xy +y^2)\|u\|^2\\
&\geq (x^2- |xy| +y^2)\|u\|^2
\end{aligned}
$$
Either $x^2$ or $y^2$ is larger than $|xy|$. Therefore $(x^2- |xy| +y^2)$ is larger than $\min(x^2;y^2)$. If that minimum is $0$, $xu+yv$ is a non-zero integer combination of $u$ or $v$, and therefore larger than $u$ (contradiction). If it isn't, then $(x^2- |xy| +y^2)>0$ and we get a contradiction.
Therefore $u$ is the shortest vector.
For LLL in larger dimension, you do the same two vectors at a time, selecting the two vectors using a process that's similar to a bubble sort : once a pair of vectors is done you get to the next vector, and you move down each time you swap vectors. For the $\mu$ coefficient you use the Gram-Schmidt Orthogonalisation.
Best Answer
For cubic cells, I think the way to do it is to calculate distance for just one point component-wise, and if the distance in any component is greater than 1/2, switch to a different point.
I thought about this in 2-d: There's two points, $A$ and $B$. Since both $A$ and $B$ form a grid, because they're periodic, you can forget about the original unit grid and simply imagine that $A$ is randomly placed within a square created by the points of $B$. Now if I draw lines that represent the points equidistant from any two corner points of this $B$ square, they will form a cross in the middle of the square. So the maximum distance that $A$ can be from $B$, component-wise, is 1/2. If I take any $B$ and calculate the distance to $A$, that will tell me what to do:
In 3-d it's the same idea, but you can move in the z-direction.
For non-cubic cells, I can't think of a way that's as quick and simple. One thing you could do is pre-compute the distances between different points in the lattice, so that if you calculate the distance to $A$ and $B_1$ and find that it's less than half of the distance between $B_1$ and $B_i$, then you know that $B_i$ will not be the closest point to $A$, and you didn't have to calculate $B_i$'s position explicitly. There might be a way to do a similar trick as the 3-d one with planes: if you pre-compute vectors $d_{ij}=B_i-B_j$ and then if $d_{ij}\cdot A - \frac{d_{ij}^2}{2} >0$, then $a$ is closer to $B_i$ than $B_j$. This should work because $d_{ij}\cdot r = c$ defines a plane parallel to the plane of points equidistant between $B_i$ and $B_j$; since it should be half-way between $B_i$ and $B_j$ substitute $d_{ij}/2$: $d_{ij}\cdot d_{ij}/2=c$. The problem here is that unlike the simple 3-d case, knowing that $A$ is closer to $B_i$ than $B_j$ doesn't help for any other point $B_k$. Still, 27+ dot products with pre-computed vectors might still be a big speed-up.