Gaussian integer division reminder
Python 2 (PyPy), O(1) operations



def cmod(a1, a2, b1, b2):
&#9;# Compute p = a * conj(b)
&#9;p1 =   a1 * b1 + a2 * b2
&#9;p2 = - a1 * b2 + a2 * b1

&#9;# Compute the norm-squared of b
&#9;b_nsq = b1 ** 2 + b2 ** 2
&#9;b_nsq_half = b_nsq // 2

&#9;# Do a symmetric mod of d1 and d2 by b_nsq
&#9;# That is, into the interval [-b_nsq_half, +b_nsq_half)
&#9;d1 = (p1 + b_nsq_half) % b_nsq - b_nsq_half
&#9;d2 = (p2 + b_nsq_half) % b_nsq - b_nsq_half

&#9;# Compute r = d / b
&#9;r1 = (d1 * b1 - d2 * b2) // b_nsq&#32;
&#9;r2 = (d1 * b2 + d2 * b1) // b_nsq&#32;

&#9;return r1, r2


Try it online!

Uses pure integer arithmetic, with no loops or branching or built-in complex numbers. This should easily complete the benchmark of 100,000 random inputs within a small fraction of a second. Porting this to a faster language would be straightforward and likely improve performance.

One possible low-level optimization is the Karatsuba algorithm trick to do complex multiplication using 3 rather than 4 integer multiplications. But, this might not be worth it for such "small" values.