The Sigma Delta Argument Test
I came up with my own solution with the premise of resolving maximum vector magnitude (including equality) by testing the angle for quadrature between the sum and difference of the two vectors:
For the sum $\Sigma$ and difference $\Delta$ of $z_1$ and $z_2$ given as (which is a 2 point DFT)
$\Sigma = z_1 + z_2$
$\Delta = z_1 - z_2$
The angle $\phi$ between $z_1$ and $z_2$ given as the argument of the complex conjugate multiplication of $\Sigma$ and $\Delta$ as $arg(\Sigma\cdot \Delta^*)$ has the following properties (See derivation at bottom of this answer):
For $z_2 < z_1, |\phi| < \pi/2$
For $z_2 = z_1, |\phi| = \pi/2$
For $z_2 > z_1, |\phi| > \pi/2$
Given the convenience of $\pi/2$ boundaries we never need to compute the argument!
So the efficiency in this approach comes down to computing the sum and difference for the two vectors and then being able to efficiently test whether then phase between them is greater than or less than $\pi/2$.
If multipliers were allowed this would be easily resolved by evaluating the real part of the complex conjugate result, thus here is the complete algorithm using a multiplier, and then to meet the objectives of the question, the approach with no multipliers follows.
If Multiplier Can Be Used
First to introduce a baseline algorithm allowing for multipliers:
1) Step 1: Sum $z_1 = I_1+jQ_1$, $z_2 = I_2+jQ_2$:
$\Sigma = I_{\Sigma} + jQ_{\Sigma} = (I_1+I_2) + j(Q_1+Q_2)$
$\Delta = I_{\Delta} + jQ_{\Delta} = (I_1-I_2) + j(Q_1-Q_2)$
2) Step 2: Compute the Real of the complex conjugate product: $\Sigma\cdot\Delta^*$. This is the dot product and the MSB of the result (the sign bit) is the binary answer directly!
$q = I_{\Sigma}I_{\Delta}+Q_{\Sigma}Q_{\Delta}$
3) Step 3: For a ternary result test q:
$q<0 \rightarrow z_2>z_1$
$q=0 \rightarrow z_2=z_1$
$q>0 \rightarrow z_2<z_1$
So this approach provides a binary > or < result with 2 real multipliers and 5 real sums, resulting in a savings compared to comparing to squared magnitudes which requires 4 real multipliers and 3 read adds. This on its own is not notable as a similar mathematical reduction could be directly achieved as the equations are similar (as already pointed out by @Cedron, @MattL, @Olli in their answers at https://dsp.stackexchange.com/questions/62893/efficient-magnitude-comparison-for-complex-numbers/), but included to show its relation to the Sigma Delta Argument Test: The magnitude test directly in similar form is to compare $I^2+Q^2$:
$$q = (I_1I_1+Q_1Q_1)-(I_2I_2+Q_2Q_2)$$
Which can be rewritten as follows to reduce multipliers (or reordered to directly match the equations above):
$$q = (I_1+Q_1)(I_1-Q_1)-(I_2+Q_2)(I_2-Q_2)$$
The No Multiplier Solution
The no multiplier solution is done by efficiently determining the location of an arbitrary complex point on a plane that is bisected by a line that crosses through the origin. Indeed with similar mapping as shown in the multiplier example above, and actual sum and difference computations may not be needed; this is done to introduce the approach and develop the steps prior to further simplification, which would ultimately map to the same process if expanded. With this approach, the objective is simplified to determining if the point is above or to the left of the line, below or to the right of the line or on the line.
This test can be visualized by rotating D by -$\pi/2$ radians which then changes the test for the boundary between S and D to be $0$ and $\pi$. This rotation si done by swapping I and Q and then change the sign on I: $-j(I+jQ) = Q-jI$
In this case, the test is to see if the point given by D lies above the line y = mx where m is the ratio of the imaginary and real terms of S. (where y is the imaginary axis Q and x is the real axis I). Thus a polar coordinate threshold of constant radius is converted to a rectangular coordinate threshold as a linear line going through the origin, providing for a simpler algorithm. The rotation although trivial can likely be eliminated by modifying the test accordingly. (Although helps in visualizing solutions).
The quadrants are rotationaly invariant to the test so I will refer to Quadrant 1 as whatever quadrant S is in to be as shown in the graphic below (there is no need to actually rotate the axis further, we just rotate our own reference). Q2 and Q4 are trivial, if D is in either of these quadrants a decision can be easily made. When D is in Q3 the test is the negative of when D is in Q1, so the alogirthm is now down to the most efficient way to determine if D is above the y=mx dashed line, below the dashed line, or on the dashed line when both D and S are in Quadrant 1.
Thus tests I have in mind to do that would involve comparisons to the $y= I_S$ and $x=Q_S$ boundaries with efficient binary search for the boundary along the D vector by adding and subtracting $I_D/2^n$ and $Q_D/2^n$ to D. It is to be seen if this or other similar efficient line comparison approaches can beat the other algorithms.

I am exploring some alternate approaches to testing the above angle criteria that may be even more efficient, and will detail the pseudo-code - stay tuned...
Mathematical Derivation
Here is the derivation on how the sum and difference leads to an angle test and provides for the more detailed mathematical relationship (to help with sensitivity testing etc):
consider
$$z_1 = A_1e^{j\phi_1}$$
$$z_2 = A_2e^{j\phi_2}$$
Where $A_1$ and $A_2$ are positive real quantities representing the magnitude of $z_1$ and $z_2$ and $\phi_1$ and $\phi_2$ are the phase in radians.
Divide both by $z_1$ to have expression for $z_2$ relative to $z_1$
$$z_1' = \frac{z_1}{z_1} = 1$$
$$z_2' = \frac{z_2}{z_2} = \frac{A_2}{A_1}e^{j\phi_2-\phi_1} = Ke^{j\phi}$$
Such that if $K>1$ then $z_2>z_1$
The sum and the difference of the $z_1'$ and $z_2'$ would be:
$$\Sigma = z_1' + z_2' = 1 + Ke^{j\phi}$$
$$\Delta = z_1' - z_2' = 1 - Ke^{j\phi}$$
The complex conjugate multiplication of two vectors provides for the angle difference between the two; for example:
Given
$$v_1= V_1e^{j\theta_1}$$
$$v_2= V_2e^{j\theta_2}$$
The complex conjugate product is:
$$v_1v_2^*= V_1e^{j\theta_1}V_2e^{-j\theta_2}= V_1V_2e^{j(\theta_1-\theta_2)}$$
So the complex conjugate product of $\Sigma$ and $\Delta$ with a result $Ae^{j\theta}$ is:
$$
\begin{aligned}
Ae^{j\theta} &= \Sigma \cdot \Delta^* \\
&= (1+Ke^{j\phi})(1-Ke^{j\phi})^* \\
&= (1+Ke^{j\phi})(1-Ke^{-j\phi)}) \\
&= 1 +K(2jsin(\phi))-K^2 \\
&= (1 - K^2) +j2Ksin(\phi) \\
\end{aligned}
$$
Noting that the above reduces to $2jsin(\phi)$ when K= 1, and when K < 1 the real component is always positive and when K > 1 the real component is always negative such that:
for $K < 1, |\theta| < \pi/2$
for $K = 1, |\theta| = \pi/2$
for $K > 1, |\theta| > \pi/2$
Below shows the results of a quick simulation to demonstrate the result summarized above where a uniformly random selection of complex $z_1$, $z_2$ pairs as plotted in the upper plot as red and blue dots, and the resulting mapping to the angle between the sum and difference of $z_1$ and $z_2$.

