+
+ if (yl > MPW_MAX) {
+ yh += yl >> MPW_BITS;
+ yl &= MPW_MAX;
+ }
+
+ while (yh > rh || (yh == rh && yl > rrr)) {
+ q--;
+ yh -= d;
+ if (yl < dd) {
+ yh++;
+ yl += MPW_MAX;
+ }
+ yl -= dd;
+ }
+ }
+
+ /* --- Remove a chunk from the dividend --- */
+
+ {
+ mpw *svv;
+ const mpw *dvv;
+ mpw mc = 0, sc = 0;
+
+ /* --- Calculate the size of the chunk --- *
+ *
+ * This does the whole job of calculating @r >> scale - qd@.
+ */
+
+ for (svv = rv + scale, dvv = dv;
+ dvv < dvl && svv < rvl;
+ svv++, dvv++) {
+ mpd x = (mpd)*dvv * (mpd)q + mc;
+ mc = x >> MPW_BITS;
+ x = (mpd)*svv - MPW(x) - sc;
+ *svv = MPW(x);
+ if (x >> MPW_BITS)
+ sc = 1;
+ else
+ sc = 0;
+ }
+
+ if (svv < rvl) {
+ mpd x = (mpd)*svv - mc - sc;
+ *svv++ = MPW(x);
+ if (x >> MPW_BITS)
+ sc = MPW_MAX;
+ else
+ sc = 0;
+ while (svv < rvl)
+ *svv++ = sc;
+ }
+
+ /* --- Fix if the quotient was too large --- *
+ *
+ * This doesn't seem to happen very often.
+ */
+
+ if (rvl[-1] > MPW_MAX / 2) {
+ mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
+ q--;
+ }
+ }
+
+ /* --- Done for another iteration --- */
+
+ if (qvl - qv > scale)
+ qv[scale] = q;
+ r = rr;
+ rr = rrr;
+ }
+ }
+
+ /* --- Now fiddle with unnormalizing and things --- */
+
+ mpx_lsr(rv, rvl, rv, rvl, norm);