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                 return 1 / dsize;
 221
 222         case 1:
 223                 {
 224                         mpi_size_t i;
 225                         mpi_limb_t n1;
 226                         mpi_limb_t d;
 227
 228                         d = dp[0];
 229                         n1 = np[nsize - 1];
 230
 231                         if (n1 >= d) {
 232                                 n1 -= d;
 233                                 most_significant_q_limb = 1;
 234                         }
 235
 236                         qp += qextra_limbs;
 237                         for (i = nsize - 2; i >= 0; i--)
 238                                 udiv_qrnnd(qp[i], n1, n1, np[i], d);
 239                         qp -= qextra_limbs;
 240
 241                         for (i = qextra_limbs - 1; i >= 0; i--)
 242                                 udiv_qrnnd(qp[i], n1, n1, 0, d);
 243
 244                         np[0] = n1;
 245                 }
 246                 break;
 247
 248         case 2:
 249                 {
 250                         mpi_size_t i;
 251                         mpi_limb_t n1, n0, n2;
 252                         mpi_limb_t d1, d0;
 253
 254                         np += nsize - 2;
 255                         d1 = dp[1];
 256                         d0 = dp[0];
 257                         n1 = np[1];
 258                         n0 = np[0];
 259
 260                         if (n1 >= d1 && (n1 > d1 || n0 >= d0)) {
 261                                 sub_ddmmss(n1, n0, n1, n0, d1, d0);
 262                                 most_significant_q_limb = 1;
 263                         }
 264
 265                         for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) {
 266                                 mpi_limb_t q;
 267                                 mpi_limb_t r;
 268
 269                                 if (i >= qextra_limbs)
 270                                         np--;
 271                                 else
 272                                         np[0] = 0;
 273
 274                                 if (n1 == d1) {
 275                                         /* Q should be either 111..111 or 111..110.  Need special
 276                                          * treatment of this rare case as normal division would
 277                                          * give overflow.  */
 278                                         q = ~(mpi_limb_t) 0;
 279
 280                                         r = n0 + d1;
 281                                         if (r < d1) {   /* Carry in the addition? */
 282                                                 add_ssaaaa(n1, n0, r - d0,
 283                                                            np[0], 0, d0);
 284                                                 qp[i] = q;
 285                                                 continue;
 286                                         }
 287                                         n1 = d0 - (d0 != 0 ? 1 : 0);
 288                                         n0 = -d0;
 289                                 } else {
 290                                         udiv_qrnnd(q, r, n1, n0, d1);
 291                                         umul_ppmm(n1, n0, d0, q);
 292                                 }
 293
 294                                 n2 = np[0];
 295 q_test:
 296                                 if (n1 > r || (n1 == r && n0 > n2)) {
 297                                         /* The estimated Q was too large.  */
 298                                         q--;
 299                                         sub_ddmmss(n1, n0, n1, n0, 0, d0);
 300                                         r += d1;
 301                                         if (r >= d1)    /* If not carry, test Q again.  */
 302                                                 goto q_test;
 303                                 }
 304
 305                                 qp[i] = q;
 306                                 sub_ddmmss(n1, n0, r, n2, n1, n0);
 307                         }
 308                         np[1] = n1;
 309                         np[0] = n0;
 310                 }
 311                 break;
 312
 313         default:
 314                 {
 315                         mpi_size_t i;
 316                         mpi_limb_t dX, d1, n0;
 317
 318                         np += nsize - dsize;
 319                         dX = dp[dsize - 1];
 320                         d1 = dp[dsize - 2];
 321                         n0 = np[dsize - 1];
 322
 323                         if (n0 >= dX) {
 324                                 if (n0 > dX
 325                                     || mpihelp_cmp(np, dp, dsize - 1) >= 0) {
 326                                         mpihelp_sub_n(np, np, dp, dsize);
 327                                         n0 = np[dsize - 1];
 328                                         most_significant_q_limb = 1;
 329                                 }
 330                         }
 331
 332                         for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) {
 333                                 mpi_limb_t q;
 334                                 mpi_limb_t n1, n2;
 335                                 mpi_limb_t cy_limb;
 336
 337                                 if (i >= qextra_limbs) {
 338                                         np--;
 339                                         n2 = np[dsize];
 340                                 } else {
 341                                         n2 = np[dsize - 1];
 342                                         MPN_COPY_DECR(np + 1, np, dsize - 1);
 343                                         np[0] = 0;
 344                                 }
 345
 346                                 if (n0 == dX) {
 347                                         /* This might over-estimate q, but it's probably not worth
 348                                          * the extra code here to find out.  */
 349                                         q = ~(mpi_limb_t) 0;
 350                                 } else {
 351                                         mpi_limb_t r;
 352
 353                                         udiv_qrnnd(q, r, n0, np[dsize - 1], dX);
 354                                         umul_ppmm(n1, n0, d1, q);
 355
 356                                         while (n1 > r
 357                                                || (n1 == r
 358                                                    && n0 > np[dsize - 2])) {
 359                                                 q--;
 360                                                 r += dX;
 361                                                 if (r < dX)     /* I.e. "carry in previous addition?" */
 362                                                         break;
 363                                                 n1 -= n0 < d1;
 364                                                 n0 -= d1;
 365                                         }
 366                                 }
 367
 368                                 /* Possible optimization: We already have (q * n0) and (1 * n1)
 369                                  * after the calculation of q.  Taking advantage of that, we
 370                                  * could make this loop make two iterations less.  */
 371                                 cy_limb = mpihelp_submul_1(np, dp, dsize, q);
 372
 373                                 if (n2 != cy_limb) {
 374                                         mpihelp_add_n(np, np, dp, dsize);
 375                                         q--;
 376                                 }
 377
 378                                 qp[i] = q;
 379                                 n0 = np[dsize - 1];
 380                         }
 381                 }
 382         }
 383
 384         return most_significant_q_limb;
 385 }
 386
 387 /****************
 388  * Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
 389  * Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
 390  * Return the single-limb remainder.
 391  * There are no constraints on the value of the divisor.
 392  *
 393  * QUOT_PTR and DIVIDEND_PTR might point to the same limb.
 394  */
 395
 396 mpi_limb_t
 397 mpihelp_divmod_1(mpi_ptr_t quot_ptr,
 398                  mpi_ptr_t dividend_ptr, mpi_size_t dividend_size,
 399                  mpi_limb_t divisor_limb)
 400 {
 401         mpi_size_t i;
 402         mpi_limb_t n1, n0, r;
 403         int dummy;
 404
 405         if (!dividend_size)
 406                 return 0;
 407
 408         /* If multiplication is much faster than division, and the
 409          * dividend is large, pre-invert the divisor, and use
 410          * only multiplications in the inner loop.
 411          *
 412          * This test should be read:
 413          * Does it ever help to use udiv_qrnnd_preinv?
 414          * && Does what we save compensate for the inversion overhead?
 415          */
 416         if (UDIV_TIME > (2 * UMUL_TIME + 6)
 417             && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) {
 418                 int normalization_steps;
 419
 420                 count_leading_zeros(normalization_steps, divisor_limb);
 421                 if (normalization_steps) {
 422                         mpi_limb_t divisor_limb_inverted;
 423
 424                         divisor_limb <<= normalization_steps;
 425
 426                         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
 427                          * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
 428                          * most significant bit (with weight 2**N) implicit.
 429                          */
 430                         /* Special case for DIVISOR_LIMB == 100...000.  */
 431                         if (!(divisor_limb << 1))
 432                                 divisor_limb_inverted = ~(mpi_limb_t) 0;
 433                         else
 434                                 udiv_qrnnd(divisor_limb_inverted, dummy,
 435                                            -divisor_limb, 0, divisor_limb);
 436
 437                         n1 = dividend_ptr[dividend_size - 1];
 438                         r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
 439
 440                         /* Possible optimization:
 441                          * if (r == 0
 442                          * && divisor_limb > ((n1 << normalization_steps)
 443                          *                 | (dividend_ptr[dividend_size - 2] >> ...)))
 444                          * ...one division less...
 445                          */
 446                         for (i = dividend_size - 2; i >= 0; i--) {
 447                                 n0 = dividend_ptr[i];
 448                                 UDIV_QRNND_PREINV(quot_ptr[i + 1], r, r,
 449                                                   ((n1 << normalization_steps)
 450                                                    | (n0 >>
 451                                                       (BITS_PER_MPI_LIMB -
 452                                                        normalization_steps))),
 453                                                   divisor_limb,
 454                                                   divisor_limb_inverted);
 455                                 n1 = n0;
 456                         }
 457                         UDIV_QRNND_PREINV(quot_ptr[0], r, r,
 458                                           n1 << normalization_steps,
 459                                           divisor_limb, divisor_limb_inverted);
 460                         return r >> normalization_steps;
 461                 } else {
 462                         mpi_limb_t divisor_limb_inverted;
 463
 464                         /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
 465                          * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
 466                          * most significant bit (with weight 2**N) implicit.
 467                          */
 468                         /* Special case for DIVISOR_LIMB == 100...000.  */
 469                         if (!(divisor_limb << 1))
 470                                 divisor_limb_inverted = ~(mpi_limb_t) 0;
 471                         else
 472                                 udiv_qrnnd(divisor_limb_inverted, dummy,
 473                                            -divisor_limb, 0, divisor_limb);
 474
 475                         i = dividend_size - 1;
 476                         r = dividend_ptr[i];
 477
 478                         if (r >= divisor_limb)
 479                                 r = 0;
 480                         else
 481                                 quot_ptr[i--] = 0;
 482
 483                         for (; i >= 0; i--) {
 484                                 n0 = dividend_ptr[i];
 485                                 UDIV_QRNND_PREINV(quot_ptr[i], r, r,
 486                                                   n0, divisor_limb,
 487                                                   divisor_limb_inverted);
 488                         }
 489                         return r;
 490                 }
 491         } else {
 492                 if (UDIV_NEEDS_NORMALIZATION) {
 493                         int normalization_steps;
 494
 495                         count_leading_zeros(normalization_steps, divisor_limb);
 496                         if (normalization_steps) {
 497                                 divisor_limb <<= normalization_steps;
 498
 499                                 n1 = dividend_ptr[dividend_size - 1];
 500                                 r = n1 >> (BITS_PER_MPI_LIMB -
 501                                            normalization_steps);
 502
 503                                 /* Possible optimization:
 504                                  * if (r == 0
 505                                  * && divisor_limb > ((n1 << normalization_steps)
 506                                  *                 | (dividend_ptr[dividend_size - 2] >> ...)))
 507                                  * ...one division less...
 508                                  */
 509                                 for (i = dividend_size - 2; i >= 0; i--) {
 510                                         n0 = dividend_ptr[i];
 511                                         udiv_qrnnd(quot_ptr[i + 1], r, r,
 512                                                    ((n1 << normalization_steps)
 513                                                     | (n0 >>
 514                                                        (BITS_PER_MPI_LIMB -
 515                                                         normalization_steps))),
 516                                                    divisor_limb);
 517                                         n1 = n0;
 518                                 }
 519                                 udiv_qrnnd(quot_ptr[0], r, r,
 520                                            n1 << normalization_steps,
 521                                            divisor_limb);
 522                                 return r >> normalization_steps;
 523                         }
 524                 }
 525                 /* No normalization needed, either because udiv_qrnnd doesn't require
 526                  * it, or because DIVISOR_LIMB is already normalized.  */
 527                 i = dividend_size - 1;
 528                 r = dividend_ptr[i];
 529
 530                 if (r >= divisor_limb)
 531                         r = 0;
 532                 else
 533                         quot_ptr[i--] = 0;
 534
 535                 for (; i >= 0; i--) {
 536                         n0 = dividend_ptr[i];
 537                         udiv_qrnnd(quot_ptr[i], r, r, n0, divisor_limb);
 538                 }
 539                 return r;
 540         }
 541 }