def compute_max_saturation(a: float, b: float) -> float:
    """
    Finds the maximum saturation possible for a given hue that fits in sRGB.

    Saturation here is defined as `S = C/L`.
    `a` and `b` must be normalized so `a^2 + b^2 == 1`.
    """

    # Max saturation will be when one of r, g or b goes below zero.

    # Select different coefficients depending on which component goes below zero first.

    if (-1.88170328 * a - 0.80936493 * b) > 1:
        # Red component
        k0 = 1.19086277
        k1 = 1.76576728
        k2 = 0.59662641
        k3 = 0.75515197
        k4 = 0.56771245
        wl = 4.0767416621
        wm = -3.3077115913
        ws = 0.2309699292

    elif (1.81444104 * a - 1.19445276 * b) > 1:
        # Green component
        k0 = 0.73956515
        k1 = -0.45954404
        k2 = 0.08285427
        k3 = 0.12541070
        k4 = 0.14503204
        wl = -1.2684380046
        wm = 2.6097574011
        ws = -0.3413193965

    else:
        # Blue component
        k0 = 1.35733652
        k1 = -0.00915799
        k2 = -1.15130210
        k3 = -0.50559606
        k4 = 0.00692167
        wl = -0.0041960863
        wm = -0.7034186147
        ws = 1.7076147010

    # Approximate max saturation using a polynomial:
    sat = k0 + k1 * a + k2 * b + k3 * (a ** 2) + k4 * a * b

    # Do one step Halley's method to get closer.
    # This gives an error less than 10e6, except for some blue hues where the `dS/dh` is close to infinite.
    # This should be sufficient for most applications, otherwise do two/three steps.

    k_l = 0.3963377774 * a + 0.2158037573 * b
    k_m = -0.1055613458 * a - 0.0638541728 * b
    k_s = -0.0894841775 * a - 1.2914855480 * b

    l_ = 1.0 + sat * k_l
    m_ = 1.0 + sat * k_m
    s_ = 1.0 + sat * k_s

    l = l_ ** 3
    m = m_ ** 3
    s = s_ ** 3

    l_ds = 3.0 * k_l * (l_ ** 2)
    m_ds = 3.0 * k_m * (m_ ** 2)
    s_ds = 3.0 * k_s * (s_ ** 2)

    l_ds2 = 6.0 * (k_l ** 2) * l_
    m_ds2 = 6.0 * (k_m ** 2) * m_
    s_ds2 = 6.0 * (k_s ** 2) * s_

    f = wl * l + wm * m + ws * s
    f1 = wl * l_ds + wm * m_ds + ws * s_ds
    f2 = wl * l_ds2 + wm * m_ds2 + ws * s_ds2

    sat = sat - f * f1 / ((f1 ** 2) - 0.5 * f * f2)

    return sat