p = p_series(a, x, with_error)
            p, p_err = p if with_error
            result  = 1.0 - p
            return with_error ? [result, p_err + 2.0 * Float::EPSILON * result.abs] : result
          elsif a >= 1.0e+06 && (x - a) * (x - a) < a # difficult asymptotic regime, only way to do this region