Q:

Are there 2D analogues for integer division and modular arithmetic?

Let's say you have a "parallelogram" of points $P = \{(0, 0), (0, 1), (1, 1), (0, 2), (1, 2)\}$. This parallelogram lies between $u = (2, 1)$ and $v = (-1, 2)$.
Then for any point $n \in \mathbb{Z}^2$, you can write
$n = dQ + r$,
where $Q = \begin{pmatrix} u\\ v\end{pmatrix} = \begin{pmatrix}2 &1 \\ -1 &2\end{pmatrix}$, and $r \in P$.
As an example, I try to show how these things relate in the following image:

$$(3, 5) = (2, 1)\begin{pmatrix}2 &1 \\ -1 &2\end{pmatrix} + (0, 1)$$
Here is another image, showing that $Q$ need not be in the form $\begin{pmatrix}a & -b \\b&a  \end{pmatrix}$.

$$Q = \begin{pmatrix}-4 & 2 \\ 1 & 2\end{pmatrix}$$
My question is: is type of thing studied (for general parallelograms), and what is it called?
I have many questions about it, and do not know where I can find more information. For example: how do you find $d$ and $r$, does it make sense (and what is the use) of defining $n_1 \equiv n_2 \pmod Q$ if $n_1$ and $n_2$ have the same "remainder" when "divided" by $Q$, and if so, what rules govern this system; what is the situation on a hexagonal lattice; what if $P$ is not a parallelogram.

Background
The reason I'm interested in this is that I use it in a large variety of algorithms, including sub-sampling of grids, maze generation, and using grids to represent more complicated things (for example, using a hex grid to represent a triangular grid + vertices).
Although I have ways to calculate $d$ and $r$ given $n$ and $Q$, my algorithms are rather clumsy, and "hacked together". I'd like to back this with some theory, perhaps simplify and optimize the algorithms, and be in a better position to design related algorithms and be able to justify it mathematically. 

A:

Rather than using matrices in the quotient there, it may be a better idea to use complex numbers. There, the division algorithm uses complex multiplication instead of matrix multiplication as you have used. For the square lattice, the relevant object is the domain of Gaussian Integers. For the hexagonal (triangular) lattice, you need the Eisenstein integers. Both are Euclidean domains, which means that division with remainder as you have described is possible.
In order to actually compute the division in the Gaussian integers (square lattice), the following is a start: 
Let $z,w \in \mathbb{Z}[i], w\not=0$. We want to find $q,r\in \mathbb{Q}$ such that $z = qw + r$ and $|r| < |w|$, or equivalently, $q$ such that $|\frac{z}{w} - q|< 1$.
If $\frac{z}{w} = x + iy$ where $x,y \in \mathbb{Q}$, take $q = [x] + i[y]$ (where $[x]$ denotes the integer closest to $x$), then indeed we have that $|\frac{z}{w} - q|\leq |(x - [x]) + i(y - [y])|\leq \sqrt{\frac{1}{2^2} + \frac{1}{2^2}} = \frac{1}{\sqrt{2}} < 1$. Hence $q = [x] + i[y]$ and $r = z - qw$ work.
For the Eisenstein integers, you should be able to devise something very similar.
If for some reason you need to have this in the language of matrices, just remember that multiplication of a complex number $x+iy$ by the complex number $a + bi$ is equivalent to multiplication of the vector $\begin{pmatrix}x \\ y\end{pmatrix}$ by the matrix $\begin{pmatrix} a & -b \\ b & a\end{pmatrix}$, so the above method can be converted to vectors in $\mathbb{Z}^2$ and matrices if you prefer.

A:

Finding $d$ and $r$: I don't know whether this is more or less clumsy than your existing algorithm: to write $n$ as $dQ+r$ for a given $Q,$ you can compute $x=nQ^{−1}.$ Then compute the floor of the elements of $x$ to get $d,$ and compute $n−dQ$ to get $r.$ To keep the arithmetic integer, if desired, you can write $Q^{−1}=\mathrm{adj}(Q)/\det Q.$ Since the adjugate matrix, $\mathrm{adj}(Q),$ is integer, the vector $y=n\,\mathrm{adj}(Q)$ is integer, and $d$ can be obtained by dividing $y$ by $\det Q$ using integer division. This algorithm applies in dimensions higher than $2$ as well.
To verify that this works, the first key idea is that $x\mapsto xQ$ is a linear map.  So $(ax+by)Q=a(vQ)+b(uQ)$ for any vectors $x$ and $y$ and scalars $a$ and $b.$  This means that the unit square spanned by $e_1=(1,0)$ and $e_2=(0,1),$ that is, the set $\{ae_1+be_2\mid 0\le a,b<1\},$ is sent to your parallelogram, that is, to $\{au+bv\mid 0\le a,b<1\},$ under this map.
The second key idea is that $Q$ is invertible.  Consider a point $n=dQ+r,$ where $d$ is an integer vector and $r$ is in your parallelogram.  Multiplying  by $Q^{-1}$ from the right on both sides, we get $nQ^{-1}=d+rQ^{-1}.$  To show that the algorithm works, that is, to show that the floor of this expression equals $d,$ we need to show that all components of $rQ^{-1}$ lie in the interval $[0,1)$.  But we have established that $Q$ sends the unit square to your parallelogram.  Since $Q^{-1}$ is invertible, it reverses the process, sending points in the parallelogram to the unit square, which is what we needed to show.
The intuition is that $dQ,$ with $d$ an integer vector, is a lattice point.  If your parallelogram is translated to this lattice point, the other three vertices of the translated parallelogram are $(d+e_1)Q,$ $(d+e_2)Q,$ and $(d+e_1+e_2)Q.$  Points in the interior of the translated parallelogram are $xQ,$ where, writing $d=(d_1,d_2),$ $x=(x_1,x_2),$ we need $d_1<x_1<d_1+1,$ $d_2<x_2<d_2+1.$  So the floor of $x_1$ and $x_2$ are $d_1$ and $d_2.$
Where do $P$ and $Q$ come from? I'm not sure from your question whether you are starting from the set $P$ and obtaining $Q$ from it, or starting from $Q$ and computing $P$ from it, or whether both $P$ and $Q$ are given data.  It seems to me that if you need to compute $P$ from $Q,$ that's a slightly tricky problem—harder, in my opinion than computing $d$ and $r$ from $Q.$  (One can always use the algorithm for finding $r,$ combined with brute force, but for a very long and thin parallelogram this might not be efficient.)
Note that the same lattice can be described by infinitely many different matrices $Q.$  Specifically, a matrix $Q'=QU,$ where $U$ is an integer matrix with determinant $\pm1,$ produces the same lattice as $Q$ does.  Matrices like $U$ are called unimodular matrices.  The product of two unimodular matrices is a unimodular matrix, so unimodular matrices form a group.  This is essentially the modular group (up to questions about whether determinant $-1$ is allowed, and whether $U$ and $-U$ should be considered equivalent).
In your position, however, I suspect a study of the modular group might not be the most direct approach to answering the questions you're interested in.  I would focus more on the computational, algorithmic literature.  I don't have a good recommendation, unfortunately, but you might start by reading up on the LLL algorithm.  Even if you don't have an immediate use for that algorithm, it is likely that someday you will.  But more to the point, you'll probably pick up a lot of useful insights about how people who compute with lattices think about things.
To finish the discussion above: if you use $Q'$ instead of $Q,$ you have a different parallelogram, and therefore a different set of points—call it $P'.$  You can create a dictionary between $P$ and $P',$ using the algorithm for finding $r.$  (Note that $P$ and $P'$ will have the same size; correspondingly, $Q$ and $Q'$ will have the same determinant, up to sign.)
Fundamental region not a parallelogram: If you want to do the same thing for a hexagonal tiling, say, you can proceed in a similar way.  (There's a bit of a terminology issue here: a hexagonal tiling is not a lattice: mathematically, a lattice is a group, which requires that the sum of two lattice vectors be a lattice vector.  This is not true of the vertices of a hexagonal tiling.)  Note that in a parallelogram tiling, the four vertices are equivalent.  So $P$ contains one vertex and a bunch of interior points.  In the hexagonal case, the six vertices are not all equivalent; there are two equivalence classes, consisting of three vertices each.  So $P$ will contain two vertices and a bunch of interior points.
You can always define a parallelogram-shaped fundamental region that gives the same thing as the hexagonal one.  Then use the algorithm above for finding $r.$  Finally you will need a dictionary mapping the points of the fundamental parallelogram to the points of the fundamental hexagon.  Devising an algorithm to produce such a dictionary may be a bit awkward, but once you have the dictionary, things are straightforward.

