1 /* Mulders' short product, square and division. 2 3 Copyright 2005-2018 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 /* References: 24 [1] Short Division of Long Integers, David Harvey and Paul Zimmermann, 25 Proceedings of the 20th Symposium on Computer Arithmetic (ARITH-20), 26 July 25-27, 2011, pages 7-14. 27 */ 28 29 /* Defines it to 1 to use short div (or 0 for FoldDiv(K)) */ 30 #define USE_SHORT_DIV 1 31 32 #define MPFR_NEED_LONGLONG_H 33 #include "mpfr-impl.h" 34 35 #ifndef MUL_FFT_THRESHOLD 36 #define MUL_FFT_THRESHOLD 8448 37 #endif 38 39 /* Don't use MPFR_MULHIGH_SIZE since it is handled by tuneup */ 40 #ifdef MPFR_MULHIGH_TAB_SIZE 41 static short mulhigh_ktab[MPFR_MULHIGH_TAB_SIZE]; 42 #else 43 static short mulhigh_ktab[] = {MPFR_MULHIGH_TAB}; 44 #define MPFR_MULHIGH_TAB_SIZE (numberof_const (mulhigh_ktab)) 45 #endif 46 47 /* Put in rp[n..2n-1] an approximation of the n high limbs 48 of {up, n} * {vp, n}. The error is less than n ulps of rp[n] (and the 49 approximation is always less or equal to the truncated full product). 50 Assume 2n limbs are allocated at rp. 51 52 Implements Algorithm ShortMulNaive from [1]. 53 */ 54 static void 55 mpfr_mulhigh_n_basecase (mpfr_limb_ptr rp, mpfr_limb_srcptr up, 56 mpfr_limb_srcptr vp, mp_size_t n) 57 { 58 mp_size_t i; 59 60 rp += n - 1; 61 umul_ppmm (rp[1], rp[0], up[n-1], vp[0]); /* we neglect up[0..n-2]*vp[0], 62 which is less than B^n */ 63 for (i = 1 ; i < n ; i++) 64 /* here, we neglect up[0..n-i-2] * vp[i], which is less than B^n too */ 65 rp[i + 1] = mpn_addmul_1 (rp, up + (n - i - 1), i + 1, vp[i]); 66 /* in total, we neglect less than n*B^n, i.e., n ulps of rp[n]. */ 67 } 68 69 /* Put in rp[n..2n-1] an approximation of the n high limbs 70 of {np, n} * {mp, n}. The error is less than n ulps of rp[n] (and the 71 approximation is always less or equal to the truncated full product). 72 73 Implements Algorithm ShortMul from [1]. 74 */ 75 void 76 mpfr_mulhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp, 77 mp_size_t n) 78 { 79 mp_size_t k; 80 81 MPFR_STAT_STATIC_ASSERT (MPFR_MULHIGH_TAB_SIZE >= 8); /* so that 3*(n/4) > n/2 */ 82 k = MPFR_LIKELY (n < MPFR_MULHIGH_TAB_SIZE) ? mulhigh_ktab[n] : 3*(n/4); 83 /* Algorithm ShortMul from [1] requires k >= (n+3)/2, which translates 84 into k >= (n+4)/2 in the C language. */ 85 MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n)); 86 if (k < 0) 87 mpn_mul_basecase (rp, np, n, mp, n); /* result is exact, no error */ 88 else if (k == 0) 89 mpfr_mulhigh_n_basecase (rp, np, mp, n); /* basecase error < n ulps */ 90 else if (n > MUL_FFT_THRESHOLD) 91 mpn_mul_n (rp, np, mp, n); /* result is exact, no error */ 92 else 93 { 94 mp_size_t l = n - k; 95 mp_limb_t cy; 96 97 mpn_mul_n (rp + 2 * l, np + l, mp + l, k); /* fills rp[2l..2n-1] */ 98 mpfr_mulhigh_n (rp, np + k, mp, l); /* fills rp[l-1..2l-1] */ 99 cy = mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 100 mpfr_mulhigh_n (rp, np, mp + k, l); /* fills rp[l-1..2l-1] */ 101 cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 102 mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */ 103 } 104 } 105 106 #if USE_SHORT_DIV == 0 107 /* Put in rp[0..n] the n+1 low limbs of {up, n} * {vp, n}. 108 Assume 2n limbs are allocated at rp. */ 109 static void 110 mpfr_mullow_n_basecase (mpfr_limb_ptr rp, mpfr_limb_srcptr up, 111 mpfr_limb_srcptr vp, mp_size_t n) 112 { 113 mp_size_t i; 114 115 rp[n] = mpn_mul_1 (rp, up, n, vp[0]); 116 for (i = 1 ; i < n ; i++) 117 mpn_addmul_1 (rp + i, up, n - i + 1, vp[i]); 118 } 119 120 /* Put in rp[0..n] the n+1 low limbs of {np, n} * {mp, n}. 121 Assume 2n limbs are allocated at rp. */ 122 void 123 mpfr_mullow_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp, 124 mp_size_t n) 125 { 126 mp_size_t k; 127 128 MPFR_STAT_STATIC_ASSERT (MPFR_MULHIGH_TAB_SIZE >= 8); /* so that 3*(n/4) > n/2 */ 129 k = MPFR_LIKELY (n < MPFR_MULHIGH_TAB_SIZE) ? mulhigh_ktab[n] : 3*(n/4); 130 MPFR_ASSERTD (k == -1 || k == 0 || (2 * k >= n && k < n)); 131 if (k < 0) 132 mpn_mul_basecase (rp, np, n, mp, n); 133 else if (k == 0) 134 mpfr_mullow_n_basecase (rp, np, mp, n); 135 else if (n > MUL_FFT_THRESHOLD) 136 mpn_mul_n (rp, np, mp, n); 137 else 138 { 139 mp_size_t l = n - k; 140 141 mpn_mul_n (rp, np, mp, k); /* fills rp[0..2k] */ 142 mpfr_mullow_n (rp + n, np + k, mp, l); /* fills rp[n..n+2l] */ 143 mpn_add_n (rp + k, rp + k, rp + n, l + 1); 144 mpfr_mullow_n (rp + n, np, mp + k, l); /* fills rp[n..n+2l] */ 145 mpn_add_n (rp + k, rp + k, rp + n, l + 1); 146 } 147 } 148 #endif 149 150 #ifdef MPFR_SQRHIGH_TAB_SIZE 151 static short sqrhigh_ktab[MPFR_SQRHIGH_TAB_SIZE]; 152 #else 153 static short sqrhigh_ktab[] = {MPFR_SQRHIGH_TAB}; 154 #define MPFR_SQRHIGH_TAB_SIZE (numberof_const (sqrhigh_ktab)) 155 #endif 156 157 /* Put in rp[n..2n-1] an approximation of the n high limbs 158 of {np, n}^2. The error is less than n ulps of rp[n]. */ 159 void 160 mpfr_sqrhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mp_size_t n) 161 { 162 mp_size_t k; 163 164 MPFR_STAT_STATIC_ASSERT (MPFR_SQRHIGH_TAB_SIZE > 2); /* ensures k < n */ 165 k = MPFR_LIKELY (n < MPFR_SQRHIGH_TAB_SIZE) ? sqrhigh_ktab[n] 166 : (n+4)/2; /* ensures that k >= (n+3)/2 */ 167 MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n)); 168 if (k < 0) 169 /* we can't use mpn_sqr_basecase here, since it requires 170 n <= SQR_KARATSUBA_THRESHOLD, where SQR_KARATSUBA_THRESHOLD 171 is not exported by GMP */ 172 mpn_sqr (rp, np, n); 173 else if (k == 0) 174 mpfr_mulhigh_n_basecase (rp, np, np, n); 175 else 176 { 177 mp_size_t l = n - k; 178 mp_limb_t cy; 179 180 mpn_sqr (rp + 2 * l, np + l, k); /* fills rp[2l..2n-1] */ 181 mpfr_mulhigh_n (rp, np, np + k, l); /* fills rp[l-1..2l-1] */ 182 /* {rp+n-1,l+1} += 2 * {rp+l-1,l+1} */ 183 cy = mpn_lshift (rp + l - 1, rp + l - 1, l + 1, 1); 184 cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 185 mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */ 186 } 187 } 188 189 #if USE_SHORT_DIV == 1 190 191 #ifdef MPFR_DIVHIGH_TAB_SIZE 192 static short divhigh_ktab[MPFR_DIVHIGH_TAB_SIZE]; 193 #else 194 static short divhigh_ktab[] = {MPFR_DIVHIGH_TAB}; 195 #define MPFR_DIVHIGH_TAB_SIZE (numberof_const (divhigh_ktab)) 196 #endif 197 198 #if !(defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_SBPI1_DIVAPPR_Q)) 199 /* Put in Q={qp, n} an approximation of N={np, 2*n} divided by D={dp, n}, 200 with the most significant limb of the quotient as return value (0 or 1). 201 Assumes the most significant bit of D is set. Clobbers N. 202 203 The approximate quotient Q satisfies - 2(n-1) < N/D - Q <= 4. 204 */ 205 static mp_limb_t 206 mpfr_divhigh_n_basecase (mpfr_limb_ptr qp, mpfr_limb_ptr np, 207 mpfr_limb_srcptr dp, mp_size_t n) 208 { 209 mp_limb_t qh, d1, d0, dinv, q2, q1, q0; 210 mpfr_pi1_t dinv2; 211 212 np += n; 213 214 if ((qh = (mpn_cmp (np, dp, n) >= 0))) 215 mpn_sub_n (np, np, dp, n); 216 217 /* now {np, n} is less than D={dp, n}, which implies np[n-1] <= dp[n-1] */ 218 219 d1 = dp[n - 1]; 220 221 if (n == 1) 222 { 223 invert_limb (dinv, d1); 224 umul_ppmm (q1, q0, np[0], dinv); 225 qp[0] = np[0] + q1; 226 return qh; 227 } 228 229 /* now n >= 2 */ 230 d0 = dp[n - 2]; 231 invert_pi1 (dinv2, d1, d0); 232 /* dinv2.inv32 = floor ((B^3 - 1) / (d0 + d1 B)) - B */ 233 while (n > 1) 234 { 235 /* Invariant: it remains to reduce n limbs from N (in addition to the 236 initial low n limbs). 237 Since n >= 2 here, necessarily we had n >= 2 initially, which means 238 that in addition to the limb np[n-1] to reduce, we have at least 2 239 extra limbs, thus accessing np[n-3] is valid. */ 240 241 /* Warning: we can have np[n-1]>d1 or (np[n-1]=d1 and np[n-2]>=d0) here, 242 since we truncate the divisor at each step, but since {np,n} < D 243 originally, the largest possible partial quotient is B-1. */ 244 if (MPFR_UNLIKELY(np[n-1] > d1 || (np[n-1] == d1 && np[n-2] >= d0))) 245 q2 = MPFR_LIMB_MAX; 246 else 247 udiv_qr_3by2 (q2, q1, q0, np[n - 1], np[n - 2], np[n - 3], 248 d1, d0, dinv2.inv32); 249 /* since q2 = floor((np[n-1]*B^2+np[n-2]*B+np[n-3])/(d1*B+d0)), 250 we have q2 <= (np[n-1]*B^2+np[n-2]*B+np[n-3])/(d1*B+d0), 251 thus np[n-1]*B^2+np[n-2]*B+np[n-3] >= q2*(d1*B+d0) 252 and {np-1, n} >= q2*D - q2*B^(n-2) >= q2*D - B^(n-1) 253 thus {np-1, n} - (q2-1)*D >= D - B^(n-1) >= 0 254 which proves that at most one correction is needed */ 255 q0 = mpn_submul_1 (np - 1, dp, n, q2); 256 if (MPFR_UNLIKELY(q0 > np[n - 1])) 257 { 258 mpn_add_n (np - 1, np - 1, dp, n); 259 q2 --; 260 } 261 qp[--n] = q2; 262 dp ++; 263 } 264 265 /* we have B+dinv2 = floor((B^3-1)/(d1*B+d0)) < B^2/d1 266 q1 = floor(np[0]*(B+dinv2)/B) <= floor(np[0]*B/d1) 267 <= floor((np[0]*B+np[1])/d1) 268 thus q1 is not larger than the true quotient. 269 q1 > np[0]*(B+dinv2)/B - 1 > np[0]*(B^3-1)/(d1*B+d0)/B - 2 270 For d1*B+d0 <> B^2/2, we have B+dinv2 = floor(B^3/(d1*B+d0)) 271 thus q1 > np[0]*B^2/(d1*B+d0) - 2, i.e., 272 (d1*B+d0)*q1 > np[0]*B^2 - 2*(d1*B+d0) 273 d1*B*q1 > np[0]*B^2 - 2*d1*B - 2*d0 - d0*q1 >= np[0]*B^2 - 2*d1*B - B^2 274 thus q1 > np[0]*B/d1 - 2 - B/d1 > np[0]*B/d1 - 4. 275 276 For d1*B+d0 = B^2/2, dinv2 = B-1 thus q1 > np[0]*(2B-1)/B - 1 > 277 np[0]*B/d1 - 2. 278 279 In all cases, if q = floor((np[0]*B+np[1])/d1), we have: 280 q - 4 <= q1 <= q 281 */ 282 umul_ppmm (q1, q0, np[0], dinv2.inv32); 283 qp[0] = np[0] + q1; 284 285 return qh; 286 } 287 #endif 288 289 /* Put in {qp, n} an approximation of N={np, 2*n} divided by D={dp, n}, 290 with the most significant limb of the quotient as return value (0 or 1). 291 Assumes the most significant bit of D is set. Clobbers N. 292 293 This implements the ShortDiv algorithm from reference [1]. 294 */ 295 mp_limb_t 296 mpfr_divhigh_n (mpfr_limb_ptr qp, mpfr_limb_ptr np, mpfr_limb_ptr dp, 297 mp_size_t n) 298 { 299 mp_size_t k, l; 300 mp_limb_t qh, cy; 301 mpfr_limb_ptr tp; 302 MPFR_TMP_DECL(marker); 303 304 MPFR_STAT_STATIC_ASSERT (MPFR_DIVHIGH_TAB_SIZE >= 15); /* so that 2*(n/3) >= (n+4)/2 */ 305 k = MPFR_LIKELY (n < MPFR_DIVHIGH_TAB_SIZE) ? divhigh_ktab[n] : 2*(n/3); 306 307 if (k == 0) 308 #if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_SBPI1_DIVAPPR_Q) 309 { 310 mpfr_pi1_t dinv2; 311 invert_pi1 (dinv2, dp[n - 1], dp[n - 2]); 312 return __gmpn_sbpi1_divappr_q (qp, np, n + n, dp, n, dinv2.inv32); 313 } 314 #else /* use our own code for base-case short division */ 315 return mpfr_divhigh_n_basecase (qp, np, dp, n); 316 #endif 317 else if (k == n) 318 /* for k=n, we use a division with remainder (mpn_divrem), 319 which computes the exact quotient */ 320 return mpn_divrem (qp, 0, np, 2 * n, dp, n); 321 322 MPFR_ASSERTD ((n+4)/2 <= k && k < n); /* bounds from [1] */ 323 MPFR_TMP_MARK (marker); 324 l = n - k; 325 /* first divide the most significant 2k limbs from N by the most significant 326 k limbs of D */ 327 qh = mpn_divrem (qp + l, 0, np + 2 * l, 2 * k, dp + l, k); /* exact */ 328 329 /* it remains {np,2l+k} = {np,n+l} as remainder */ 330 331 /* now we have to subtract high(Q1)*D0 where Q1=qh*B^k+{qp+l,k} and 332 D0={dp,l} */ 333 tp = MPFR_TMP_LIMBS_ALLOC (2 * l); 334 mpfr_mulhigh_n (tp, qp + k, dp, l); 335 /* we are only interested in the upper l limbs from {tp,2l} */ 336 cy = mpn_sub_n (np + n, np + n, tp + l, l); 337 if (qh) 338 cy += mpn_sub_n (np + n, np + n, dp, l); 339 while (cy > 0) /* Q1 was too large: subtract 1 to Q1 and add D to np+l */ 340 { 341 qh -= mpn_sub_1 (qp + l, qp + l, k, MPFR_LIMB_ONE); 342 cy -= mpn_add_n (np + l, np + l, dp, n); 343 } 344 345 /* now it remains {np,n+l} to divide by D */ 346 cy = mpfr_divhigh_n (qp, np + k, dp + k, l); 347 qh += mpn_add_1 (qp + l, qp + l, k, cy); 348 MPFR_TMP_FREE(marker); 349 350 return qh; 351 } 352 353 #else /* below is the FoldDiv(K) algorithm from [1] */ 354 355 mp_limb_t 356 mpfr_divhigh_n (mpfr_limb_ptr qp, mpfr_limb_ptr np, mpfr_limb_ptr dp, 357 mp_size_t n) 358 { 359 mp_size_t k, r; 360 mpfr_limb_ptr ip, tp, up; 361 mp_limb_t qh = 0, cy, cc; 362 int count; 363 MPFR_TMP_DECL(marker); 364 365 #define K 3 366 if (n < K) 367 return mpn_divrem (qp, 0, np, 2 * n, dp, n); 368 369 k = (n - 1) / K + 1; /* ceil(n/K) */ 370 371 MPFR_TMP_MARK (marker); 372 ip = MPFR_TMP_LIMBS_ALLOC (k + 1); 373 tp = MPFR_TMP_LIMBS_ALLOC (n + k); 374 up = MPFR_TMP_LIMBS_ALLOC (2 * (k + 1)); 375 mpn_invert (ip, dp + n - (k + 1), k + 1, NULL); /* takes about 13% for n=1000 */ 376 /* {ip, k+1} = floor((B^(2k+2)-1)/D - B^(k+1) where D = {dp+n-(k+1),k+1} */ 377 for (r = n, cc = 0UL; r > 0;) 378 { 379 /* cc is the carry at np[n+r] */ 380 MPFR_ASSERTD(cc <= 1); 381 /* FIXME: why can we have cc as large as say 8? */ 382 count = 0; 383 while (cc > 0) 384 { 385 count ++; 386 MPFR_ASSERTD(count <= 1); 387 /* subtract {dp+n-r,r} from {np+n,r} */ 388 cc -= mpn_sub_n (np + n, np + n, dp + n - r, r); 389 /* add 1 at qp[r] */ 390 qh += mpn_add_1 (qp + r, qp + r, n - r, 1UL); 391 } 392 /* it remains r limbs to reduce, i.e., the remainder is {np, n+r} */ 393 if (r < k) 394 { 395 ip += k - r; 396 k = r; 397 } 398 /* now r >= k */ 399 /* qp + r - 2 * k -> up */ 400 mpfr_mulhigh_n (up, np + n + r - (k + 1), ip, k + 1); 401 /* take into account the term B^k in the inverse: B^k * {np+n+r-k, k} */ 402 cy = mpn_add_n (qp + r - k, up + k + 2, np + n + r - k, k); 403 /* since we need only r limbs of tp (below), it suffices to consider 404 r high limbs of dp */ 405 if (r > k) 406 { 407 #if 0 408 mpn_mul (tp, dp + n - r, r, qp + r - k, k); 409 #else /* use a short product for the low k x k limbs */ 410 /* we know the upper k limbs of the r-limb product cancel with the 411 remainder, thus we only need to compute the low r-k limbs */ 412 if (r - k >= k) 413 mpn_mul (tp + k, dp + n - r + k, r - k, qp + r - k, k); 414 else /* r-k < k */ 415 { 416 /* #define LOW */ 417 #ifndef LOW 418 mpn_mul (tp + k, qp + r - k, k, dp + n - r + k, r - k); 419 #else 420 mpfr_mullow_n_basecase (tp + k, qp + r - k, dp + n - r + k, r - k); 421 /* take into account qp[2r-2k] * dp[n - r + k] */ 422 tp[r] += qp[2*r-2*k] * dp[n - r + k]; 423 #endif 424 /* tp[k..r] is filled */ 425 } 426 #if 0 427 mpfr_mulhigh_n (up, dp + n - r, qp + r - k, k); 428 #else /* compute one more limb. FIXME: we could add one limb of dp in the 429 above, to save one mpn_addmul_1 call */ 430 mpfr_mulhigh_n (up, dp + n - r, qp + r - k, k - 1); /* {up,2k-2} */ 431 /* add {qp + r - k, k - 1} * dp[n-r+k-1] */ 432 up[2*k-2] = mpn_addmul_1 (up + k - 1, qp + r - k, k-1, dp[n-r+k-1]); 433 /* add {dp+n-r, k} * qp[r-1] */ 434 up[2*k-1] = mpn_addmul_1 (up + k - 1, dp + n - r, k, qp[r-1]); 435 #endif 436 #ifndef LOW 437 cc = mpn_add_n (tp + k, tp + k, up + k, k); 438 mpn_add_1 (tp + 2 * k, tp + 2 * k, r - k, cc); 439 #else 440 /* update tp[k..r] */ 441 if (r - k + 1 <= k) 442 mpn_add_n (tp + k, tp + k, up + k, r - k + 1); 443 else /* r - k >= k */ 444 { 445 cc = mpn_add_n (tp + k, tp + k, up + k, k); 446 mpn_add_1 (tp + 2 * k, tp + 2 * k, r - 2 * k + 1, cc); 447 } 448 #endif 449 #endif 450 } 451 else /* last step: since we only want the quotient, no need to update, 452 just propagate the carry cy */ 453 { 454 MPFR_ASSERTD(r < n); 455 if (cy > 0) 456 qh += mpn_add_1 (qp + r, qp + r, n - r, cy); 457 break; 458 } 459 /* subtract {tp, n+k} from {np+r-k, n+k}; however we only want to 460 update {np+n, n} */ 461 /* we should have tp[r] = np[n+r-k] up to 1 */ 462 MPFR_ASSERTD(tp[r] == np[n + r - k] || tp[r] + 1 == np[n + r - k]); 463 #ifndef LOW 464 cc = mpn_sub_n (np + n - 1, np + n - 1, tp + k - 1, r + 1); /* borrow at np[n+r] */ 465 #else 466 cc = mpn_sub_n (np + n - 1, np + n - 1, tp + k - 1, r - k + 2); 467 #endif 468 /* if cy = 1, subtract {dp, n} from {np+r, n}, thus 469 {dp+n-r,r} from {np+n,r} */ 470 if (cy) 471 { 472 if (r < n) 473 cc += mpn_sub_n (np + n - 1, np + n - 1, dp + n - r - 1, r + 1); 474 else 475 cc += mpn_sub_n (np + n, np + n, dp + n - r, r); 476 /* propagate cy */ 477 if (r == n) 478 qh = cy; 479 else 480 qh += mpn_add_1 (qp + r, qp + r, n - r, cy); 481 } 482 /* cc is the borrow at np[n+r] */ 483 count = 0; 484 while (cc > 0) /* quotient was too large */ 485 { 486 count++; 487 MPFR_ASSERTD (count <= 1); 488 cy = mpn_add_n (np + n, np + n, dp + n - (r - k), r - k); 489 cc -= mpn_add_1 (np + n + r - k, np + n + r - k, k, cy); 490 qh -= mpn_sub_1 (qp + r - k, qp + r - k, n - (r - k), 1UL); 491 } 492 r -= k; 493 cc = np[n + r]; 494 } 495 MPFR_TMP_FREE(marker); 496 497 return qh; 498 } 499 #endif 500