1 /* Mulders' short product, square and division. 2 3 Copyright 2005-2023 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 https://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 #define MPFR_NEED_LONGLONG_H 30 #include "mpfr-impl.h" 31 32 /* Don't use MPFR_MULHIGH_SIZE since it is handled by tuneup */ 33 #ifdef MPFR_MULHIGH_TAB_SIZE 34 static short mulhigh_ktab[MPFR_MULHIGH_TAB_SIZE]; 35 #else 36 static short mulhigh_ktab[] = {MPFR_MULHIGH_TAB}; 37 #define MPFR_MULHIGH_TAB_SIZE (numberof_const (mulhigh_ktab)) 38 #endif 39 40 /* Put in rp[n..2n-1] an approximation of the n high limbs 41 of {up, n} * {vp, n}. The error is less than n ulps of rp[n] (and the 42 approximation is always less or equal to the truncated full product). 43 Assume 2n limbs are allocated at rp. 44 45 Implements Algorithm ShortMulNaive from [1]. 46 */ 47 static void 48 mpfr_mulhigh_n_basecase (mpfr_limb_ptr rp, mpfr_limb_srcptr up, 49 mpfr_limb_srcptr vp, mp_size_t n) 50 { 51 mp_size_t i; 52 53 rp += n - 1; 54 umul_ppmm (rp[1], rp[0], up[n-1], vp[0]); /* we neglect up[0..n-2]*vp[0], 55 which is less than B^n */ 56 for (i = 1 ; i < n ; i++) 57 /* here, we neglect up[0..n-i-2] * vp[i], which is less than B^n too */ 58 rp[i + 1] = mpn_addmul_1 (rp, up + (n - i - 1), i + 1, vp[i]); 59 /* in total, we neglect less than n*B^n, i.e., n ulps of rp[n]. */ 60 } 61 62 /* Put in rp[n..2n-1] an approximation of the n high limbs 63 of {np, n} * {mp, n}. The error is less than n ulps of rp[n] (and the 64 approximation is always less or equal to the truncated full product). 65 66 Implements Algorithm ShortMul from [1]. 67 */ 68 void 69 mpfr_mulhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mpfr_limb_srcptr mp, 70 mp_size_t n) 71 { 72 mp_size_t k; 73 74 MPFR_STAT_STATIC_ASSERT (MPFR_MULHIGH_TAB_SIZE >= 8); /* so that 3*(n/4) > n/2 */ 75 k = MPFR_LIKELY (n < MPFR_MULHIGH_TAB_SIZE) ? mulhigh_ktab[n] : 3*(n/4); 76 /* Algorithm ShortMul from [1] requires k >= (n+3)/2, which translates 77 into k >= (n+4)/2 in the C language. */ 78 MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n)); 79 if (k < 0) 80 mpn_mul_basecase (rp, np, n, mp, n); /* result is exact, no error */ 81 else if (k == 0) 82 mpfr_mulhigh_n_basecase (rp, np, mp, n); /* basecase error < n ulps */ 83 else if (n > MUL_FFT_THRESHOLD) 84 mpn_mul_n (rp, np, mp, n); /* result is exact, no error */ 85 else 86 { 87 mp_size_t l = n - k; 88 mp_limb_t cy; 89 90 mpn_mul_n (rp + 2 * l, np + l, mp + l, k); /* fills rp[2l..2n-1] */ 91 mpfr_mulhigh_n (rp, np + k, mp, l); /* fills rp[l-1..2l-1] */ 92 cy = mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 93 mpfr_mulhigh_n (rp, np, mp + k, l); /* fills rp[l-1..2l-1] */ 94 cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 95 mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */ 96 } 97 } 98 99 #ifdef MPFR_SQRHIGH_TAB_SIZE 100 static short sqrhigh_ktab[MPFR_SQRHIGH_TAB_SIZE]; 101 #else 102 static short sqrhigh_ktab[] = {MPFR_SQRHIGH_TAB}; 103 #define MPFR_SQRHIGH_TAB_SIZE (numberof_const (sqrhigh_ktab)) 104 #endif 105 106 /* Put in rp[n..2n-1] an approximation of the n high limbs 107 of {np, n}^2. The error is less than n ulps of rp[n]. */ 108 void 109 mpfr_sqrhigh_n (mpfr_limb_ptr rp, mpfr_limb_srcptr np, mp_size_t n) 110 { 111 mp_size_t k; 112 113 MPFR_STAT_STATIC_ASSERT (MPFR_SQRHIGH_TAB_SIZE > 2); /* ensures k < n */ 114 k = MPFR_LIKELY (n < MPFR_SQRHIGH_TAB_SIZE) ? sqrhigh_ktab[n] 115 : (n+4)/2; /* ensures that k >= (n+3)/2 */ 116 MPFR_ASSERTD (k == -1 || k == 0 || (k >= (n+4)/2 && k < n)); 117 if (k < 0) 118 /* we can't use mpn_sqr_basecase here, since it requires 119 n <= SQR_KARATSUBA_THRESHOLD, where SQR_KARATSUBA_THRESHOLD 120 is not exported by GMP */ 121 mpn_sqr (rp, np, n); 122 else if (k == 0) 123 mpfr_mulhigh_n_basecase (rp, np, np, n); 124 else 125 { 126 mp_size_t l = n - k; 127 mp_limb_t cy; 128 129 mpn_sqr (rp + 2 * l, np + l, k); /* fills rp[2l..2n-1] */ 130 mpfr_mulhigh_n (rp, np, np + k, l); /* fills rp[l-1..2l-1] */ 131 /* {rp+n-1,l+1} += 2 * {rp+l-1,l+1} */ 132 cy = mpn_lshift (rp + l - 1, rp + l - 1, l + 1, 1); 133 cy += mpn_add_n (rp + n - 1, rp + n - 1, rp + l - 1, l + 1); 134 mpn_add_1 (rp + n + l, rp + n + l, k, cy); /* propagate carry */ 135 } 136 } 137 138 #ifdef MPFR_DIVHIGH_TAB_SIZE 139 static short divhigh_ktab[MPFR_DIVHIGH_TAB_SIZE]; 140 #else 141 static short divhigh_ktab[] = {MPFR_DIVHIGH_TAB}; 142 #define MPFR_DIVHIGH_TAB_SIZE (numberof_const (divhigh_ktab)) 143 #endif 144 145 /* Put in Q={qp, n} an approximation of N={np, 2*n} divided by D={dp, n}, 146 with the most significant limb of the quotient as return value (0 or 1). 147 Assumes the most significant bit of D is set. Clobbers N. 148 149 The approximate quotient Q satisfies - 2(n-1) < N/D - Q <= 4. 150 151 Assumes n >= 2. 152 */ 153 static mp_limb_t 154 mpfr_divhigh_n_basecase (mpfr_limb_ptr qp, mpfr_limb_ptr np, 155 mpfr_limb_srcptr dp, mp_size_t n) 156 { 157 mp_limb_t qh, d1, d0, q2, q1, q0; 158 mpfr_pi1_t dinv2; 159 160 MPFR_ASSERTD(n >= 2); 161 162 np += n; 163 164 if ((qh = (mpn_cmp (np, dp, n) >= 0))) 165 mpn_sub_n (np, np, dp, n); 166 167 /* now {np, n} is less than D={dp, n}, which implies np[n-1] <= dp[n-1] */ 168 169 d1 = dp[n - 1]; 170 171 /* we assumed n >= 2 */ 172 d0 = dp[n - 2]; 173 invert_pi1 (dinv2, d1, d0); 174 /* dinv2.inv32 = floor ((B^3 - 1) / (d0 + d1 B)) - B */ 175 while (n > 1) 176 { 177 /* Invariant: it remains to reduce n limbs from N (in addition to the 178 initial low n limbs). 179 Since n >= 2 here, necessarily we had n >= 2 initially, which means 180 that in addition to the limb np[n-1] to reduce, we have at least 2 181 extra limbs, thus accessing np[n-3] is valid. */ 182 183 /* Warning: we can have np[n-1]>d1 or (np[n-1]=d1 and np[n-2]>=d0) here, 184 since we truncate the divisor at each step, but since {np,n} < D 185 originally, the largest possible partial quotient is B-1. */ 186 if (MPFR_UNLIKELY(np[n-1] > d1 || (np[n-1] == d1 && np[n-2] >= d0))) 187 q2 = MPFR_LIMB_MAX; 188 else 189 udiv_qr_3by2 (q2, q1, q0, np[n - 1], np[n - 2], np[n - 3], 190 d1, d0, dinv2.inv32); 191 /* since q2 = floor((np[n-1]*B^2+np[n-2]*B+np[n-3])/(d1*B+d0)), 192 we have q2 <= (np[n-1]*B^2+np[n-2]*B+np[n-3])/(d1*B+d0), 193 thus np[n-1]*B^2+np[n-2]*B+np[n-3] >= q2*(d1*B+d0) 194 and {np-1, n} >= q2*D - q2*B^(n-2) >= q2*D - B^(n-1) 195 thus {np-1, n} - (q2-1)*D >= D - B^(n-1) >= 0 196 which proves that at most one correction is needed */ 197 q0 = mpn_submul_1 (np - 1, dp, n, q2); 198 if (MPFR_UNLIKELY(q0 > np[n - 1])) 199 { 200 mpn_add_n (np - 1, np - 1, dp, n); 201 q2 --; 202 } 203 qp[--n] = q2; 204 dp ++; 205 } 206 207 /* we have B+dinv2 = floor((B^3-1)/(d1*B+d0)) < B^2/d1 208 q1 = floor(np[0]*(B+dinv2)/B) <= floor(np[0]*B/d1) 209 <= floor((np[0]*B+np[1])/d1) 210 thus q1 is not larger than the true quotient. 211 q1 > np[0]*(B+dinv2)/B - 1 > np[0]*(B^3-1)/(d1*B+d0)/B - 2 212 For d1*B+d0 <> B^2/2, we have B+dinv2 = floor(B^3/(d1*B+d0)) 213 thus q1 > np[0]*B^2/(d1*B+d0) - 2, i.e., 214 (d1*B+d0)*q1 > np[0]*B^2 - 2*(d1*B+d0) 215 d1*B*q1 > np[0]*B^2 - 2*d1*B - 2*d0 - d0*q1 >= np[0]*B^2 - 2*d1*B - B^2 216 thus q1 > np[0]*B/d1 - 2 - B/d1 > np[0]*B/d1 - 4. 217 218 For d1*B+d0 = B^2/2, dinv2 = B-1 thus q1 > np[0]*(2B-1)/B - 1 > 219 np[0]*B/d1 - 2. 220 221 In all cases, if q = floor((np[0]*B+np[1])/d1), we have: 222 q - 4 <= q1 <= q 223 */ 224 umul_ppmm (q1, q0, np[0], dinv2.inv32); 225 qp[0] = np[0] + q1; 226 227 return qh; 228 } 229 230 /* Put in {qp, n} an approximation of N={np, 2*n} divided by D={dp, n}, 231 with the most significant limb of the quotient as return value (0 or 1). 232 Assumes the most significant bit of D is set. Clobbers N. 233 234 This implements the ShortDiv algorithm from reference [1]. 235 236 Assumes n >= 2 (which should be fulfilled also in the recursive calls). 237 */ 238 mp_limb_t 239 mpfr_divhigh_n (mpfr_limb_ptr qp, mpfr_limb_ptr np, mpfr_limb_ptr dp, 240 mp_size_t n) 241 { 242 mp_size_t k, l; 243 mp_limb_t qh, cy; 244 mpfr_limb_ptr tp; 245 MPFR_TMP_DECL(marker); 246 247 MPFR_STAT_STATIC_ASSERT (MPFR_DIVHIGH_TAB_SIZE >= 15); /* so that 2*(n/3) >= (n+4)/2 */ 248 MPFR_ASSERTD(n >= 2); 249 k = MPFR_LIKELY (n < MPFR_DIVHIGH_TAB_SIZE) ? divhigh_ktab[n] : 2*(n/3); 250 251 if (k == 0) 252 { 253 #if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_SBPI1_DIVAPPR_Q) 254 mpfr_pi1_t dinv2; 255 invert_pi1 (dinv2, dp[n - 1], dp[n - 2]); 256 if (n > 2) /* sbpi1_divappr_q wants n > 2 */ 257 return __gmpn_sbpi1_divappr_q (qp, np, n + n, dp, n, dinv2.inv32); 258 else 259 return mpfr_divhigh_n_basecase (qp, np, dp, n); 260 #else /* use our own code for base-case short division */ 261 return mpfr_divhigh_n_basecase (qp, np, dp, n); 262 #endif 263 } 264 265 /* Check the bounds from [1]. In addition, we forbid k=n-1, which would 266 give l=1 in the recursive call. It follows n >= 5. */ 267 MPFR_ASSERTD ((n+4)/2 <= k && k < n-1); 268 269 MPFR_TMP_MARK (marker); 270 l = n - k; 271 /* first divide the most significant 2k limbs from N by the most significant 272 k limbs of D */ 273 qh = mpn_divrem (qp + l, 0, np + 2 * l, 2 * k, dp + l, k); /* exact */ 274 275 /* it remains {np,2l+k} = {np,n+l} as remainder */ 276 277 /* now we have to subtract high(Q1)*D0 where Q1=qh*B^k+{qp+l,k} and 278 D0={dp,l} */ 279 tp = MPFR_TMP_LIMBS_ALLOC (2 * l); 280 mpfr_mulhigh_n (tp, qp + k, dp, l); 281 /* we are only interested in the upper l limbs from {tp,2l} */ 282 cy = mpn_sub_n (np + n, np + n, tp + l, l); 283 if (qh) 284 cy += mpn_sub_n (np + n, np + n, dp, l); 285 while (cy > 0) /* Q1 was too large: subtract 1 to Q1 and add D to np+l */ 286 { 287 qh -= mpn_sub_1 (qp + l, qp + l, k, MPFR_LIMB_ONE); 288 cy -= mpn_add_n (np + l, np + l, dp, n); 289 } 290 291 /* now it remains {np,n+l} to divide by D */ 292 cy = mpfr_divhigh_n (qp, np + k, dp + k, l); 293 qh += mpn_add_1 (qp + l, qp + l, k, cy); 294 MPFR_TMP_FREE(marker); 295 296 return qh; 297 } 298
Indexes created Fri Jun 05 00:26:10 UTC 2026