chiark - git - mdw - catacomb/blob - mpx.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mpx.c,v 1.7 1999/12/22 15:49:07 mdw Exp $
   4  *
   5  * Low-level multiprecision arithmetic
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mpx.c,v $
  33  * Revision 1.7  1999/12/22 15:49:07  mdw
  34  * New function for division by a small integer.
  35  *
  36  * Revision 1.6  1999/11/20 22:43:44  mdw
  37  * Integrate testing for MPX routines.
  38  *
  39  * Revision 1.5  1999/11/20 22:23:27  mdw
  40  * Add function versions of some low-level macros with wider use.
  41  *
  42  * Revision 1.4  1999/11/17 18:04:09  mdw
  43  * Add two's-complement functionality.  Improve mpx_udiv a little by
  44  * performing the multiplication of the divisor by q with the subtraction
  45  * from r.
  46  *
  47  * Revision 1.3  1999/11/13 01:57:31  mdw
  48  * Remove stray debugging code.
  49  *
  50  * Revision 1.2  1999/11/13 01:50:59  mdw
  51  * Multiprecision routines finished and tested.
  52  *
  53  * Revision 1.1  1999/09/03 08:41:12  mdw
  54  * Initial import.
  55  *
  56  */
  57
  58 /*----- Header files ------------------------------------------------------*/
  59
  60 #include <assert.h>
  61 #include <stdio.h>
  62 #include <stdlib.h>
  63 #include <string.h>
  64
  65 #include <mLib/bits.h>
  66
  67 #include "mptypes.h"
  68 #include "mpx.h"
  69
  70 /*----- Loading and storing -----------------------------------------------*/
  71
  72 /* --- @mpx_storel@ --- *
  73  *
  74  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
  75  *              @void *pp@ = pointer to octet array
  76  *              @size_t sz@ = size of octet array
  77  *
  78  * Returns:     ---
  79  *
  80  * Use:         Stores an MP in an octet array, least significant octet
  81  *              first.  High-end octets are silently discarded if there
  82  *              isn't enough space for them.
  83  */
  84
  85 void mpx_storel(const mpw *v, const mpw *vl, void *pp, size_t sz)
  86 {
  87   mpw n, w = 0;
  88   octet *p = pp, *q = p + sz;
  89   unsigned bits = 0;
  90
  91   while (p < q) {
  92     if (bits < 8) {
  93       if (v >= vl) {
  94         *p++ = U8(w);
  95         break;
  96       }
  97       n = *v++;
  98       *p++ = U8(w | n << bits);
  99       w = n >> (8 - bits);
 100       bits += MPW_BITS - 8;
 101     } else {
 102       *p++ = U8(w);
 103       w >>= 8;
 104       bits -= 8;
 105     }
 106   }
 107   memset(p, 0, q - p);
 108 }
 109
 110 /* --- @mpx_loadl@ --- *
 111  *
 112  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 113  *              @const void *pp@ = pointer to octet array
 114  *              @size_t sz@ = size of octet array
 115  *
 116  * Returns:     ---
 117  *
 118  * Use:         Loads an MP in an octet array, least significant octet
 119  *              first.  High-end octets are ignored if there isn't enough
 120  *              space for them.
 121  */
 122
 123 void mpx_loadl(mpw *v, mpw *vl, const void *pp, size_t sz)
 124 {
 125   unsigned n;
 126   mpw w = 0;
 127   const octet *p = pp, *q = p + sz;
 128   unsigned bits = 0;
 129
 130   if (v >= vl)
 131     return;
 132   while (p < q) {
 133     n = U8(*p++);
 134     w |= n << bits;
 135     bits += 8;
 136     if (bits >= MPW_BITS) {
 137       *v++ = MPW(w);
 138       w = n >> (MPW_BITS - bits + 8);
 139       bits -= MPW_BITS;
 140       if (v >= vl)
 141         return;
 142     }
 143   }
 144   *v++ = w;
 145   MPX_ZERO(v, vl);
 146 }
 147
 148 /* --- @mpx_storeb@ --- *
 149  *
 150  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 151  *              @void *pp@ = pointer to octet array
 152  *              @size_t sz@ = size of octet array
 153  *
 154  * Returns:     ---
 155  *
 156  * Use:         Stores an MP in an octet array, most significant octet
 157  *              first.  High-end octets are silently discarded if there
 158  *              isn't enough space for them.
 159  */
 160
 161 void mpx_storeb(const mpw *v, const mpw *vl, void *pp, size_t sz)
 162 {
 163   mpw n, w = 0;
 164   octet *p = pp, *q = p + sz;
 165   unsigned bits = 0;
 166
 167   while (q > p) {
 168     if (bits < 8) {
 169       if (v >= vl) {
 170         *--q = U8(w);
 171         break;
 172       }
 173       n = *v++;
 174       *--q = U8(w | n << bits);
 175       w = n >> (8 - bits);
 176       bits += MPW_BITS - 8;
 177     } else {
 178       *--q = U8(w);
 179       w >>= 8;
 180       bits -= 8;
 181     }
 182   }
 183   memset(p, 0, q - p);
 184 }
 185
 186 /* --- @mpx_loadb@ --- *
 187  *
 188  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 189  *              @const void *pp@ = pointer to octet array
 190  *              @size_t sz@ = size of octet array
 191  *
 192  * Returns:     ---
 193  *
 194  * Use:         Loads an MP in an octet array, most significant octet
 195  *              first.  High-end octets are ignored if there isn't enough
 196  *              space for them.
 197  */
 198
 199 void mpx_loadb(mpw *v, mpw *vl, const void *pp, size_t sz)
 200 {
 201   unsigned n;
 202   mpw w = 0;
 203   const octet *p = pp, *q = p + sz;
 204   unsigned bits = 0;
 205
 206   if (v >= vl)
 207     return;
 208   while (q > p) {
 209     n = U8(*--q);
 210     w |= n << bits;
 211     bits += 8;
 212     if (bits >= MPW_BITS) {
 213       *v++ = MPW(w);
 214       w = n >> (MPW_BITS - bits + 8);
 215       bits -= MPW_BITS;
 216       if (v >= vl)
 217         return;
 218     }
 219   }
 220   *v++ = w;
 221   MPX_ZERO(v, vl);
 222 }
 223
 224 /*----- Logical shifting --------------------------------------------------*/
 225
 226 /* --- @mpx_lsl@ --- *
 227  *
 228  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 229  *              @const mpw *av, *avl@ = source vector base and limit
 230  *              @size_t n@ = number of bit positions to shift by
 231  *
 232  * Returns:     ---
 233  *
 234  * Use:         Performs a logical shift left operation on an integer.
 235  */
 236
 237 void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 238 {
 239   size_t nw;
 240   unsigned nb;
 241
 242   /* --- Trivial special case --- */
 243
 244   if (n == 0)
 245     MPX_COPY(dv, dvl, av, avl);
 246
 247   /* --- Single bit shifting --- */
 248
 249   else if (n == 1) {
 250     mpw w = 0;
 251     while (av < avl) {
 252       mpw t;
 253       if (dv >= dvl)
 254         goto done;
 255       t = *av++;
 256       *dv++ = MPW((t << 1) | w);
 257       w = t >> (MPW_BITS - 1);
 258     }
 259     if (dv >= dvl)
 260       goto done;
 261     *dv++ = MPW(w);
 262     MPX_ZERO(dv, dvl);
 263     goto done;
 264   }
 265
 266   /* --- Break out word and bit shifts for more sophisticated work --- */
 267
 268   nw = n / MPW_BITS;
 269   nb = n % MPW_BITS;
 270
 271   /* --- Handle a shift by a multiple of the word size --- */
 272
 273   if (nb == 0) {
 274     MPX_COPY(dv + nw, dvl, av, avl);
 275     memset(dv, 0, MPWS(nw));
 276   }
 277
 278   /* --- And finally the difficult case --- *
 279    *
 280    * This is a little convoluted, because I have to start from the end and
 281    * work backwards to avoid overwriting the source, if they're both the same
 282    * block of memory.
 283    */
 284
 285   else {
 286     mpw w;
 287     size_t nr = MPW_BITS - nb;
 288     size_t dvn = dvl - dv;
 289     size_t avn = avl - av;
 290
 291     if (dvn <= nw) {
 292       MPX_ZERO(dv, dvl);
 293       goto done;
 294     }
 295
 296     if (dvn > avn + nw) {
 297       size_t off = avn + nw + 1;
 298       MPX_ZERO(dv + off, dvl);
 299       dvl = dv + off;
 300       w = 0;
 301     } else {
 302       avl = av + dvn - nw;
 303       w = *--avl << nb;
 304     }
 305
 306     while (avl > av) {
 307       mpw t = *--avl;
 308       *--dvl = (t >> nr) | w;
 309       w = t << nb;
 310     }
 311
 312     *--dvl = w;
 313     MPX_ZERO(dv, dvl);
 314   }
 315
 316 done:;
 317 }
 318
 319 /* --- @mpx_lsr@ --- *
 320  *
 321  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 322  *              @const mpw *av, *avl@ = source vector base and limit
 323  *              @size_t n@ = number of bit positions to shift by
 324  *
 325  * Returns:     ---
 326  *
 327  * Use:         Performs a logical shift right operation on an integer.
 328  */
 329
 330 void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 331 {
 332   size_t nw;
 333   unsigned nb;
 334
 335   /* --- Trivial special case --- */
 336
 337   if (n == 0)
 338     MPX_COPY(dv, dvl, av, avl);
 339
 340   /* --- Single bit shifting --- */
 341
 342   else if (n == 1) {
 343     mpw w = *av++ >> 1;
 344     while (av < avl) {
 345       mpw t;
 346       if (dv >= dvl)
 347         goto done;
 348       t = *av++;
 349       *dv++ = MPW((t << (MPW_BITS - 1)) | w);
 350       w = t >> 1;
 351     }
 352     if (dv >= dvl)
 353       goto done;
 354     *dv++ = MPW(w);
 355     MPX_ZERO(dv, dvl);
 356     goto done;
 357   }
 358
 359   /* --- Break out word and bit shifts for more sophisticated work --- */
 360
 361   nw = n / MPW_BITS;
 362   nb = n % MPW_BITS;
 363
 364   /* --- Handle a shift by a multiple of the word size --- */
 365
 366   if (nb == 0)
 367     MPX_COPY(dv, dvl, av + nw, avl);
 368
 369   /* --- And finally the difficult case --- */
 370
 371   else {
 372     mpw w;
 373     size_t nr = MPW_BITS - nb;
 374
 375     av += nw;
 376     w = *av++;
 377     while (av < avl) {
 378       mpw t;
 379       if (dv >= dvl)
 380         goto done;
 381       t = *av++;
 382       *dv++ = MPW((w >> nb) | (t << nr));
 383       w = t;
 384     }
 385     if (dv < dvl) {
 386       *dv++ = MPW(w >> nb);
 387       MPX_ZERO(dv, dvl);
 388     }
 389   }
 390
 391 done:;
 392 }
 393
 394 /*----- Unsigned arithmetic -----------------------------------------------*/
 395
 396 /* --- @mpx_2c@ --- *
 397  *
 398  * Arguments:   @mpw *dv, *dvl@ = destination vector
 399  *              @const mpw *v, *vl@ = source vector
 400  *
 401  * Returns:     ---
 402  *
 403  * Use:         Calculates the two's complement of @v@.
 404  */
 405
 406 void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
 407 {
 408   mpw c = 0;
 409   while (dv < dvl && v < vl)
 410     *dv++ = c = MPW(~*v++);
 411   if (dv < dvl) {
 412     if (c > MPW_MAX / 2)
 413       c = MPW(~0);
 414     while (dv < dvl)
 415       *dv++ = c;
 416   }
 417   MPX_UADDN(dv, dvl, 1);
 418 }
 419
 420 /* --- @mpx_ucmp@ --- *
 421  *
 422  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 423  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 424  *
 425  * Returns:     Less than, equal to, or greater than zero depending on
 426  *              whether @a@ is less than, equal to or greater than @b@,
 427  *              respectively.
 428  *
 429  * Use:         Performs an unsigned integer comparison.
 430  */
 431
 432 int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
 433 {
 434   MPX_SHRINK(av, avl);
 435   MPX_SHRINK(bv, bvl);
 436
 437   if (avl - av > bvl - bv)
 438     return (+1);
 439   else if (avl - av < bvl - bv)
 440     return (-1);
 441   else while (avl > av) {
 442     mpw a = *--avl, b = *--bvl;
 443     if (a > b)
 444       return (+1);
 445     else if (a < b)
 446       return (-1);
 447   }
 448   return (0);
 449 }
 450
 451 /* --- @mpx_uadd@ --- *
 452  *
 453  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 454  *              @const mpw *av, *avl@ = first addend vector base and limit
 455  *              @const mpw *bv, *bvl@ = second addend vector base and limit
 456  *
 457  * Returns:     ---
 458  *
 459  * Use:         Performs unsigned integer addition.  If the result overflows
 460  *              the destination vector, high-order bits are discarded.  This
 461  *              means that two's complement addition happens more or less for
 462  *              free, although that's more a side-effect than anything else.
 463  *              The result vector may be equal to either or both source
 464  *              vectors, but may not otherwise overlap them.
 465  */
 466
 467 void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 468               const mpw *bv, const mpw *bvl)
 469 {
 470   mpw c = 0;
 471
 472   while (av < avl || bv < bvl) {
 473     mpw a, b;
 474     mpd x;
 475     if (dv >= dvl)
 476       return;
 477     a = (av < avl) ? *av++ : 0;
 478     b = (bv < bvl) ? *bv++ : 0;
 479     x = (mpd)a + (mpd)b + c;
 480     *dv++ = MPW(x);
 481     c = x >> MPW_BITS;
 482   }
 483   if (dv < dvl) {
 484     *dv++ = c;
 485     MPX_ZERO(dv, dvl);
 486   }
 487 }
 488
 489 /* --- @mpx_uaddn@ --- *
 490  *
 491  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 492  *              @mpw n@ = other addend
 493  *
 494  * Returns:     ---
 495  *
 496  * Use:         Adds a small integer to a multiprecision number.
 497  */
 498
 499 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
 500
 501 /* --- @mpx_usub@ --- *
 502  *
 503  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 504  *              @const mpw *av, *avl@ = first argument vector base and limit
 505  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 506  *
 507  * Returns:     ---
 508  *
 509  * Use:         Performs unsigned integer subtraction.  If the result
 510  *              overflows the destination vector, high-order bits are
 511  *              discarded.  This means that two's complement subtraction
 512  *              happens more or less for free, althuogh that's more a side-
 513  *              effect than anything else.  The result vector may be equal to
 514  *              either or both source vectors, but may not otherwise overlap
 515  *              them.
 516  */
 517
 518 void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 519               const mpw *bv, const mpw *bvl)
 520 {
 521   mpw c = 0;
 522
 523   while (av < avl || bv < bvl) {
 524     mpw a, b;
 525     mpd x;
 526     if (dv >= dvl)
 527       return;
 528     a = (av < avl) ? *av++ : 0;
 529     b = (bv < bvl) ? *bv++ : 0;
 530     x = (mpd)a - (mpd)b - c;
 531     *dv++ = MPW(x);
 532     if (x >> MPW_BITS)
 533       c = 1;
 534     else
 535       c = 0;
 536   }
 537   if (c)
 538     c = MPW_MAX;
 539   while (dv < dvl)
 540     *dv++ = c;
 541 }
 542
 543 /* --- @mpx_usubn@ --- *
 544  *
 545  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 546  *              @n@ = subtrahend
 547  *
 548  * Returns:     ---
 549  *
 550  * Use:         Subtracts a small integer from a multiprecision number.
 551  */
 552
 553 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
 554
 555 /* --- @mpx_umul@ --- *
 556  *
 557  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 558  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 559  *              @const mpw *bv, *bvl@ = multiplier vector base and limit
 560  *
 561  * Returns:     ---
 562  *
 563  * Use:         Performs unsigned integer multiplication.  If the result
 564  *              overflows the desination vector, high-order bits are
 565  *              discarded.  The result vector may not overlap the argument
 566  *              vectors in any way.
 567  */
 568
 569 void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 570               const mpw *bv, const mpw *bvl)
 571 {
 572   /* --- This is probably worthwhile on a multiply --- */
 573
 574   MPX_SHRINK(av, avl);
 575   MPX_SHRINK(bv, bvl);
 576
 577   /* --- Deal with a multiply by zero --- */
 578
 579   if (bv == bvl) {
 580     MPX_ZERO(dv, dvl);
 581     return;
 582   }
 583
 584   /* --- Do the initial multiply and initialize the accumulator --- */
 585
 586   MPX_UMULN(dv, dvl, av, avl, *bv++);
 587
 588   /* --- Do the remaining multiply/accumulates --- */
 589
 590   while (dv < dvl && bv < bvl) {
 591     mpw m = *bv++;
 592     mpw c = 0;
 593     const mpw *avv = av;
 594     mpw *dvv = ++dv;
 595
 596     while (avv < avl) {
 597       mpd x;
 598       if (dvv >= dvl)
 599         goto next;
 600       x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
 601       *dvv++ = MPW(x);
 602       c = x >> MPW_BITS;
 603     }
 604     MPX_UADDN(dvv, dvl, c);
 605   next:;
 606   }
 607 }
 608
 609 /* --- @mpx_umuln@ --- *
 610  *
 611  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 612  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 613  *              @mpw m@ = multiplier
 614  *
 615  * Returns:     ---
 616  *
 617  * Use:         Multiplies a multiprecision integer by a single-word value.
 618  *              The destination and source may be equal.  The destination
 619  *              is completely cleared after use.
 620  */
 621
 622 void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
 623 {
 624   MPX_UMULN(dv, dvl, av, avl, m);
 625 }
 626
 627 /* --- @mpx_umlan@ --- *
 628  *
 629  * Arguments:   @mpw *dv, *dvl@ = destination/accumulator base and limit
 630  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 631  *              @mpw m@ = multiplier
 632  *
 633  * Returns:     ---
 634  *
 635  * Use:         Multiplies a multiprecision integer by a single-word value
 636  *              and adds the result to an accumulator.
 637  */
 638
 639 void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
 640 {
 641   MPX_UMLAN(dv, dvl, av, avl, m);
 642 }
 643
 644 /* --- @mpx_usqr@ --- *
 645  *
 646  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 647  *              @const mpw *av, *av@ = source vector base and limit
 648  *
 649  * Returns:     ---
 650  *
 651  * Use:         Performs unsigned integer squaring.  The result vector must
 652  *              not overlap the source vector in any way.
 653  */
 654
 655 void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
 656 {
 657   MPX_ZERO(dv, dvl);
 658
 659   /* --- Main loop --- */
 660
 661   while (av < avl) {
 662     const mpw *avv = av;
 663     mpw *dvv = dv;
 664     mpw a = *av;
 665     mpd c;
 666
 667     /* --- Stop if I've run out of destination --- */
 668
 669     if (dvv >= dvl)
 670       break;
 671
 672     /* --- Work out the square at this point in the proceedings --- */
 673
 674     {
 675       mpd x = (mpd)a * (mpd)a + *dvv;
 676       *dvv++ = MPW(x);
 677       c = MPW(x >> MPW_BITS);
 678     }
 679
 680     /* --- Now fix up the rest of the vector upwards --- */
 681
 682     avv++;
 683     while (dvv < dvl && avv < avl) {
 684       mpd x = (mpd)a * (mpd)*avv++;
 685       mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
 686       c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
 687       *dvv++ = MPW(y);
 688     }
 689     while (dvv < dvl && c) {
 690       mpd x = c + *dvv;
 691       *dvv++ = MPW(x);
 692       c = x >> MPW_BITS;
 693     }
 694
 695     /* --- Get ready for the next round --- */
 696
 697     av++;
 698     dv += 2;
 699   }
 700 }
 701
 702 /* --- @mpx_udiv@ --- *
 703  *
 704  * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
 705  *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
 706  *              @const mpw *dv, *dvl@ = divisor vector base and limit
 707  *              @mpw *sv, *svl@ = scratch workspace
 708  *
 709  * Returns:     ---
 710  *
 711  * Use:         Performs unsigned integer division.  If the result overflows
 712  *              the quotient vector, high-order bits are discarded.  (Clearly
 713  *              the remainder vector can't overflow.)  The various vectors
 714  *              may not overlap in any way.  Yes, I know it's a bit odd
 715  *              requiring the dividend to be in the result position but it
 716  *              does make some sense really.  The remainder must have
 717  *              headroom for at least two extra words.  The scratch space
 718  *              must be at least one word larger than the divisor.
 719  */
 720
 721 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
 722               const mpw *dv, const mpw *dvl,
 723               mpw *sv, mpw *svl)
 724 {
 725   unsigned norm = 0;
 726   size_t scale;
 727   mpw d, dd;
 728
 729   /* --- Initialize the quotient --- */
 730
 731   MPX_ZERO(qv, qvl);
 732
 733   /* --- Perform some sanity checks --- */
 734
 735   MPX_SHRINK(dv, dvl);
 736   assert(((void)"division by zero in mpx_udiv", dv < dvl));
 737
 738   /* --- Normalize the divisor --- *
 739    *
 740    * The algorithm requires that the divisor be at least two digits long.
 741    * This is easy to fix.
 742    */
 743
 744   {
 745     unsigned b;
 746
 747     d = dvl[-1];
 748     for (b = MPW_BITS / 2; b; b >>= 1) {
 749       if (d < (MPW_MAX >> b)) {
 750         d <<= b;
 751         norm += b;
 752       }
 753     }
 754     if (dv + 1 == dvl)
 755       norm += MPW_BITS;
 756   }
 757
 758   /* --- Normalize the dividend/remainder to match --- */
 759
 760   if (norm) {
 761     mpx_lsl(rv, rvl, rv, rvl, norm);
 762     mpx_lsl(sv, svl, dv, dvl, norm);
 763     dv = sv;
 764     dvl = svl;
 765     MPX_SHRINK(dv, dvl);
 766   }
 767
 768   MPX_SHRINK(rv, rvl);
 769   d = dvl[-1];
 770   dd = dvl[-2];
 771
 772   /* --- Work out the relative scales --- */
 773
 774   {
 775     size_t rvn = rvl - rv;
 776     size_t dvn = dvl - dv;
 777
 778     /* --- If the divisor is clearly larger, notice this --- */
 779
 780     if (dvn > rvn) {
 781       mpx_lsr(rv, rvl, rv, rvl, norm);
 782       return;
 783     }
 784
 785     scale = rvn - dvn;
 786   }
 787
 788   /* --- Calculate the most significant quotient digit --- *
 789    *
 790    * Because the divisor has its top bit set, this can only happen once.  The
 791    * pointer arithmetic is a little contorted, to make sure that the
 792    * behaviour is defined.
 793    */
 794
 795   if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
 796     mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
 797     if (qvl - qv > scale)
 798       qv[scale] = 1;
 799   }
 800
 801   /* --- Now for the main loop --- */
 802
 803   {
 804     mpw *rvv = rvl - 2;
 805
 806     while (scale) {
 807       mpw q;
 808       mpd rh;
 809
 810       /* --- Get an estimate for the next quotient digit --- */
 811
 812       mpw r = rvv[1];
 813       mpw rr = rvv[0];
 814       mpw rrr = *--rvv;
 815
 816       scale--;
 817       rh = ((mpd)r << MPW_BITS) | rr;
 818       if (r == d)
 819         q = MPW_MAX;
 820       else
 821         q = MPW(rh / d);
 822
 823       /* --- Refine the estimate --- */
 824
 825       {
 826         mpd yh = (mpd)d * q;
 827         mpd yl = (mpd)dd * q;
 828
 829         if (yl > MPW_MAX) {
 830           yh += yl >> MPW_BITS;
 831           yl &= MPW_MAX;
 832         }
 833
 834         while (yh > rh || (yh == rh && yl > rrr)) {
 835           q--;
 836           yh -= d;
 837           if (yl < dd) {
 838             yh++;
 839             yl += MPW_MAX;
 840           }
 841           yl -= dd;
 842         }
 843       }
 844
 845       /* --- Remove a chunk from the dividend --- */
 846
 847       {
 848         mpw *svv;
 849         const mpw *dvv;
 850         mpw mc = 0, sc = 0;
 851
 852         /* --- Calculate the size of the chunk --- *
 853          *
 854          * This does the whole job of calculating @r >> scale - qd@.
 855          */
 856
 857         for (svv = rv + scale, dvv = dv;
 858              dvv < dvl && svv < rvl;
 859              svv++, dvv++) {
 860           mpd x = (mpd)*dvv * (mpd)q + mc;
 861           mc = x >> MPW_BITS;
 862           x = (mpd)*svv - MPW(x) - sc;
 863           *svv = MPW(x);
 864           if (x >> MPW_BITS)
 865             sc = 1;
 866           else
 867             sc = 0;
 868         }
 869
 870         if (svv < rvl) {
 871           mpd x = (mpd)*svv - mc - sc;
 872           *svv++ = MPW(x);
 873           if (x >> MPW_BITS)
 874             sc = MPW_MAX;
 875           else
 876             sc = 0;
 877           while (svv < rvl)
 878             *svv++ = sc;
 879         }
 880
 881         /* --- Fix if the quotient was too large --- *
 882          *
 883          * This doesn't seem to happen very often.
 884          */
 885
 886         if (rvl[-1] > MPW_MAX / 2) {
 887           mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
 888           q--;
 889         }
 890       }
 891
 892       /* --- Done for another iteration --- */
 893
 894       if (qvl - qv > scale)
 895         qv[scale] = q;
 896       r = rr;
 897       rr = rrr;
 898     }
 899   }
 900
 901   /* --- Now fiddle with unnormalizing and things --- */
 902
 903   mpx_lsr(rv, rvl, rv, rvl, norm);
 904 }
 905
 906 /* --- @mpx_udivn@ --- *
 907  *
 908  * Arguments:   @mpw *qv, *qvl@ = storage for the quotient (may overlap
 909  *                      dividend)
 910  *              @const mpw *rv, *rvl@ = dividend
 911  *              @mpw d@ = single-precision divisor
 912  *
 913  * Returns:     Remainder after divison.
 914  *
 915  * Use:         Performs a single-precision division operation.
 916  */
 917
 918 mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
 919 {
 920   size_t i;
 921   size_t ql = qvl - qv;
 922   mpd r = 0;
 923
 924   i = rvl - rv;
 925   while (i > 0) {
 926     i--;
 927     r = (r << MPW_BITS) | rv[i];
 928     if (i < ql)
 929       qv[i] = r / d;
 930     r %= d;
 931   }
 932   return (MPW(r));
 933 }
 934
 935 /*----- Test rig ----------------------------------------------------------*/
 936
 937 #ifdef TEST_RIG
 938
 939 #include <mLib/alloc.h>
 940 #include <mLib/dstr.h>
 941 #include <mLib/quis.h>
 942 #include <mLib/testrig.h>
 943
 944 #include "mpscan.h"
 945
 946 #define ALLOC(v, vl, sz) do {                                           \
 947   size_t _sz = (sz);                                                    \
 948   mpw *_vv = xmalloc(MPWS(_sz));                                        \
 949   mpw *_vvl = _vv + _sz;                                                \
 950   (v) = _vv;                                                            \
 951   (vl) = _vvl;                                                          \
 952 } while (0)
 953
 954 #define LOAD(v, vl, d) do {                                             \
 955   const dstr *_d = (d);                                                 \
 956   mpw *_v, *_vl;                                                        \
 957   ALLOC(_v, _vl, MPW_RQ(_d->len));                                      \
 958   mpx_loadb(_v, _vl, _d->buf, _d->len);                                 \
 959   (v) = _v;                                                             \
 960   (vl) = _vl;                                                           \
 961 } while (0)
 962
 963 #define MAX(x, y) ((x) > (y) ? (x) : (y))
 964
 965 static void dumpbits(const char *msg, const void *pp, size_t sz)
 966 {
 967   const octet *p = pp;
 968   fputs(msg, stderr);
 969   for (; sz; sz--)
 970     fprintf(stderr, " %02x", *p++);
 971   fputc('\n', stderr);
 972 }
 973
 974 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
 975 {
 976   fputs(msg, stderr);
 977   MPX_SHRINK(v, vl);
 978   while (v < vl)
 979     fprintf(stderr, " %08lx", (unsigned long)*--vl);
 980   fputc('\n', stderr);
 981 }
 982
 983 static int chkscan(const mpw *v, const mpw *vl,
 984                    const void *pp, size_t sz, int step)
 985 {
 986   mpscan mps;
 987   const octet *p = pp;
 988   unsigned bit = 0;
 989   int ok = 1;
 990
 991   mpscan_initx(&mps, v, vl);
 992   while (sz) {
 993     unsigned x = *p;
 994     int i;
 995     p += step;
 996     for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
 997       if (MPSCAN_BIT(&mps) != (x & 1)) {
 998         fprintf(stderr,
 999                 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1000                 step, bit, x & 1, MPSCAN_BIT(&mps));
1001         ok = 0;
1002       }
1003       x >>= 1;
1004       bit++;
1005     }
1006     sz--;
1007   }
1008
1009   return (ok);
1010 }
1011
1012 static int loadstore(dstr *v)
1013 {
1014   dstr d = DSTR_INIT;
1015   size_t sz = MPW_RQ(v->len) * 2, diff;
1016   mpw *m, *ml;
1017   int ok = 1;
1018
1019   dstr_ensure(&d, v->len);
1020   m = xmalloc(MPWS(sz));
1021
1022   for (diff = 0; diff < sz; diff += 5) {
1023     size_t oct;
1024
1025     ml = m + sz - diff;
1026
1027     mpx_loadl(m, ml, v->buf, v->len);
1028     if (!chkscan(m, ml, v->buf, v->len, +1))
1029       ok = 0;
1030     MPX_OCTETS(oct, m, ml);
1031     mpx_storel(m, ml, d.buf, d.sz);
1032     if (memcmp(d.buf, v->buf, oct) != 0) {
1033       dumpbits("\n*** storel failed", d.buf, d.sz);
1034       ok = 0;
1035     }
1036
1037     mpx_loadb(m, ml, v->buf, v->len);
1038     if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1039       ok = 0;
1040     MPX_OCTETS(oct, m, ml);
1041     mpx_storeb(m, ml, d.buf, d.sz);
1042     if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1043       dumpbits("\n*** storeb failed", d.buf, d.sz);
1044       ok = 0;
1045     }
1046   }
1047
1048   if (!ok)
1049     dumpbits("input data", v->buf, v->len);
1050
1051   free(m);
1052   dstr_destroy(&d);
1053   return (ok);
1054 }
1055
1056 static int lsl(dstr *v)
1057 {
1058   mpw *a, *al;
1059   int n = *(int *)v[1].buf;
1060   mpw *c, *cl;
1061   mpw *d, *dl;
1062   int ok = 1;
1063
1064   LOAD(a, al, &v[0]);
1065   LOAD(c, cl, &v[2]);
1066   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1067
1068   mpx_lsl(d, dl, a, al, n);
1069   if (MPX_UCMP(d, dl, !=, c, cl)) {
1070     fprintf(stderr, "\n*** lsl(%i) failed\n", n);
1071     dumpmp("       a", a, al);
1072     dumpmp("expected", c, cl);
1073     dumpmp("  result", d, dl);
1074     ok = 0;
1075   }
1076
1077   free(a); free(c); free(d);
1078   return (ok);
1079 }
1080
1081 static int lsr(dstr *v)
1082 {
1083   mpw *a, *al;
1084   int n = *(int *)v[1].buf;
1085   mpw *c, *cl;
1086   mpw *d, *dl;
1087   int ok = 1;
1088
1089   LOAD(a, al, &v[0]);
1090   LOAD(c, cl, &v[2]);
1091   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1092
1093   mpx_lsr(d, dl, a, al, n);
1094   if (MPX_UCMP(d, dl, !=, c, cl)) {
1095     fprintf(stderr, "\n*** lsr(%i) failed\n", n);
1096     dumpmp("       a", a, al);
1097     dumpmp("expected", c, cl);
1098     dumpmp("  result", d, dl);
1099     ok = 0;
1100   }
1101
1102   free(a); free(c); free(d);
1103   return (ok);
1104 }
1105
1106 static int uadd(dstr *v)
1107 {
1108   mpw *a, *al;
1109   mpw *b, *bl;
1110   mpw *c, *cl;
1111   mpw *d, *dl;
1112   int ok = 1;
1113
1114   LOAD(a, al, &v[0]);
1115   LOAD(b, bl, &v[1]);
1116   LOAD(c, cl, &v[2]);
1117   ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1118
1119   mpx_uadd(d, dl, a, al, b, bl);
1120   if (MPX_UCMP(d, dl, !=, c, cl)) {
1121     fprintf(stderr, "\n*** uadd failed\n");
1122     dumpmp("       a", a, al);
1123     dumpmp("       b", b, bl);
1124     dumpmp("expected", c, cl);
1125     dumpmp("  result", d, dl);
1126     ok = 0;
1127   }
1128
1129   free(a); free(b); free(c); free(d);
1130   return (ok);
1131 }
1132
1133 static int usub(dstr *v)
1134 {
1135   mpw *a, *al;
1136   mpw *b, *bl;
1137   mpw *c, *cl;
1138   mpw *d, *dl;
1139   int ok = 1;
1140
1141   LOAD(a, al, &v[0]);
1142   LOAD(b, bl, &v[1]);
1143   LOAD(c, cl, &v[2]);
1144   ALLOC(d, dl, al - a);
1145
1146   mpx_usub(d, dl, a, al, b, bl);
1147   if (MPX_UCMP(d, dl, !=, c, cl)) {
1148     fprintf(stderr, "\n*** usub failed\n");
1149     dumpmp("       a", a, al);
1150     dumpmp("       b", b, bl);
1151     dumpmp("expected", c, cl);
1152     dumpmp("  result", d, dl);
1153     ok = 0;
1154   }
1155
1156   free(a); free(b); free(c); free(d);
1157   return (ok);
1158 }
1159
1160 static int umul(dstr *v)
1161 {
1162   mpw *a, *al;
1163   mpw *b, *bl;
1164   mpw *c, *cl;
1165   mpw *d, *dl;
1166   int ok = 1;
1167
1168   LOAD(a, al, &v[0]);
1169   LOAD(b, bl, &v[1]);
1170   LOAD(c, cl, &v[2]);
1171   ALLOC(d, dl, (al - a) + (bl - b));
1172
1173   mpx_umul(d, dl, a, al, b, bl);
1174   if (MPX_UCMP(d, dl, !=, c, cl)) {
1175     fprintf(stderr, "\n*** umul failed\n");
1176     dumpmp("       a", a, al);
1177     dumpmp("       b", b, bl);
1178     dumpmp("expected", c, cl);
1179     dumpmp("  result", d, dl);
1180     ok = 0;
1181   }
1182
1183   free(a); free(b); free(c); free(d);
1184   return (ok);
1185 }
1186
1187 static int usqr(dstr *v)
1188 {
1189   mpw *a, *al;
1190   mpw *c, *cl;
1191   mpw *d, *dl;
1192   int ok = 1;
1193
1194   LOAD(a, al, &v[0]);
1195   LOAD(c, cl, &v[1]);
1196   ALLOC(d, dl, 2 * (al - a));
1197
1198   mpx_usqr(d, dl, a, al);
1199   if (MPX_UCMP(d, dl, !=, c, cl)) {
1200     fprintf(stderr, "\n*** usqr failed\n");
1201     dumpmp("       a", a, al);
1202     dumpmp("expected", c, cl);
1203     dumpmp("  result", d, dl);
1204     ok = 0;
1205   }
1206
1207   free(a); free(c); free(d);
1208   return (ok);
1209 }
1210
1211 static int udiv(dstr *v)
1212 {
1213   mpw *a, *al;
1214   mpw *b, *bl;
1215   mpw *q, *ql;
1216   mpw *r, *rl;
1217   mpw *qq, *qql;
1218   mpw *s, *sl;
1219   int ok = 1;
1220
1221   ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1222   LOAD(b, bl, &v[1]);
1223   LOAD(q, ql, &v[2]);
1224   LOAD(r, rl, &v[3]);
1225   ALLOC(qq, qql, al - a);
1226   ALLOC(s, sl, (bl - b) + 1);
1227
1228   mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1229   if (MPX_UCMP(qq, qql, !=, q, ql) ||
1230       MPX_UCMP(a, al, !=, r, rl)) {
1231     fprintf(stderr, "\n*** udiv failed\n");
1232     dumpmp(" divisor", b, bl);
1233     dumpmp("expect r", r, rl);
1234     dumpmp("result r", a, al);
1235     dumpmp("expect q", q, ql);
1236     dumpmp("result q", qq, qql);
1237     ok = 0;
1238   }
1239
1240   free(a); free(b); free(r); free(q); free(s); free(qq);
1241   return (ok);
1242 }
1243
1244 static test_chunk defs[] = {
1245   { "load-store", loadstore, { &type_hex, 0 } },
1246   { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
1247   { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1248   { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1249   { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1250   { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1251   { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1252   { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1253   { 0, 0, { 0 } }
1254 };
1255
1256 int main(int argc, char *argv[])
1257 {
1258   test_run(argc, argv, defs, SRCDIR"/tests/mpx");
1259   return (0);
1260 }
1261
1262
1263 #endif
1264
1265 /*----- That's all, folks -------------------------------------------------*/