This paper presents the results of a comprehensive investigation of complex linear physical-layer network (PNC) in two-way relay channels (TWRC). A critical question at relay R is as follows: "Given channel gain ratio $η= 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 $(α,β)$ that minimizes the symbol error rate (SER) of $w_N=αw_A+βw_B$ when we attempt 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 $α,β,w_A, w_B$, and $w_N$ are 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$-dimensional vectors and $α$ and $β$ were represented by $2\times 2$ matrices, corresponds to a subcase of our Gaussian-integer formulation where $q$ is real prime only. Extension to Gaussian prime $q$, where $q$ can be complex, gives us a larger set of signal constellations to achieve different rates at different SNR. (2) We show how to divide the complex plane of $η$ into different Voronoi regions such that the $η$ within each Voronoi region share the same optimal PNC mapping $(α_{opt},β_{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 $(α_{opt},β_{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.