1 /* mpn_dcpi1_bdiv_qr -- divide-and-conquer Hensel division with precomputed
2 inverse, returning quotient and remainder.
4 Contributed to the GNU project by Niels Möller and Torbjorn Granlund.
6 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
7 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
8 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
10 Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
12 This file is part of the GNU MP Library.
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
19 The GNU MP Library is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
21 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22 License for more details.
24 You should have received a copy of the GNU Lesser General Public License
25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
31 /* Computes Hensel binary division of {np, 2*n} by {dp, n}.
35 q = n * d^{-1} mod 2^{qn * GMP_NUMB_BITS},
37 r = (n - q * d) * 2^{-qn * GMP_NUMB_BITS}
39 Stores q at qp. Stores the n least significant limbs of r at the high half
40 of np, and returns the borrow from the subtraction n - q*d.
42 d must be odd. dinv is (-d)^-1 mod 2^GMP_NUMB_BITS. */
45 mpn_dcpi1_bdiv_qr_n_itch (mp_size_t n)
51 mpn_dcpi1_bdiv_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
52 mp_limb_t dinv, mp_ptr tp)
58 lo = n >> 1; /* floor(n/2) */
59 hi = n - lo; /* ceil(n/2) */
61 if (BELOW_THRESHOLD (lo, DC_BDIV_QR_THRESHOLD))
62 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * lo, dp, lo, dinv);
64 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp);
66 mpn_mul (tp, dp + lo, hi, qp, lo);
68 mpn_incr_u (tp + lo, cy);
69 rh = mpn_sub (np + lo, np + lo, n + hi, tp, n);
71 if (BELOW_THRESHOLD (hi, DC_BDIV_QR_THRESHOLD))
72 cy = mpn_sbpi1_bdiv_qr (qp + lo, np + lo, 2 * hi, dp, hi, dinv);
74 cy = mpn_dcpi1_bdiv_qr_n (qp + lo, np + lo, dp, hi, dinv, tp);
76 mpn_mul (tp, qp + lo, hi, dp + hi, lo);
78 mpn_incr_u (tp + hi, cy);
79 rh += mpn_sub_n (np + n, np + n, tp, n);
85 mpn_dcpi1_bdiv_qr (mp_ptr qp, mp_ptr np, mp_size_t nn,
86 mp_srcptr dp, mp_size_t dn, mp_limb_t dinv)
95 ASSERT (dn >= 2); /* to adhere to mpn_sbpi1_div_qr's limits */
96 ASSERT (nn - dn >= 1); /* to adhere to mpn_sbpi1_div_qr's limits */
99 tp = TMP_SALLOC_LIMBS (dn);
105 /* Reduce qn mod dn without division, optimizing small operations. */
110 /* Perform the typically smaller block first. */
111 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
112 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
114 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
120 mpn_mul (tp, qp, qn, dp + qn, dn - qn);
122 mpn_mul (tp, dp + qn, dn - qn, qp, qn);
123 mpn_incr_u (tp + qn, cy);
125 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn);
135 rr += mpn_sub_1 (np + dn, np + dn, qn, cy);
136 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp);
146 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
147 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
149 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
155 mpn_mul (tp, qp, qn, dp + qn, dn - qn);
157 mpn_mul (tp, dp + qn, dn - qn, qp, qn);
158 mpn_incr_u (tp + qn, cy);
160 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn);