For discrete-time systems, stability_margins with method='poly', and so with the default 'best', returns gain crossovers where |L| is nowhere near 1, and with them a wrong phase margin. It starts at modest orders. For a moving average, L = T/(1 - T) with T(z) = (1/N)(z^-1 + ... + z^-N):

(control 0.10.2, numpy 2.5.4. "True" is from a dense grid with the one gain crossing refined by bisection.)

import warnings
import numpy as np
import control as ct

warnings.simplefilter("ignore")
for N in (8, 17, 144):
    num = np.r_[0.0, np.full(N, 1.0 / N)]          # T over z^N
    L = ct.tf(num, np.polysub(np.r_[1.0, np.zeros(N)], num), dt=1)
    for method in ("poly", "best", "frd"):
        gm, pm, sm, wpc, wgc, wms = ct.stability_margins(L, method=method)
        print(N, method, round(pm, 3), "|L| at wgc:", round(abs(complex(L(np.exp(1j * wgc)))), 4))

At N = 144, poly reports wgc ≈ 3.098 rad/sample, where |L| ≈ 0.0001.

Cause. _poly_z_mag1_crossing and _poly_z_real_crossing keep a root of num(z)num(1/z) - den(z)den(1/z) if ||z| - 1| < eps**(1/len(p2)) (margins.py lines 119 and 135 on main). That tolerance grows with degree: 0.36 at N = 17, 0.88 at N = 144. At N = 144, 143 roots with |z| from 0.962 to 1.040 are kept, and nothing checks |L| at them. best doesn't fall back, because _likely_numerical_inaccuracy compares coefficient norms, not degree.

Suggested fix. Accept a candidate only if abs(abs(L(exp(1j*w))) - 1) (or the imaginary part, for phase crossings) is small, to a tolerance independent of degree. Making _likely_numerical_inaccuracy consider degree would also send best to frd.

Workaround. Pass an FRD on a grid that covers the crossover. Plain method='frd' can miss a low-frequency one: for 0.02 / (z^5 (z - 1)) it returns pm = inf, where poly is right.

Related but different: #523, fixed for a low-order example.