Repository: sympy/sympy
File: sympy/polys/rootconditions.py (Python)

Code before:
            factors -= factors_prev

        nonzeros -= factors_prev

    return (sign, factors), nonzeros


def dup_schur_conditions(f: dup[Er], K: Domain[Er]) -> list[Er]:
    # Check if -1 is a root, since the transformation is not defined in this case
    evaluation = K.to_sympy(dup_eval(f, -K.one, K))
    if evaluation.is_number and math.isclose(evaluation, 0.0):
        return [-K.one]

    excluded_domains = [EX, EXRAW]
    if not K in excluded_domains and I in K:
        raise NotImplementedError(
            "Schur conditions is not implemented for complex domains"
        )
    if K.is_RR: #XXX: Precision problems with things like RR[x] still remain
        pq = dup_convert(f, K, QQ)
        _, pz = dup_clear_denoms(pq, QQ, convert=True)
        p_transformed: list = dup_transform(pz, [ZZ.one, ZZ.one],
                                            [-ZZ.one, ZZ.one], ZZ)
        conds: list = _dup_routh_hurwitz_fraction_free(p_transformed, ZZ)

        return [K.convert_from(c, ZZ) for c in conds]

    p_transformed = dup_transform(f, [K.one, K.one], [-K.one, K.one], K)
    return dup_routh_hurwitz(p_transformed, K)


Review comment (bug):
There are some precision problems:

```python
In [7]: Poly((s - 1)*(s - 0.3)*(s+0.4), s, domain = EXRAW).schur_conditions()
Out[7]: [6.9888, 8.7808, 1.08801856413265e-16] # wrong result

In [8]: Poly((s - 1)*(s - 0.3)*(s+0.4), s).schur_conditions()
Out[8]: [17472.0, 21952.0, 0.0]

In [9]: Poly((s - 1)*(s - 0.3)*(s+0.5), s, domain = EXRAW).schur_conditions()
Out[9]: [5.98, 9.66, -1]

```

This seems to occour specifically with floats like 0.4, 0.2, ..., which have a theoretically infinite binary representation.

We could consider applying the transformation only after converting every domain into `ZZ`, but I’m not sure if that’s always possible, and it could come with a significant computational cost.

Code after:

        nonzeros -= factors_prev

    return (sign, factors), nonzeros


def dup_schur_conditions(f: dup[Er], K: Domain[Er]) -> list[Er]:
    """
    Computes the schur stability conditions of ``f``.

    Explanation
    ===========

    Returns expressions ``[e1, e2, ...]`` such that all roots of the
    polynomial lie inside the unit circle if and only if ``ei > 0``
    for all ``i``.

    References
    ==========

    .. [1] https://faculty.washington.edu/chx/teaching/me547/2_1_stability.pdf#:~:text=2.6%20Routh,plane%20Real

    """
    if not f:
        raise ValueError("zero polynomial")

    excluded_domains = [EX, EXRAW]
    if not K in excluded_domains and I in K:
        raise NotImplementedError(
            "Schur conditions is not implemented for complex domains"
        )

    if K.is_Exact:
        return _dup_schur_conditions(f, K)
    else:
        K_exact = K.get_exact()
        pe = dup_convert(f, K, K_exact)

        conds = _dup_schur_conditions(pe, K_exact)
        return [K.convert_from(c, K_exact) for c in conds]


def _dup_schur_conditions(f: dup[Er], K: Domain[Er]) -> list[Er]:
    """
    Computes the Schur stability conditions, mapping the unit circle to the left
    half plane, transforming the problem into a Routh-Hurwitz stability problem.

    The transformation is done via ``z = (1 + s)/(1 - s)`` to the polynomial
    ``f``.

    """