1 /* mpn_dcpi1_bdiv_qr -- divide-and-conquer Hensel division with precomputed 2 inverse, returning quotient and remainder. 3 4 Contributed to the GNU project by Niels M�ller and Torbjorn Granlund. 5 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. 9 10 Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc. 11 12 This file is part of the GNU MP Library. 13 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. 18 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. 23 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/. */ 26 27 #include "gmp.h" 28 #include "gmp-impl.h" 29 30 31 /* Computes Hensel binary division of {np, 2*n} by {dp, n}. 32 33 Output: 34 35 q = n * d^{-1} mod 2^{qn * GMP_NUMB_BITS}, 36 37 r = (n - q * d) * 2^{-qn * GMP_NUMB_BITS} 38 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. 41 42 d must be odd. dinv is (-d)^-1 mod 2^GMP_NUMB_BITS. */ 43 44 mp_size_t 45 mpn_dcpi1_bdiv_qr_n_itch (mp_size_t n) 46 { 47 return n; 48 } 49 50 mp_limb_t 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) 53 { 54 mp_size_t lo, hi; 55 mp_limb_t cy; 56 mp_limb_t rh; 57 58 lo = n >> 1; /* floor(n/2) */ 59 hi = n - lo; /* ceil(n/2) */ 60 61 if (BELOW_THRESHOLD (lo, DC_BDIV_QR_THRESHOLD)) 62 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * lo, dp, lo, dinv); 63 else 64 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp); 65 66 mpn_mul (tp, dp + lo, hi, qp, lo); 67 68 mpn_incr_u (tp + lo, cy); 69 rh = mpn_sub (np + lo, np + lo, n + hi, tp, n); 70 71 if (BELOW_THRESHOLD (hi, DC_BDIV_QR_THRESHOLD)) 72 cy = mpn_sbpi1_bdiv_qr (qp + lo, np + lo, 2 * hi, dp, hi, dinv); 73 else 74 cy = mpn_dcpi1_bdiv_qr_n (qp + lo, np + lo, dp, hi, dinv, tp); 75 76 mpn_mul (tp, qp + lo, hi, dp + hi, lo); 77 78 mpn_incr_u (tp + hi, cy); 79 rh += mpn_sub_n (np + n, np + n, tp, n); 80 81 return rh; 82 } 83 84 mp_limb_t 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) 87 { 88 mp_size_t qn; 89 mp_limb_t rr, cy; 90 mp_ptr tp; 91 TMP_DECL; 92 93 TMP_MARK; 94 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 */ 97 ASSERT (dp[0] & 1); 98 99 tp = TMP_SALLOC_LIMBS (dn); 100 101 qn = nn - dn; 102 103 if (qn > dn) 104 { 105 /* Reduce qn mod dn without division, optimizing small operations. */ 106 do 107 qn -= dn; 108 while (qn > dn); 109 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); 113 else 114 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp); 115 116 rr = 0; 117 if (qn != dn) 118 { 119 if (qn > dn - qn) 120 mpn_mul (tp, qp, qn, dp + qn, dn - qn); 121 else 122 mpn_mul (tp, dp + qn, dn - qn, qp, qn); 123 mpn_incr_u (tp + qn, cy); 124 125 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn); 126 cy = 0; 127 } 128 129 np += qn; 130 qp += qn; 131 132 qn = nn - dn - qn; 133 do 134 { 135 rr += mpn_sub_1 (np + dn, np + dn, qn, cy); 136 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp); 137 qp += dn; 138 np += dn; 139 qn -= dn; 140 } 141 while (qn > 0); 142 TMP_FREE; 143 return rr + cy; 144 } 145 146 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD)) 147 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv); 148 else 149 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp); 150 151 rr = 0; 152 if (qn != dn) 153 { 154 if (qn > dn - qn) 155 mpn_mul (tp, qp, qn, dp + qn, dn - qn); 156 else 157 mpn_mul (tp, dp + qn, dn - qn, qp, qn); 158 mpn_incr_u (tp + qn, cy); 159 160 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn); 161 cy = 0; 162 } 163 164 TMP_FREE; 165 return rr + cy; 166 } 167