docs: def mpfr_floordiv(rop, x, y, rnd):
      """
      Given two MPFR numbers x and y, compute floor(x / y),
      rounded if necessary using the given rounding mode.
      The result is placed in 'rop'.
      """
      # Algorithm notes
      # ---------------
      # A simple and obvious approach is to compute floor(x / y) exactly, and
      # then round to the nearest representable value using the given rounding
      # mode.  This requires computing x / y to a precision sufficient to ensure
      # that floor(x / y) is exactly representable.  If abs(x / y) < 2**r, then
      # abs(floor(x / y)) <= 2**r, and so r bits of precision is enough.
      # However, for large quotients this is impractical, and we need some other
      # method.  For x / y sufficiently large, it's possible to show that x / y
      # and floor(x / y) are indistinguishable, in the sense that both quantities
      # round to the same value.  More precisely, we have the following theorem:
      #
      # Theorem.  Suppose that x and y are nonzero finite binary floats
      # representable with p and q bits of precision, respectively.  Let R be any
      # of the IEEE 754 standard rounding modes, and choose a target precision r.
      # Write rnd for the rounding operation from Q to precision-r binary floats
      # with rounding mode R.  Write bin(x) for the binade of a nonzero float x.
      #
      # If R is a round-to-nearest rounding mode, and either
      #
      # (1) p <= q + r and |x / y| >= 2^(q + r), or
      # (2) p > q + r and bin(x) - bin(y) >= p
      #
      # then
      #
      #    rnd(floor(x / y)) == rnd(x / y)
      #
      # Conversely, if R is a directed rounding mode, and either
      #
      # (1) p < q + r and |x / y| >= 2^(q + r - 1), or
      # (2) p >= q + r and bin(x) - bin(y) >= p
      #
      # then again
      #
      #    rnd(floor(x / y)) == rnd(x / y).
      #
      # Proof.  See separate notes and Coq proof in the float-proofs
      # repository.
      #
      # Rather than distinguish between the various cases (R directed
      # or not, p large versus p small) above, we use a weaker but
      # simpler amalgamation of the above result:
      #
      # Corollary 1. With x, y, p, q, R, r and rnd as above, if
      #
      #     |x / y| >= 2^max(q + r, p)
      #
      # then
      #
      #     rnd(floor(x / y)) == rnd(x / y).
      #
      # Proof. Note that |x / y| >= 2^p implies bin(x) - bin(y) >= p,
      # so it's enough that |x / y| >= 2^max(p, q + r) in the case of
      # a round-to-nearest mode, and that |x / y| >= 2^max(p, q + r - 1)
      # in the case of a directed rounding mode.
  
      # In special cases, it's safe to defer to mpfr_div: the result in
      # these cases is always 0, infinity, or nan.
      if not mpfr.mpfr_regular_p(x) or not mpfr.mpfr_regular_p(y):
          return mpfr.mpfr_div(rop, x, y, rnd)
  
      e = _quotient_exponent(x, y)
  
      p = mpfr.mpfr_get_prec(x)
      q = mpfr.mpfr_get_prec(y)
      r = mpfr.mpfr_get_prec(rop)
  
      # If e - 1 >= max(p, q+r) then |x / y| >= 2^(e-1) >= 2^max(p, q+r),
      # so by the above theorem, round(floordiv(x, y)) == round(div(x, y)).
      if e - 1 >= max(p, q + r):
          return mpfr.mpfr_div(rop, x, y, rnd)
  
      # Slow version: compute to sufficient bits to get integer precision.  Given
      # that 2**(e-1) <= x / y < 2**e, need >= e bits of precision.
      z_prec = max(e, 2)
      z = mpfr.Mpfr_t()
      mpfr.mpfr_init2(z, z_prec)
  
      # Compute the floor exactly. The division may set the
      # inexact flag, so we save its state first.
      old_inexact = mpfr.mpfr_inexflag_p()
      mpfr.mpfr_div(z, x, y, mpfr.MPFR_RNDD)
      if not old_inexact:
          mpfr.mpfr_clear_inexflag()
  
      # Floor result should be exactly representable, so any rounding mode will
      # do.
      ternary = mpfr.mpfr_rint_floor(z, z, rnd)
      assert ternary == 0
  
      # ... and round to the given rounding mode.
      return mpfr.mpfr_set(rop, z, rnd)