1 // SPDX-License-Identifier: GPL-2.0-or-later 2 /* mpihelp-div.c - MPI helper functions 3 * Copyright (C) 1994, 1996 Free Software Foundation, Inc. 4 * Copyright (C) 1998, 1999 Free Software Foundation, Inc. 5 * 6 * This file is part of GnuPG. 7 * 8 * Note: This code is heavily based on the GNU MP Library. 9 * Actually it's the same code with only minor changes in the 10 * way the data is stored; this is to support the abstraction 11 * of an optional secure memory allocation which may be used 12 * to avoid revealing of sensitive data due to paging etc. 13 * The GNU MP Library itself is published under the LGPL; 14 * however I decided to publish this code under the plain GPL. 15 */ 16 17 #include "mpi-internal.h" 18 #include "longlong.h" 19 20 #ifndef UMUL_TIME 21 #define UMUL_TIME 1 22 #endif 23 #ifndef UDIV_TIME 24 #define UDIV_TIME UMUL_TIME 25 #endif 26 27 28 mpi_limb_t 29 mpihelp_mod_1(mpi_ptr_t dividend_ptr, mpi_size_t dividend_size, 30 mpi_limb_t divisor_limb) 31 { 32 mpi_size_t i; 33 mpi_limb_t n1, n0, r; 34 mpi_limb_t dummy __maybe_unused; 35 36 /* Botch: Should this be handled at all? Rely on callers? */ 37 if (!dividend_size) 38 return 0; 39 40 /* If multiplication is much faster than division, and the 41 * dividend is large, pre-invert the divisor, and use 42 * only multiplications in the inner loop. 43 * 44 * This test should be read: 45 * Does it ever help to use udiv_qrnnd_preinv? 46 * && Does what we save compensate for the inversion overhead? 47 */ 48 if (UDIV_TIME > (2 * UMUL_TIME + 6) 49 && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) { 50 int normalization_steps; 51 52 normalization_steps = count_leading_zeros(divisor_limb); 53 if (normalization_steps) { 54 mpi_limb_t divisor_limb_inverted; 55 56 divisor_limb <<= normalization_steps; 57 58 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The 59 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the 60 * most significant bit (with weight 2**N) implicit. 61 * 62 * Special case for DIVISOR_LIMB == 100...000. 63 */ 64 if (!(divisor_limb << 1)) 65 divisor_limb_inverted = ~(mpi_limb_t)0; 66 else 67 udiv_qrnnd(divisor_limb_inverted, dummy, 68 -divisor_limb, 0, divisor_limb); 69 70 n1 = dividend_ptr[dividend_size - 1]; 71 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps); 72 73 /* Possible optimization: 74 * if (r == 0 75 * && divisor_limb > ((n1 << normalization_steps) 76 * | (dividend_ptr[dividend_size - 2] >> ...))) 77 * ...one division less... 78 */ 79 for (i = dividend_size - 2; i >= 0; i--) { 80 n0 = dividend_ptr[i]; 81 UDIV_QRNND_PREINV(dummy, r, r, 82 ((n1 << normalization_steps) 83 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))), 84 divisor_limb, divisor_limb_inverted); 85 n1 = n0; 86 } 87 UDIV_QRNND_PREINV(dummy, r, r, 88 n1 << normalization_steps, 89 divisor_limb, divisor_limb_inverted); 90 return r >> normalization_steps; 91 } else { 92 mpi_limb_t divisor_limb_inverted; 93 94 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The 95 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the 96 * most significant bit (with weight 2**N) implicit. 97 * 98 * Special case for DIVISOR_LIMB == 100...000. 99 */ 100 if (!(divisor_limb << 1)) 101 divisor_limb_inverted = ~(mpi_limb_t)0; 102 else 103 udiv_qrnnd(divisor_limb_inverted, dummy, 104 -divisor_limb, 0, divisor_limb); 105 106 i = dividend_size - 1; 107 r = dividend_ptr[i]; 108 109 if (r >= divisor_limb) 110 r = 0; 111 else 112 i--; 113 114 for ( ; i >= 0; i--) { 115 n0 = dividend_ptr[i]; 116 UDIV_QRNND_PREINV(dummy, r, r, 117 n0, divisor_limb, divisor_limb_inverted); 118 } 119 return r; 120 } 121 } else { 122 if (UDIV_NEEDS_NORMALIZATION) { 123 int normalization_steps; 124 125 normalization_steps = count_leading_zeros(divisor_limb); 126 if (normalization_steps) { 127 divisor_limb <<= normalization_steps; 128 129 n1 = dividend_ptr[dividend_size - 1]; 130 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps); 131 132 /* Possible optimization: 133 * if (r == 0 134 * && divisor_limb > ((n1 << normalization_steps) 135 * | (dividend_ptr[dividend_size - 2] >> ...))) 136 * ...one division less... 137 */ 138 for (i = dividend_size - 2; i >= 0; i--) { 139 n0 = dividend_ptr[i]; 140 udiv_qrnnd(dummy, r, r, 141 ((n1 << normalization_steps) 142 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))), 143 divisor_limb); 144 n1 = n0; 145 } 146 udiv_qrnnd(dummy, r, r, 147 n1 << normalization_steps, 148 divisor_limb); 149 return r >> normalization_steps; 150 } 151 } 152 /* No normalization needed, either because udiv_qrnnd doesn't require 153 * it, or because DIVISOR_LIMB is already normalized. 154 */ 155 i = dividend_size - 1; 156 r = dividend_ptr[i]; 157 158 if (r >= divisor_limb) 159 r = 0; 160 else 161 i--; 162 163 for (; i >= 0; i--) { 164 n0 = dividend_ptr[i]; 165 udiv_qrnnd(dummy, r, r, n0, divisor_limb); 166 } 167 return r; 168 } 169 } 170 171 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write 172 * the NSIZE-DSIZE least significant quotient limbs at QP 173 * and the DSIZE long remainder at NP. If QEXTRA_LIMBS is 174 * non-zero, generate that many fraction bits and append them after the 175 * other quotient limbs. 176 * Return the most significant limb of the quotient, this is always 0 or 1. 177 * 178 * Preconditions: 179 * 0. NSIZE >= DSIZE. 180 * 1. The most significant bit of the divisor must be set. 181 * 2. QP must either not overlap with the input operands at all, or 182 * QP + DSIZE >= NP must hold true. (This means that it's 183 * possible to put the quotient in the high part of NUM, right after the 184 * remainder in NUM. 185 * 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. 186 */ 187 188 mpi_limb_t 189 mpihelp_divrem(mpi_ptr_t qp, mpi_size_t qextra_limbs, 190 mpi_ptr_t np, mpi_size_t nsize, mpi_ptr_t dp, mpi_size_t dsize) 191 { 192 mpi_limb_t most_significant_q_limb = 0; 193 194 switch (dsize) { 195 case 0: 196 /* We are asked to divide by zero, so go ahead and do it! (To make 197 the compiler not remove this statement, return the value.) */ 198 /* 199 * existing clients of this function have been modified 200 * not to call it with dsize == 0, so this should not happen 201 */ 202 return 1 / dsize; 203 204 case 1: 205 { 206 mpi_size_t i; 207 mpi_limb_t n1; 208 mpi_limb_t d; 209 210 d = dp[0]; 211 n1 = np[nsize - 1]; 212 213 if (n1 >= d) { 214 n1 -= d; 215 most_significant_q_limb = 1; 216 } 217 218 qp += qextra_limbs; 219 for (i = nsize - 2; i >= 0; i--) 220 udiv_qrnnd(qp[i], n1, n1, np[i], d); 221 qp -= qextra_limbs; 222 223 for (i = qextra_limbs - 1; i >= 0; i--) 224 udiv_qrnnd(qp[i], n1, n1, 0, d); 225 226 np[0] = n1; 227 } 228 break; 229 230 case 2: 231 { 232 mpi_size_t i; 233 mpi_limb_t n1, n0, n2; 234 mpi_limb_t d1, d0; 235 236 np += nsize - 2; 237 d1 = dp[1]; 238 d0 = dp[0]; 239 n1 = np[1]; 240 n0 = np[0]; 241 242 if (n1 >= d1 && (n1 > d1 || n0 >= d0)) { 243 sub_ddmmss(n1, n0, n1, n0, d1, d0); 244 most_significant_q_limb = 1; 245 } 246 247 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) { 248 mpi_limb_t q; 249 mpi_limb_t r; 250 251 if (i >= qextra_limbs) 252 np--; 253 else 254 np[0] = 0; 255 256 if (n1 == d1) { 257 /* Q should be either 111..111 or 111..110. Need special 258 * treatment of this rare case as normal division would 259 * give overflow. */ 260 q = ~(mpi_limb_t) 0; 261 262 r = n0 + d1; 263 if (r < d1) { /* Carry in the addition? */ 264 add_ssaaaa(n1, n0, r - d0, 265 np[0], 0, d0); 266 qp[i] = q; 267 continue; 268 } 269 n1 = d0 - (d0 != 0 ? 1 : 0); 270 n0 = -d0; 271 } else { 272 udiv_qrnnd(q, r, n1, n0, d1); 273 umul_ppmm(n1, n0, d0, q); 274 } 275 276 n2 = np[0]; 277 q_test: 278 if (n1 > r || (n1 == r && n0 > n2)) { 279 /* The estimated Q was too large. */ 280 q--; 281 sub_ddmmss(n1, n0, n1, n0, 0, d0); 282 r += d1; 283 if (r >= d1) /* If not carry, test Q again. */ 284 goto q_test; 285 } 286 287 qp[i] = q; 288 sub_ddmmss(n1, n0, r, n2, n1, n0); 289 } 290 np[1] = n1; 291 np[0] = n0; 292 } 293 break; 294 295 default: 296 { 297 mpi_size_t i; 298 mpi_limb_t dX, d1, n0; 299 300 np += nsize - dsize; 301 dX = dp[dsize - 1]; 302 d1 = dp[dsize - 2]; 303 n0 = np[dsize - 1]; 304 305 if (n0 >= dX) { 306 if (n0 > dX 307 || mpihelp_cmp(np, dp, dsize - 1) >= 0) { 308 mpihelp_sub_n(np, np, dp, dsize); 309 n0 = np[dsize - 1]; 310 most_significant_q_limb = 1; 311 } 312 } 313 314 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) { 315 mpi_limb_t q; 316 mpi_limb_t n1, n2; 317 mpi_limb_t cy_limb; 318 319 if (i >= qextra_limbs) { 320 np--; 321 n2 = np[dsize]; 322 } else { 323 n2 = np[dsize - 1]; 324 MPN_COPY_DECR(np + 1, np, dsize - 1); 325 np[0] = 0; 326 } 327 328 if (n0 == dX) { 329 /* This might over-estimate q, but it's probably not worth 330 * the extra code here to find out. */ 331 q = ~(mpi_limb_t) 0; 332 } else { 333 mpi_limb_t r; 334 335 udiv_qrnnd(q, r, n0, np[dsize - 1], dX); 336 umul_ppmm(n1, n0, d1, q); 337 338 while (n1 > r 339 || (n1 == r 340 && n0 > np[dsize - 2])) { 341 q--; 342 r += dX; 343 if (r < dX) /* I.e. "carry in previous addition?" */ 344 break; 345 n1 -= n0 < d1; 346 n0 -= d1; 347 } 348 } 349 350 /* Possible optimization: We already have (q * n0) and (1 * n1) 351 * after the calculation of q. Taking advantage of that, we 352 * could make this loop make two iterations less. */ 353 cy_limb = mpihelp_submul_1(np, dp, dsize, q); 354 355 if (n2 != cy_limb) { 356 mpihelp_add_n(np, np, dp, dsize); 357 q--; 358 } 359 360 qp[i] = q; 361 n0 = np[dsize - 1]; 362 } 363 } 364 } 365 366 return most_significant_q_limb; 367 } 368 369 /**************** 370 * Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB. 371 * Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR. 372 * Return the single-limb remainder. 373 * There are no constraints on the value of the divisor. 374 * 375 * QUOT_PTR and DIVIDEND_PTR might point to the same limb. 376 */ 377 378 mpi_limb_t 379 mpihelp_divmod_1(mpi_ptr_t quot_ptr, 380 mpi_ptr_t dividend_ptr, mpi_size_t dividend_size, 381 mpi_limb_t divisor_limb) 382 { 383 mpi_size_t i; 384 mpi_limb_t n1, n0, r; 385 mpi_limb_t dummy __maybe_unused; 386 387 if (!dividend_size) 388 return 0; 389 390 /* If multiplication is much faster than division, and the 391 * dividend is large, pre-invert the divisor, and use 392 * only multiplications in the inner loop. 393 * 394 * This test should be read: 395 * Does it ever help to use udiv_qrnnd_preinv? 396 * && Does what we save compensate for the inversion overhead? 397 */ 398 if (UDIV_TIME > (2 * UMUL_TIME + 6) 399 && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) { 400 int normalization_steps; 401 402 normalization_steps = count_leading_zeros(divisor_limb); 403 if (normalization_steps) { 404 mpi_limb_t divisor_limb_inverted; 405 406 divisor_limb <<= normalization_steps; 407 408 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The 409 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the 410 * most significant bit (with weight 2**N) implicit. 411 */ 412 /* Special case for DIVISOR_LIMB == 100...000. */ 413 if (!(divisor_limb << 1)) 414 divisor_limb_inverted = ~(mpi_limb_t)0; 415 else 416 udiv_qrnnd(divisor_limb_inverted, dummy, 417 -divisor_limb, 0, divisor_limb); 418 419 n1 = dividend_ptr[dividend_size - 1]; 420 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps); 421 422 /* Possible optimization: 423 * if (r == 0 424 * && divisor_limb > ((n1 << normalization_steps) 425 * | (dividend_ptr[dividend_size - 2] >> ...))) 426 * ...one division less... 427 */ 428 for (i = dividend_size - 2; i >= 0; i--) { 429 n0 = dividend_ptr[i]; 430 UDIV_QRNND_PREINV(quot_ptr[i + 1], r, r, 431 ((n1 << normalization_steps) 432 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))), 433 divisor_limb, divisor_limb_inverted); 434 n1 = n0; 435 } 436 UDIV_QRNND_PREINV(quot_ptr[0], r, r, 437 n1 << normalization_steps, 438 divisor_limb, divisor_limb_inverted); 439 return r >> normalization_steps; 440 } else { 441 mpi_limb_t divisor_limb_inverted; 442 443 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The 444 * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the 445 * most significant bit (with weight 2**N) implicit. 446 */ 447 /* Special case for DIVISOR_LIMB == 100...000. */ 448 if (!(divisor_limb << 1)) 449 divisor_limb_inverted = ~(mpi_limb_t) 0; 450 else 451 udiv_qrnnd(divisor_limb_inverted, dummy, 452 -divisor_limb, 0, divisor_limb); 453 454 i = dividend_size - 1; 455 r = dividend_ptr[i]; 456 457 if (r >= divisor_limb) 458 r = 0; 459 else 460 quot_ptr[i--] = 0; 461 462 for ( ; i >= 0; i--) { 463 n0 = dividend_ptr[i]; 464 UDIV_QRNND_PREINV(quot_ptr[i], r, r, 465 n0, divisor_limb, divisor_limb_inverted); 466 } 467 return r; 468 } 469 } else { 470 if (UDIV_NEEDS_NORMALIZATION) { 471 int normalization_steps; 472 473 normalization_steps = count_leading_zeros(divisor_limb); 474 if (normalization_steps) { 475 divisor_limb <<= normalization_steps; 476 477 n1 = dividend_ptr[dividend_size - 1]; 478 r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps); 479 480 /* Possible optimization: 481 * if (r == 0 482 * && divisor_limb > ((n1 << normalization_steps) 483 * | (dividend_ptr[dividend_size - 2] >> ...))) 484 * ...one division less... 485 */ 486 for (i = dividend_size - 2; i >= 0; i--) { 487 n0 = dividend_ptr[i]; 488 udiv_qrnnd(quot_ptr[i + 1], r, r, 489 ((n1 << normalization_steps) 490 | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))), 491 divisor_limb); 492 n1 = n0; 493 } 494 udiv_qrnnd(quot_ptr[0], r, r, 495 n1 << normalization_steps, 496 divisor_limb); 497 return r >> normalization_steps; 498 } 499 } 500 /* No normalization needed, either because udiv_qrnnd doesn't require 501 * it, or because DIVISOR_LIMB is already normalized. 502 */ 503 i = dividend_size - 1; 504 r = dividend_ptr[i]; 505 506 if (r >= divisor_limb) 507 r = 0; 508 else 509 quot_ptr[i--] = 0; 510 511 for (; i >= 0; i--) { 512 n0 = dividend_ptr[i]; 513 udiv_qrnnd(quot_ptr[i], r, r, n0, divisor_limb); 514 } 515 return r; 516 } 517 } 518