mpn/generic/tdiv_qr.c

1 /* mpn_tdiv_qr -- Divide the numerator (np,nn) by the denominator (dp,dn) and
2    write the nn-dn+1 quotient limbs at qp and the dn remainder limbs at rp.  If
9    2. nn >= dn, even if qxn is non-zero.  (??? relax this ???)
11    The time complexity of this is O(qn*qn+M(dn,qn)), where M(m,n) is the time
48 	     mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
53   ASSERT (dn >= 0);
54   ASSERT (dn == 0 || dp[dn - 1] != 0);
55   ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, np, nn));
56   ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, dp, dn));
58   switch (dn)
112 	adjust = np[nn - 1] >= dp[dn - 1];	/* conservative tests for quotient size */
113 	if (nn + adjust >= 2 * dn)
119 	    qp[nn - dn] = 0;			  /* zero high quotient limb */
120 	    if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0) /* normalize divisor */
122 		count_leading_zeros (cnt, dp[dn - 1]);
124 		d2p = TMP_ALLOC_LIMBS (dn);
125 		mpn_lshift (d2p, dp, dn, cnt);
141 	    invert_pi1 (dinv, d2p[dn - 1], d2p[dn - 2]);
142 	    if (BELOW_THRESHOLD (dn, DC_DIV_QR_THRESHOLD))
143 	      mpn_sbpi1_div_qr (qp, n2p, nn, d2p, dn, dinv.inv32);
144 	    else if (BELOW_THRESHOLD (dn, MUPI_DIV_QR_THRESHOLD) ||   /* fast condition */
146 		     (double) (2 * (MU_DIV_QR_THRESHOLD - MUPI_DIV_QR_THRESHOLD)) * dn /* slow... */
147 		     + (double) MUPI_DIV_QR_THRESHOLD * nn > (double) dn * nn)    /* ...condition */
148 	      mpn_dcpi1_div_qr (qp, n2p, nn, d2p, dn, &dinv);
151 		mp_size_t itch = mpn_mu_div_qr_itch (nn, dn, 0);
153 		mpn_mu_div_qr (qp, rp, n2p, nn, d2p, dn, scratch);
158 	      mpn_rshift (rp, n2p, dn, cnt);
160 	      MPN_COPY (rp, n2p, dn);
171 	       Divide a numerator N with nn limbs by a denominator D with dn
172 	       limbs forming a quotient of qn=nn-dn+1 limbs.  When qn is small
173 	       compared to dn, conventional division algorithms perform poorly.
209 	    qn = nn - dn;
215 		MPN_COPY (rp, np, dn);
220 	    in = dn - qn;		/* (at least partially) ignored # of limbs in ops */
224 	    if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0)
226 		count_leading_zeros (cnt, dp[dn - 1]);
354 	    tp = TMP_ALLOC_LIMBS (dn);
361 		    ASSERT_ALWAYS (rn == dn);
370 	    MPN_COPY (rp + in, n2p, dn - in);
379 		mpn_add_n (rp, rp, dp, dn);