lib/mpi/mpih-div.c

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