1 /* mpn_sbpi1_div_q -- Schoolbook division using the M�ller-Granlund 3/2 2 division algorithm. 3 4 Contributed to the GNU project by Torbjorn Granlund. 5 6 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY 7 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 8 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. 9 10 Copyright 2007, 2009 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 28 #include "gmp.h" 29 #include "gmp-impl.h" 30 #include "longlong.h" 31 32 mp_limb_t 33 mpn_sbpi1_div_q (mp_ptr qp, 34 mp_ptr np, mp_size_t nn, 35 mp_srcptr dp, mp_size_t dn, 36 mp_limb_t dinv) 37 { 38 mp_limb_t qh; 39 mp_size_t qn, i; 40 mp_limb_t n1, n0; 41 mp_limb_t d1, d0; 42 mp_limb_t cy, cy1; 43 mp_limb_t q; 44 mp_limb_t flag; 45 46 mp_size_t dn_orig = dn; 47 mp_srcptr dp_orig = dp; 48 mp_ptr np_orig = np; 49 50 ASSERT (dn > 2); 51 ASSERT (nn >= dn); 52 ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0); 53 54 np += nn; 55 56 qn = nn - dn; 57 if (qn + 1 < dn) 58 { 59 dp += dn - (qn + 1); 60 dn = qn + 1; 61 } 62 63 qh = mpn_cmp (np - dn, dp, dn) >= 0; 64 if (qh != 0) 65 mpn_sub_n (np - dn, np - dn, dp, dn); 66 67 qp += qn; 68 69 dn -= 2; /* offset dn by 2 for main division loops, 70 saving two iterations in mpn_submul_1. */ 71 d1 = dp[dn + 1]; 72 d0 = dp[dn + 0]; 73 74 np -= 2; 75 76 n1 = np[1]; 77 78 for (i = qn - (dn + 2); i >= 0; i--) 79 { 80 np--; 81 if (UNLIKELY (n1 == d1) && np[1] == d0) 82 { 83 q = GMP_NUMB_MASK; 84 mpn_submul_1 (np - dn, dp, dn + 2, q); 85 n1 = np[1]; /* update n1, last loop's value will now be invalid */ 86 } 87 else 88 { 89 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 90 91 cy = mpn_submul_1 (np - dn, dp, dn, q); 92 93 cy1 = n0 < cy; 94 n0 = (n0 - cy) & GMP_NUMB_MASK; 95 cy = n1 < cy1; 96 n1 -= cy1; 97 np[0] = n0; 98 99 if (UNLIKELY (cy != 0)) 100 { 101 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1); 102 q--; 103 } 104 } 105 106 *--qp = q; 107 } 108 109 flag = ~CNST_LIMB(0); 110 111 if (dn >= 0) 112 { 113 for (i = dn; i > 0; i--) 114 { 115 np--; 116 if (UNLIKELY (n1 >= (d1 & flag))) 117 { 118 q = GMP_NUMB_MASK; 119 cy = mpn_submul_1 (np - dn, dp, dn + 2, q); 120 121 if (UNLIKELY (n1 != cy)) 122 { 123 if (n1 < (cy & flag)) 124 { 125 q--; 126 mpn_add_n (np - dn, np - dn, dp, dn + 2); 127 } 128 else 129 flag = 0; 130 } 131 n1 = np[1]; 132 } 133 else 134 { 135 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 136 137 cy = mpn_submul_1 (np - dn, dp, dn, q); 138 139 cy1 = n0 < cy; 140 n0 = (n0 - cy) & GMP_NUMB_MASK; 141 cy = n1 < cy1; 142 n1 -= cy1; 143 np[0] = n0; 144 145 if (UNLIKELY (cy != 0)) 146 { 147 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1); 148 q--; 149 } 150 } 151 152 *--qp = q; 153 154 /* Truncate operands. */ 155 dn--; 156 dp++; 157 } 158 159 np--; 160 if (UNLIKELY (n1 >= (d1 & flag))) 161 { 162 q = GMP_NUMB_MASK; 163 cy = mpn_submul_1 (np, dp, 2, q); 164 165 if (UNLIKELY (n1 != cy)) 166 { 167 if (n1 < (cy & flag)) 168 { 169 q--; 170 add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]); 171 } 172 else 173 flag = 0; 174 } 175 n1 = np[1]; 176 } 177 else 178 { 179 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 180 181 np[0] = n0; 182 np[1] = n1; 183 } 184 185 *--qp = q; 186 } 187 ASSERT_ALWAYS (np[1] == n1); 188 np += 2; 189 190 191 dn = dn_orig; 192 if (UNLIKELY (n1 < (dn & flag))) 193 { 194 mp_limb_t q, x; 195 196 /* The quotient may be too large if the remainder is small. Recompute 197 for above ignored operand parts, until the remainder spills. 198 199 FIXME: The quality of this code isn't the same as the code above. 200 1. We don't compute things in an optimal order, high-to-low, in order 201 to terminate as quickly as possible. 202 2. We mess with pointers and sizes, adding and subtracting and 203 adjusting to get things right. It surely could be streamlined. 204 3. The only termination criteria are that we determine that the 205 quotient needs to be adjusted, or that we have recomputed 206 everything. We should stop when the remainder is so large 207 that no additional subtracting could make it spill. 208 4. If nothing else, we should not do two loops of submul_1 over the 209 data, instead handle both the triangularization and chopping at 210 once. */ 211 212 x = n1; 213 214 if (dn > 2) 215 { 216 /* Compensate for triangularization. */ 217 mp_limb_t y; 218 219 dp = dp_orig; 220 if (qn + 1 < dn) 221 { 222 dp += dn - (qn + 1); 223 dn = qn + 1; 224 } 225 226 y = np[-2]; 227 228 for (i = dn - 3; i >= 0; i--) 229 { 230 q = qp[i]; 231 cy = mpn_submul_1 (np - (dn - i), dp, dn - i - 2, q); 232 233 if (y < cy) 234 { 235 if (x == 0) 236 { 237 cy = mpn_sub_1 (qp, qp, qn, 1); 238 ASSERT_ALWAYS (cy == 0); 239 return qh - cy; 240 } 241 x--; 242 } 243 y -= cy; 244 } 245 np[-2] = y; 246 } 247 248 dn = dn_orig; 249 if (qn + 1 < dn) 250 { 251 /* Compensate for ignored dividend and divisor tails. */ 252 253 dp = dp_orig; 254 np = np_orig; 255 256 if (qh != 0) 257 { 258 cy = mpn_sub_n (np + qn, np + qn, dp, dn - (qn + 1)); 259 if (cy != 0) 260 { 261 if (x == 0) 262 { 263 if (qn != 0) 264 cy = mpn_sub_1 (qp, qp, qn, 1); 265 return qh - cy; 266 } 267 x--; 268 } 269 } 270 271 if (qn == 0) 272 return qh; 273 274 for (i = dn - qn - 2; i >= 0; i--) 275 { 276 cy = mpn_submul_1 (np + i, qp, qn, dp[i]); 277 cy = mpn_sub_1 (np + qn + i, np + qn + i, dn - qn - i - 1, cy); 278 if (cy != 0) 279 { 280 if (x == 0) 281 { 282 cy = mpn_sub_1 (qp, qp, qn, 1); 283 return qh; 284 } 285 x--; 286 } 287 } 288 } 289 } 290 291 return qh; 292 } 293