# Fast Low-Precision Computer-Generated Holography on GPU
## 2 Methodology
### Packed Separable PSF Phases  
The phase of the PSF \(P(x,y)\) is given by the parabola as follows:  
\[\phi(x,y)=\angle P(x,y)=\frac{\pi}{\lambda\zeta}\big{[}(x-\delta)^{2}+(y- \epsilon)^{2}\big{]}=\frac{\pi}{\lambda\zeta}(x-\delta)^{2}+\frac{\pi}{ \lambda\zeta}(y-\epsilon)^{2} \tag{3}\]  
i.e., scaled by a constant inversely proportional to the axial distance of the point emitter. Note that the phase \(\phi\) is, therefore, separable in \(\phi_{\rm x}\) and \(\phi_{\rm y}\):  
\[\phi_{\rm x}(x)=\frac{\pi}{\lambda\zeta}(x-\delta)^{2};\quad\phi_{\rm y}(y)= \frac{\pi}{\lambda\zeta}(y-\epsilon)^{2};\quad\phi(x,y)=\phi_{\rm x}(x)+\phi_ {\rm y}(y). \tag{4}\]  
This is leveraged to minimize redundant operations when calculating and packing the bits encoding the \(\phi_{\rm x}\) and \(\phi_{\rm y}\) values. When using fixed-point integers, one should carefully select the units so as to optimize the number of required bits for encoding numbers in a given application, balancing precision and computational performance constraints. For the conversion to \(n\)-bit (fixed-point) integer representations, we want to encode units as integer multiples of \(\pi\cdot 2^{1-n},n\in\mathbb{N}_{0}\); we make this choice since the exponential of purely imaginary functions is periodic (consisting of sines/cosines), so we can ignore any overflow when adding phase values together, which is some multiple of \(2\pi\). The precision of the quantization is denoted as \(s\), which should be no larger than the hologram pixel pitch \(p\) in order to ensure numerically distinguishable signal changes across subsequent pixels; otherwise, the coordinates for two neighboring pixels may be quantized (e.g. rounded) to the same value, thereby resulting in erroneous identical computed PSF amplitudes. One may want to choose \(s<p\) if finer lateral quantizations of the point cloud would be desirable, e.g., for sub-pixel lateral resolution of point coordinates for very detailed objects. These precision considerations must also be taken to determine the axial quantization of the factor \(\frac{\pi}{\lambda\zeta}\), which we denote as \(c\), i.e., the depth resolution of the point cloud. It should also be proportional to some \(\pi\cdot 2^{1-m},m\in\mathbb{N}_{0}\) so to facilitate integer multiplication, but not necessarily the same as the precision chosen for \(s\). We chose the same precision for this parabolic phase scaling factor as that for \(s\) because of GPU register size constraints, giving us the following integer:  
\[c=\left\lfloor\frac{\pi}{\lambda\zeta}\frac{s^{2}2^{n}}{2\pi}\right\rceil= \left\lfloor\frac{s^{2}}{2\lambda\zeta}\right\rfloor \tag{5}\]  
where \(\lfloor\cdot\rceil\) is the rounding operator. This process is summarized for \(\phi_{\mathrm{x}}\) in Figure 1; it is identical for \(\phi_{\mathrm{y}}\), except for the other chosen coordinate difference \((y-\epsilon)\). The final multiplication result can have a different bit-width \(k<n\) than the intermediate integers to significantly reduce the multiplier design in, for example, a customized FPGA implementation. Because typical GPUs are more constrained, namely, memory is addressed byte-wise and registers are 32-bit for most modern GPU models, we chose \(k=8\) and \(n=32\). These low-precision \(\phi_{\mathrm{x}}\) and \(\phi_{\mathrm{y}}\) phase values are tightly packed together in registers, and all pairwise combinations are added together efficiently to obtain the phase \(\phi=\phi_{\mathrm{x}}+\phi_{\mathrm{y}}\) in every pixel, as is detailed further in the \"Implementation\" section. Figure 1: Diagram of the fixed-point integer multiplications for computing the phase. (a) Fixed-point integer multiplication diagram with example bit patterns; the boxes indicate where significant bits are found (all other bits are assumed to be 0). (b) Proposed scheme for computing phase values. The bits encode the binary digits representing the negative powers of two, multiplied by \(2\pi\).