This paper presents the results of a comprehensive investigation of complex linear physical-layer network coding (PNC) in two-way relay channels. In this system, two nodes A and B communicate with each other via a relay R. Nodes A and B send complex symbols, $w_{A}$ and $w_{B}$ , simultaneously to relay R. Based on the simultaneously received signals, relay R computes a linear combination of the symbols, $w_{N}=\alpha w_{A}+\beta w_{B}$ , as a network-coded symbol and then broadcasts $w_{N}$ to nodes A and B. Node A then obtains $w_{B}$ from $w_{N}$ and its self-information $w_{A}$ by $w_{B}=\beta ^{-1}(w_{N}-\alpha w_{A})$ . Node B obtains $w_{B}$ in a similar way. A critical question at relay R is as follows: “given channel gain ratio $\eta = h_{A}/h_{B}$ , where $h_{A}$ and $h_{B}$ are the complex channel gains from nodes A and B to relay R, respectively, what is the optimal coefficients $(\alpha ,\beta )$ that minimizes the symbol error rate (SER) of $w_{N}=\alpha w_{A}+\beta w_{B}$ when the relay attempts to detect $w_{N}$ in the presence of noise?” Our contributions with respect to this question are as follows: 1) we put forth a general Gaussian-integer formulation for complex linear PNC in which $\alpha ,\beta ,w_{A}, w_{B}$ , and $w_{N}$ are the elements of a finite field of Gaussian integers, that is, the field of $\mathbb {Z}[i]/q$ , where $q$ is a Gaussian prime. Previous vector formulation, in which $w_{A}$ , $w_{B}$ , and $w_{N}$ were represented by 2-D vectors and $\alpha $ and $\beta $ were represented by $2\times 2$ matrices, corresponds to a subcase of our Gaussian-integer formulation, where $q$ is real prime only. Extension to the Gaussian prime $q$ , where $q$ can be complex, gives us a larger set of signal constellations to achieve different rates at different values of SNR; and 2) we show how to divide the complex plane of $\eta $ into different Voronoi regions, such that the $\eta $ within each Voronoi region shares the same optimal PNC mapping $(\alpha _{\mathrm{ opt}},\beta _{\mathrm{ opt}})$ . We uncover the structure of the Voronoi regions that allows us to compute a minimum-distance metric that characterizes the SER of $w_{N}$ under optimal PNC mapping $(\alpha _{\mathrm{ opt}},\beta _{\mathrm{ opt}})$ . Overall, the contributions in 1) and 2) yield a toolset for a comprehensive understanding of complex linear PNC in $\mathbb {Z}[i]/q$ . We believe investigation of linear PNC beyond $\mathbb {Z}[i]/q$ can follow the same approach.