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

   1 /* -*-c-*-
   2  *
   3  * Low-level multiprecision arithmetic
   4  *
   5  * (c) 1999 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include <assert.h>
  31 #include <stdio.h>
  32 #include <stdlib.h>
  33 #include <string.h>
  34
  35 #include <mLib/bits.h>
  36
  37 #include "mptypes.h"
  38 #include "mpx.h"
  39 #include "bitops.h"
  40
  41 /*----- Loading and storing -----------------------------------------------*/
  42
  43 /* --- @mpx_storel@ --- *
  44  *
  45  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
  46  *              @void *pp@ = pointer to octet array
  47  *              @size_t sz@ = size of octet array
  48  *
  49  * Returns:     ---
  50  *
  51  * Use:         Stores an MP in an octet array, least significant octet
  52  *              first.  High-end octets are silently discarded if there
  53  *              isn't enough space for them.
  54  */
  55
  56 void mpx_storel(const mpw *v, const mpw *vl, void *pp, size_t sz)
  57 {
  58   mpw n, w = 0;
  59   octet *p = pp, *q = p + sz;
  60   unsigned bits = 0;
  61
  62   while (p < q) {
  63     if (bits < 8) {
  64       if (v >= vl) {
  65         *p++ = U8(w);
  66         break;
  67       }
  68       n = *v++;
  69       *p++ = U8(w | n << bits);
  70       w = n >> (8 - bits);
  71       bits += MPW_BITS - 8;
  72     } else {
  73       *p++ = U8(w);
  74       w >>= 8;
  75       bits -= 8;
  76     }
  77   }
  78   memset(p, 0, q - p);
  79 }
  80
  81 /* --- @mpx_loadl@ --- *
  82  *
  83  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
  84  *              @const void *pp@ = pointer to octet array
  85  *              @size_t sz@ = size of octet array
  86  *
  87  * Returns:     ---
  88  *
  89  * Use:         Loads an MP in an octet array, least significant octet
  90  *              first.  High-end octets are ignored if there isn't enough
  91  *              space for them.
  92  */
  93
  94 void mpx_loadl(mpw *v, mpw *vl, const void *pp, size_t sz)
  95 {
  96   unsigned n;
  97   mpw w = 0;
  98   const octet *p = pp, *q = p + sz;
  99   unsigned bits = 0;
 100
 101   if (v >= vl)
 102     return;
 103   while (p < q) {
 104     n = U8(*p++);
 105     w |= n << bits;
 106     bits += 8;
 107     if (bits >= MPW_BITS) {
 108       *v++ = MPW(w);
 109       w = n >> (MPW_BITS - bits + 8);
 110       bits -= MPW_BITS;
 111       if (v >= vl)
 112         return;
 113     }
 114   }
 115   *v++ = w;
 116   MPX_ZERO(v, vl);
 117 }
 118
 119 /* --- @mpx_storeb@ --- *
 120  *
 121  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 122  *              @void *pp@ = pointer to octet array
 123  *              @size_t sz@ = size of octet array
 124  *
 125  * Returns:     ---
 126  *
 127  * Use:         Stores an MP in an octet array, most significant octet
 128  *              first.  High-end octets are silently discarded if there
 129  *              isn't enough space for them.
 130  */
 131
 132 void mpx_storeb(const mpw *v, const mpw *vl, void *pp, size_t sz)
 133 {
 134   mpw n, w = 0;
 135   octet *p = pp, *q = p + sz;
 136   unsigned bits = 0;
 137
 138   while (q > p) {
 139     if (bits < 8) {
 140       if (v >= vl) {
 141         *--q = U8(w);
 142         break;
 143       }
 144       n = *v++;
 145       *--q = U8(w | n << bits);
 146       w = n >> (8 - bits);
 147       bits += MPW_BITS - 8;
 148     } else {
 149       *--q = U8(w);
 150       w >>= 8;
 151       bits -= 8;
 152     }
 153   }
 154   memset(p, 0, q - p);
 155 }
 156
 157 /* --- @mpx_loadb@ --- *
 158  *
 159  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 160  *              @const void *pp@ = pointer to octet array
 161  *              @size_t sz@ = size of octet array
 162  *
 163  * Returns:     ---
 164  *
 165  * Use:         Loads an MP in an octet array, most significant octet
 166  *              first.  High-end octets are ignored if there isn't enough
 167  *              space for them.
 168  */
 169
 170 void mpx_loadb(mpw *v, mpw *vl, const void *pp, size_t sz)
 171 {
 172   unsigned n;
 173   mpw w = 0;
 174   const octet *p = pp, *q = p + sz;
 175   unsigned bits = 0;
 176
 177   if (v >= vl)
 178     return;
 179   while (q > p) {
 180     n = U8(*--q);
 181     w |= n << bits;
 182     bits += 8;
 183     if (bits >= MPW_BITS) {
 184       *v++ = MPW(w);
 185       w = n >> (MPW_BITS - bits + 8);
 186       bits -= MPW_BITS;
 187       if (v >= vl)
 188         return;
 189     }
 190   }
 191   *v++ = w;
 192   MPX_ZERO(v, vl);
 193 }
 194
 195 /* --- @mpx_storel2cn@ --- *
 196  *
 197  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 198  *              @void *pp@ = pointer to octet array
 199  *              @size_t sz@ = size of octet array
 200  *
 201  * Returns:     ---
 202  *
 203  * Use:         Stores a negative MP in an octet array, least significant
 204  *              octet first, as two's complement.  High-end octets are
 205  *              silently discarded if there isn't enough space for them.
 206  *              This obviously makes the output bad.
 207  */
 208
 209 void mpx_storel2cn(const mpw *v, const mpw *vl, void *pp, size_t sz)
 210 {
 211   unsigned c = 1;
 212   unsigned b = 0;
 213   mpw n, w = 0;
 214   octet *p = pp, *q = p + sz;
 215   unsigned bits = 0;
 216
 217   while (p < q) {
 218     if (bits < 8) {
 219       if (v >= vl) {
 220         b = w;
 221         break;
 222       }
 223       n = *v++;
 224       b = w | n << bits;
 225       w = n >> (8 - bits);
 226       bits += MPW_BITS - 8;
 227     } else {
 228       b = w;
 229       w >>= 8;
 230       bits -= 8;
 231     }
 232     b = U8(~b + c);
 233     c = c && !b;
 234     *p++ = b;
 235   }
 236   while (p < q) {
 237     b = U8(~b + c);
 238     c = c && !b;
 239     *p++ = b;
 240     b = 0;
 241   }
 242 }
 243
 244 /* --- @mpx_loadl2cn@ --- *
 245  *
 246  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 247  *              @const void *pp@ = pointer to octet array
 248  *              @size_t sz@ = size of octet array
 249  *
 250  * Returns:     ---
 251  *
 252  * Use:         Loads a negative MP in an octet array, least significant
 253  *              octet first, as two's complement.  High-end octets are
 254  *              ignored if there isn't enough space for them.  This probably
 255  *              means you made the wrong choice coming here.
 256  */
 257
 258 void mpx_loadl2cn(mpw *v, mpw *vl, const void *pp, size_t sz)
 259 {
 260   unsigned n;
 261   unsigned c = 1;
 262   mpw w = 0;
 263   const octet *p = pp, *q = p + sz;
 264   unsigned bits = 0;
 265
 266   if (v >= vl)
 267     return;
 268   while (p < q) {
 269     n = U8(~(*p++) + c);
 270     c = c && !n;
 271     w |= n << bits;
 272     bits += 8;
 273     if (bits >= MPW_BITS) {
 274       *v++ = MPW(w);
 275       w = n >> (MPW_BITS - bits + 8);
 276       bits -= MPW_BITS;
 277       if (v >= vl)
 278         return;
 279     }
 280   }
 281   *v++ = w;
 282   MPX_ZERO(v, vl);
 283 }
 284
 285 /* --- @mpx_storeb2cn@ --- *
 286  *
 287  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 288  *              @void *pp@ = pointer to octet array
 289  *              @size_t sz@ = size of octet array
 290  *
 291  * Returns:     ---
 292  *
 293  * Use:         Stores a negative MP in an octet array, most significant
 294  *              octet first, as two's complement.  High-end octets are
 295  *              silently discarded if there isn't enough space for them,
 296  *              which probably isn't what you meant.
 297  */
 298
 299 void mpx_storeb2cn(const mpw *v, const mpw *vl, void *pp, size_t sz)
 300 {
 301   mpw n, w = 0;
 302   unsigned b = 0;
 303   unsigned c = 1;
 304   octet *p = pp, *q = p + sz;
 305   unsigned bits = 0;
 306
 307   while (q > p) {
 308     if (bits < 8) {
 309       if (v >= vl) {
 310         b = w;
 311         break;
 312       }
 313       n = *v++;
 314       b = w | n << bits;
 315       w = n >> (8 - bits);
 316       bits += MPW_BITS - 8;
 317     } else {
 318       b = w;
 319       w >>= 8;
 320       bits -= 8;
 321     }
 322     b = U8(~b + c);
 323     c = c && !b;
 324     *--q = b;
 325   }
 326   while (q > p) {
 327     b = ~b + c;
 328     c = c && !(b & 0xff);
 329     *--q = b;
 330     b = 0;
 331   }
 332 }
 333
 334 /* --- @mpx_loadb2cn@ --- *
 335  *
 336  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 337  *              @const void *pp@ = pointer to octet array
 338  *              @size_t sz@ = size of octet array
 339  *
 340  * Returns:     ---
 341  *
 342  * Use:         Loads a negative MP in an octet array, most significant octet
 343  *              first as two's complement.  High-end octets are ignored if
 344  *              there isn't enough space for them.  This probably means you
 345  *              chose this function wrongly.
 346  */
 347
 348 void mpx_loadb2cn(mpw *v, mpw *vl, const void *pp, size_t sz)
 349 {
 350   unsigned n;
 351   unsigned c = 1;
 352   mpw w = 0;
 353   const octet *p = pp, *q = p + sz;
 354   unsigned bits = 0;
 355
 356   if (v >= vl)
 357     return;
 358   while (q > p) {
 359     n = U8(~(*--q) + c);
 360     c = c && !n;
 361     w |= n << bits;
 362     bits += 8;
 363     if (bits >= MPW_BITS) {
 364       *v++ = MPW(w);
 365       w = n >> (MPW_BITS - bits + 8);
 366       bits -= MPW_BITS;
 367       if (v >= vl)
 368         return;
 369     }
 370   }
 371   *v++ = w;
 372   MPX_ZERO(v, vl);
 373 }
 374
 375 /*----- Logical shifting --------------------------------------------------*/
 376
 377 /* --- @mpx_lsl@ --- *
 378  *
 379  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 380  *              @const mpw *av, *avl@ = source vector base and limit
 381  *              @size_t n@ = number of bit positions to shift by
 382  *
 383  * Returns:     ---
 384  *
 385  * Use:         Performs a logical shift left operation on an integer.
 386  */
 387
 388 void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 389 {
 390   size_t nw;
 391   unsigned nb;
 392
 393   /* --- Trivial special case --- */
 394
 395   if (n == 0)
 396     MPX_COPY(dv, dvl, av, avl);
 397
 398   /* --- Single bit shifting --- */
 399
 400   else if (n == 1) {
 401     mpw w = 0;
 402     while (av < avl) {
 403       mpw t;
 404       if (dv >= dvl)
 405         goto done;
 406       t = *av++;
 407       *dv++ = MPW((t << 1) | w);
 408       w = t >> (MPW_BITS - 1);
 409     }
 410     if (dv >= dvl)
 411       goto done;
 412     *dv++ = MPW(w);
 413     MPX_ZERO(dv, dvl);
 414     goto done;
 415   }
 416
 417   /* --- Break out word and bit shifts for more sophisticated work --- */
 418
 419   nw = n / MPW_BITS;
 420   nb = n % MPW_BITS;
 421
 422   /* --- Handle a shift by a multiple of the word size --- */
 423
 424   if (nb == 0) {
 425     if (nw >= dvl - dv)
 426       MPX_ZERO(dv, dvl);
 427     else {
 428       MPX_COPY(dv + nw, dvl, av, avl);
 429       memset(dv, 0, MPWS(nw));
 430     }
 431   }
 432
 433   /* --- And finally the difficult case --- *
 434    *
 435    * This is a little convoluted, because I have to start from the end and
 436    * work backwards to avoid overwriting the source, if they're both the same
 437    * block of memory.
 438    */
 439
 440   else {
 441     mpw w;
 442     size_t nr = MPW_BITS - nb;
 443     size_t dvn = dvl - dv;
 444     size_t avn = avl - av;
 445
 446     if (dvn <= nw) {
 447       MPX_ZERO(dv, dvl);
 448       goto done;
 449     }
 450
 451     if (dvn > avn + nw) {
 452       size_t off = avn + nw + 1;
 453       MPX_ZERO(dv + off, dvl);
 454       dvl = dv + off;
 455       w = 0;
 456     } else {
 457       avl = av + dvn - nw;
 458       w = *--avl << nb;
 459     }
 460
 461     while (avl > av) {
 462       mpw t = *--avl;
 463       *--dvl = MPW((t >> nr) | w);
 464       w = t << nb;
 465     }
 466
 467     *--dvl = MPW(w);
 468     MPX_ZERO(dv, dvl);
 469   }
 470
 471 done:;
 472 }
 473
 474 /* --- @mpx_lslc@ --- *
 475  *
 476  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 477  *              @const mpw *av, *avl@ = source vector base and limit
 478  *              @size_t n@ = number of bit positions to shift by
 479  *
 480  * Returns:     ---
 481  *
 482  * Use:         Performs a logical shift left operation on an integer, only
 483  *              it fills in the bits with ones instead of zeroes.
 484  */
 485
 486 void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 487 {
 488   size_t nw;
 489   unsigned nb;
 490
 491   /* --- Trivial special case --- */
 492
 493   if (n == 0)
 494     MPX_COPY(dv, dvl, av, avl);
 495
 496   /* --- Single bit shifting --- */
 497
 498   else if (n == 1) {
 499     mpw w = 1;
 500     while (av < avl) {
 501       mpw t;
 502       if (dv >= dvl)
 503         goto done;
 504       t = *av++;
 505       *dv++ = MPW((t << 1) | w);
 506       w = t >> (MPW_BITS - 1);
 507     }
 508     if (dv >= dvl)
 509       goto done;
 510     *dv++ = MPW(w);
 511     MPX_ZERO(dv, dvl);
 512     goto done;
 513   }
 514
 515   /* --- Break out word and bit shifts for more sophisticated work --- */
 516
 517   nw = n / MPW_BITS;
 518   nb = n % MPW_BITS;
 519
 520   /* --- Handle a shift by a multiple of the word size --- */
 521
 522   if (nb == 0) {
 523     if (nw >= dvl - dv)
 524       MPX_ONE(dv, dvl);
 525     else {
 526       MPX_COPY(dv + nw, dvl, av, avl);
 527       MPX_ONE(dv, dv + nw);
 528     }
 529   }
 530
 531   /* --- And finally the difficult case --- *
 532    *
 533    * This is a little convoluted, because I have to start from the end and
 534    * work backwards to avoid overwriting the source, if they're both the same
 535    * block of memory.
 536    */
 537
 538   else {
 539     mpw w;
 540     size_t nr = MPW_BITS - nb;
 541     size_t dvn = dvl - dv;
 542     size_t avn = avl - av;
 543
 544     if (dvn <= nw) {
 545       MPX_ONE(dv, dvl);
 546       goto done;
 547     }
 548
 549     if (dvn > avn + nw) {
 550       size_t off = avn + nw + 1;
 551       MPX_ZERO(dv + off, dvl);
 552       dvl = dv + off;
 553       w = 0;
 554     } else {
 555       avl = av + dvn - nw;
 556       w = *--avl << nb;
 557     }
 558
 559     while (avl > av) {
 560       mpw t = *--avl;
 561       *--dvl = MPW((t >> nr) | w);
 562       w = t << nb;
 563     }
 564
 565     *--dvl = MPW((MPW_MAX >> nr) | w);
 566     MPX_ONE(dv, dvl);
 567   }
 568
 569 done:;
 570 }
 571
 572 /* --- @mpx_lsr@ --- *
 573  *
 574  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 575  *              @const mpw *av, *avl@ = source vector base and limit
 576  *              @size_t n@ = number of bit positions to shift by
 577  *
 578  * Returns:     ---
 579  *
 580  * Use:         Performs a logical shift right operation on an integer.
 581  */
 582
 583 void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 584 {
 585   size_t nw;
 586   unsigned nb;
 587
 588   /* --- Trivial special case --- */
 589
 590   if (n == 0)
 591     MPX_COPY(dv, dvl, av, avl);
 592
 593   /* --- Single bit shifting --- */
 594
 595   else if (n == 1) {
 596     mpw w = av < avl ? *av++ >> 1 : 0;
 597     while (av < avl) {
 598       mpw t;
 599       if (dv >= dvl)
 600         goto done;
 601       t = *av++;
 602       *dv++ = MPW((t << (MPW_BITS - 1)) | w);
 603       w = t >> 1;
 604     }
 605     if (dv >= dvl)
 606       goto done;
 607     *dv++ = MPW(w);
 608     MPX_ZERO(dv, dvl);
 609     goto done;
 610   }
 611
 612   /* --- Break out word and bit shifts for more sophisticated work --- */
 613
 614   nw = n / MPW_BITS;
 615   nb = n % MPW_BITS;
 616
 617   /* --- Handle a shift by a multiple of the word size --- */
 618
 619   if (nb == 0) {
 620     if (nw >= avl - av)
 621       MPX_ZERO(dv, dvl);
 622     else
 623       MPX_COPY(dv, dvl, av + nw, avl);
 624   }
 625
 626   /* --- And finally the difficult case --- */
 627
 628   else {
 629     mpw w;
 630     size_t nr = MPW_BITS - nb;
 631
 632     av += nw;
 633     w = av < avl ? *av++ : 0;
 634     while (av < avl) {
 635       mpw t;
 636       if (dv >= dvl)
 637         goto done;
 638       t = *av++;
 639       *dv++ = MPW((w >> nb) | (t << nr));
 640       w = t;
 641     }
 642     if (dv < dvl) {
 643       *dv++ = MPW(w >> nb);
 644       MPX_ZERO(dv, dvl);
 645     }
 646   }
 647
 648 done:;
 649 }
 650
 651 /*----- Bitwise operations ------------------------------------------------*/
 652
 653 /* --- @mpx_bitop@ --- *
 654  *
 655  * Arguments:   @mpw *dv, *dvl@ = destination vector
 656  *              @const mpw *av, *avl@ = first source vector
 657  *              @const mpw *bv, *bvl@ = second source vector
 658  *
 659  * Returns:     ---
 660  *
 661  * Use;         Provides the dyadic boolean functions.
 662  */
 663
 664 /* GCC complains about the generated code, so try to silence it. */
 665 #if __GNUC__ >= 5 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
 666 #  pragma GCC diagnostic push
 667 #  pragma GCC diagnostic ignored "-Wunused-but-set-variable"
 668 #endif
 669
 670 #define MPX_BITBINOP(string)                                            \
 671                                                                         \
 672 void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,  \
 673                      const mpw *bv, const mpw *bvl)                     \
 674 {                                                                       \
 675   MPX_SHRINK(av, avl);                                                  \
 676   MPX_SHRINK(bv, bvl);                                                  \
 677                                                                         \
 678   while (dv < dvl) {                                                    \
 679     mpw a, b;                                                           \
 680     a = (av < avl) ? *av++ : 0;                                         \
 681     b = (bv < bvl) ? *bv++ : 0;                                         \
 682     *dv++ = B##string(a, b);                                            \
 683   }                                                                     \
 684 }
 685
 686 MPX_DOBIN(MPX_BITBINOP)
 687
 688 #if __GNUC__ >= 5 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)
 689 #  pragma GCC diagnostic pop
 690 #endif
 691
 692 void mpx_not(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
 693 {
 694   MPX_SHRINK(av, avl);
 695
 696   while (dv < dvl) {
 697     mpw a;
 698     a = (av < avl) ? *av++ : 0;
 699     *dv++ = ~a;
 700   }
 701 }
 702
 703 /*----- Unsigned arithmetic -----------------------------------------------*/
 704
 705 /* --- @mpx_2c@ --- *
 706  *
 707  * Arguments:   @mpw *dv, *dvl@ = destination vector
 708  *              @const mpw *v, *vl@ = source vector
 709  *
 710  * Returns:     ---
 711  *
 712  * Use:         Calculates the two's complement of @v@.
 713  */
 714
 715 void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
 716 {
 717   mpw c = 0;
 718   while (dv < dvl && v < vl)
 719     *dv++ = c = MPW(~*v++);
 720   if (dv < dvl) {
 721     if (c > MPW_MAX / 2)
 722       c = MPW(~0);
 723     while (dv < dvl)
 724       *dv++ = c;
 725   }
 726   MPX_UADDN(dv, dvl, 1);
 727 }
 728
 729 /* --- @mpx_ueq@ --- *
 730  *
 731  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 732  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 733  *
 734  * Returns:     Nonzero if the two vectors are equal.
 735  *
 736  * Use:         Performs an unsigned integer test for equality.
 737  */
 738
 739 int mpx_ueq(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
 740 {
 741   MPX_SHRINK(av, avl);
 742   MPX_SHRINK(bv, bvl);
 743   if (avl - av != bvl - bv)
 744     return (0);
 745   while (av < avl) {
 746     if (*av++ != *bv++)
 747       return (0);
 748   }
 749   return (1);
 750 }
 751
 752 /* --- @mpx_ucmp@ --- *
 753  *
 754  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 755  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 756  *
 757  * Returns:     Less than, equal to, or greater than zero depending on
 758  *              whether @a@ is less than, equal to or greater than @b@,
 759  *              respectively.
 760  *
 761  * Use:         Performs an unsigned integer comparison.
 762  */
 763
 764 int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
 765 {
 766   MPX_SHRINK(av, avl);
 767   MPX_SHRINK(bv, bvl);
 768
 769   if (avl - av > bvl - bv)
 770     return (+1);
 771   else if (avl - av < bvl - bv)
 772     return (-1);
 773   else while (avl > av) {
 774     mpw a = *--avl, b = *--bvl;
 775     if (a > b)
 776       return (+1);
 777     else if (a < b)
 778       return (-1);
 779   }
 780   return (0);
 781 }
 782
 783 /* --- @mpx_uadd@ --- *
 784  *
 785  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 786  *              @const mpw *av, *avl@ = first addend vector base and limit
 787  *              @const mpw *bv, *bvl@ = second addend vector base and limit
 788  *
 789  * Returns:     ---
 790  *
 791  * Use:         Performs unsigned integer addition.  If the result overflows
 792  *              the destination vector, high-order bits are discarded.  This
 793  *              means that two's complement addition happens more or less for
 794  *              free, although that's more a side-effect than anything else.
 795  *              The result vector may be equal to either or both source
 796  *              vectors, but may not otherwise overlap them.
 797  */
 798
 799 void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 800               const mpw *bv, const mpw *bvl)
 801 {
 802   mpw c = 0;
 803
 804   while (av < avl || bv < bvl) {
 805     mpw a, b;
 806     mpd x;
 807     if (dv >= dvl)
 808       return;
 809     a = (av < avl) ? *av++ : 0;
 810     b = (bv < bvl) ? *bv++ : 0;
 811     x = (mpd)a + (mpd)b + c;
 812     *dv++ = MPW(x);
 813     c = x >> MPW_BITS;
 814   }
 815   if (dv < dvl) {
 816     *dv++ = c;
 817     MPX_ZERO(dv, dvl);
 818   }
 819 }
 820
 821 /* --- @mpx_uaddn@ --- *
 822  *
 823  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 824  *              @mpw n@ = other addend
 825  *
 826  * Returns:     ---
 827  *
 828  * Use:         Adds a small integer to a multiprecision number.
 829  */
 830
 831 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
 832
 833 /* --- @mpx_uaddnlsl@ --- *
 834  *
 835  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 836  *              @mpw a@ = second argument
 837  *              @unsigned o@ = offset in bits
 838  *
 839  * Returns:     ---
 840  *
 841  * Use:         Computes %$d + 2^o a$%.  If the result overflows then
 842  *              high-order bits are discarded, as usual.  We must have
 843  *              @0 < o < MPW_BITS@.
 844  */
 845
 846 void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
 847 {
 848   mpd x = (mpd)a << o;
 849
 850   while (x && dv < dvl) {
 851     x += *dv;
 852     *dv++ = MPW(x);
 853     x >>= MPW_BITS;
 854   }
 855 }
 856
 857 /* --- @mpx_usub@ --- *
 858  *
 859  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 860  *              @const mpw *av, *avl@ = first argument vector base and limit
 861  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 862  *
 863  * Returns:     ---
 864  *
 865  * Use:         Performs unsigned integer subtraction.  If the result
 866  *              overflows the destination vector, high-order bits are
 867  *              discarded.  This means that two's complement subtraction
 868  *              happens more or less for free, althuogh that's more a side-
 869  *              effect than anything else.  The result vector may be equal to
 870  *              either or both source vectors, but may not otherwise overlap
 871  *              them.
 872  */
 873
 874 void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 875               const mpw *bv, const mpw *bvl)
 876 {
 877   mpw c = 0;
 878
 879   while (av < avl || bv < bvl) {
 880     mpw a, b;
 881     mpd x;
 882     if (dv >= dvl)
 883       return;
 884     a = (av < avl) ? *av++ : 0;
 885     b = (bv < bvl) ? *bv++ : 0;
 886     x = (mpd)a - (mpd)b - c;
 887     *dv++ = MPW(x);
 888     if (x >> MPW_BITS)
 889       c = 1;
 890     else
 891       c = 0;
 892   }
 893   if (c)
 894     c = MPW_MAX;
 895   while (dv < dvl)
 896     *dv++ = c;
 897 }
 898
 899 /* --- @mpx_usubn@ --- *
 900  *
 901  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 902  *              @n@ = subtrahend
 903  *
 904  * Returns:     ---
 905  *
 906  * Use:         Subtracts a small integer from a multiprecision number.
 907  */
 908
 909 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
 910
 911 /* --- @mpx_uaddnlsl@ --- *
 912  *
 913  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 914  *              @mpw a@ = second argument
 915  *              @unsigned o@ = offset in bits
 916  *
 917  * Returns:     ---
 918  *
 919  * Use:         Computes %$d + 2^o a$%.  If the result overflows then
 920  *              high-order bits are discarded, as usual.  We must have
 921  *              @0 < o < MPW_BITS@.
 922  */
 923
 924 void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
 925 {
 926   mpw b = a >> (MPW_BITS - o);
 927   a <<= o;
 928
 929   if (dv < dvl) {
 930     mpd x = (mpd)*dv - MPW(a);
 931     *dv++ = MPW(x);
 932     if (x >> MPW_BITS)
 933       b++;
 934     MPX_USUBN(dv, dvl, b);
 935   }
 936 }
 937
 938 /* --- @mpx_umul@ --- *
 939  *
 940  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 941  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 942  *              @const mpw *bv, *bvl@ = multiplier vector base and limit
 943  *
 944  * Returns:     ---
 945  *
 946  * Use:         Performs unsigned integer multiplication.  If the result
 947  *              overflows the desination vector, high-order bits are
 948  *              discarded.  The result vector may not overlap the argument
 949  *              vectors in any way.
 950  */
 951
 952 void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 953               const mpw *bv, const mpw *bvl)
 954 {
 955   /* --- This is probably worthwhile on a multiply --- */
 956
 957   MPX_SHRINK(av, avl);
 958   MPX_SHRINK(bv, bvl);
 959
 960   /* --- Deal with a multiply by zero --- */
 961
 962   if (bv == bvl) {
 963     MPX_ZERO(dv, dvl);
 964     return;
 965   }
 966
 967   /* --- Do the initial multiply and initialize the accumulator --- */
 968
 969   MPX_UMULN(dv, dvl, av, avl, *bv++);
 970
 971   /* --- Do the remaining multiply/accumulates --- */
 972
 973   while (dv < dvl && bv < bvl) {
 974     mpw m = *bv++;
 975     mpw c = 0;
 976     const mpw *avv = av;
 977     mpw *dvv = ++dv;
 978
 979     while (avv < avl) {
 980       mpd x;
 981       if (dvv >= dvl)
 982         goto next;
 983       x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
 984       *dvv++ = MPW(x);
 985       c = x >> MPW_BITS;
 986     }
 987     MPX_UADDN(dvv, dvl, c);
 988   next:;
 989   }
 990 }
 991
 992 /* --- @mpx_umuln@ --- *
 993  *
 994  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 995  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 996  *              @mpw m@ = multiplier
 997  *
 998  * Returns:     ---
 999  *
1000  * Use:         Multiplies a multiprecision integer by a single-word value.
1001  *              The destination and source may be equal.  The destination
1002  *              is completely cleared after use.
1003  */
1004
1005 void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
1006   { MPX_UMULN(dv, dvl, av, avl, m); }
1007
1008 /* --- @mpx_umlan@ --- *
1009  *
1010  * Arguments:   @mpw *dv, *dvl@ = destination/accumulator base and limit
1011  *              @const mpw *av, *avl@ = multiplicand vector base and limit
1012  *              @mpw m@ = multiplier
1013  *
1014  * Returns:     ---
1015  *
1016  * Use:         Multiplies a multiprecision integer by a single-word value
1017  *              and adds the result to an accumulator.
1018  */
1019
1020 void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
1021   { MPX_UMLAN(dv, dvl, av, avl, m); }
1022
1023 /* --- @mpx_usqr@ --- *
1024  *
1025  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
1026  *              @const mpw *av, *av@ = source vector base and limit
1027  *
1028  * Returns:     ---
1029  *
1030  * Use:         Performs unsigned integer squaring.  The result vector must
1031  *              not overlap the source vector in any way.
1032  */
1033
1034 void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
1035 {
1036   MPX_ZERO(dv, dvl);
1037
1038   /* --- Main loop --- */
1039
1040   while (av < avl) {
1041     const mpw *avv = av;
1042     mpw *dvv = dv;
1043     mpw a = *av;
1044     mpd c;
1045
1046     /* --- Stop if I've run out of destination --- */
1047
1048     if (dvv >= dvl)
1049       break;
1050
1051     /* --- Work out the square at this point in the proceedings --- */
1052
1053     {
1054       mpd x = (mpd)a * (mpd)a + *dvv;
1055       *dvv++ = MPW(x);
1056       c = MPW(x >> MPW_BITS);
1057     }
1058
1059     /* --- Now fix up the rest of the vector upwards --- */
1060
1061     avv++;
1062     while (dvv < dvl && avv < avl) {
1063       mpd x = (mpd)a * (mpd)*avv++;
1064       mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
1065       c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
1066       *dvv++ = MPW(y);
1067     }
1068     while (dvv < dvl && c) {
1069       mpd x = c + *dvv;
1070       *dvv++ = MPW(x);
1071       c = x >> MPW_BITS;
1072     }
1073
1074     /* --- Get ready for the next round --- */
1075
1076     av++;
1077     dv += 2;
1078   }
1079 }
1080
1081 /* --- @mpx_udiv@ --- *
1082  *
1083  * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
1084  *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
1085  *              @const mpw *dv, *dvl@ = divisor vector base and limit
1086  *              @mpw *sv, *svl@ = scratch workspace
1087  *
1088  * Returns:     ---
1089  *
1090  * Use:         Performs unsigned integer division.  If the result overflows
1091  *              the quotient vector, high-order bits are discarded.  (Clearly
1092  *              the remainder vector can't overflow.)  The various vectors
1093  *              may not overlap in any way.  Yes, I know it's a bit odd
1094  *              requiring the dividend to be in the result position but it
1095  *              does make some sense really.  The remainder must have
1096  *              headroom for at least two extra words.  The scratch space
1097  *              must be at least one word larger than the divisor.
1098  */
1099
1100 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
1101               const mpw *dv, const mpw *dvl,
1102               mpw *sv, mpw *svl)
1103 {
1104   unsigned norm = 0;
1105   size_t scale;
1106   mpw d, dd;
1107
1108   /* --- Initialize the quotient --- */
1109
1110   MPX_ZERO(qv, qvl);
1111
1112   /* --- Perform some sanity checks --- */
1113
1114   MPX_SHRINK(dv, dvl);
1115   assert(((void)"division by zero in mpx_udiv", dv < dvl));
1116
1117   /* --- Normalize the divisor --- *
1118    *
1119    * The algorithm requires that the divisor be at least two digits long.
1120    * This is easy to fix.
1121    */
1122
1123   {
1124     unsigned b;
1125
1126     d = dvl[-1];
1127     for (b = MPW_P2; b; b >>= 1) {
1128       if (d <= (MPW_MAX >> b)) {
1129         d <<= b;
1130         norm += b;
1131       }
1132     }
1133     if (dv + 1 == dvl)
1134       norm += MPW_BITS;
1135   }
1136
1137   /* --- Normalize the dividend/remainder to match --- */
1138
1139   if (norm) {
1140     mpx_lsl(rv, rvl, rv, rvl, norm);
1141     mpx_lsl(sv, svl, dv, dvl, norm);
1142     dv = sv;
1143     dvl = svl;
1144     MPX_SHRINK(dv, dvl);
1145   }
1146
1147   MPX_SHRINK(rv, rvl);
1148   d = dvl[-1];
1149   dd = dvl[-2];
1150
1151   /* --- Work out the relative scales --- */
1152
1153   {
1154     size_t rvn = rvl - rv;
1155     size_t dvn = dvl - dv;
1156
1157     /* --- If the divisor is clearly larger, notice this --- */
1158
1159     if (dvn > rvn) {
1160       mpx_lsr(rv, rvl, rv, rvl, norm);
1161       return;
1162     }
1163
1164     scale = rvn - dvn;
1165   }
1166
1167   /* --- Calculate the most significant quotient digit --- *
1168    *
1169    * Because the divisor has its top bit set, this can only happen once.  The
1170    * pointer arithmetic is a little contorted, to make sure that the
1171    * behaviour is defined.
1172    */
1173
1174   if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
1175     mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1176     if (qvl - qv > scale)
1177       qv[scale] = 1;
1178   }
1179
1180   /* --- Now for the main loop --- */
1181
1182   {
1183     mpw *rvv = rvl - 2;
1184
1185     while (scale) {
1186       mpw q;
1187       mpd rh;
1188
1189       /* --- Get an estimate for the next quotient digit --- */
1190
1191       mpw r = rvv[1];
1192       mpw rr = rvv[0];
1193       mpw rrr = *--rvv;
1194
1195       scale--;
1196       rh = ((mpd)r << MPW_BITS) | rr;
1197       if (r == d)
1198         q = MPW_MAX;
1199       else
1200         q = MPW(rh / d);
1201
1202       /* --- Refine the estimate --- */
1203
1204       {
1205         mpd yh = (mpd)d * q;
1206         mpd yy = (mpd)dd * q;
1207         mpw yl;
1208
1209         if (yy > MPW_MAX)
1210           yh += yy >> MPW_BITS;
1211         yl = MPW(yy);
1212
1213         while (yh > rh || (yh == rh && yl > rrr)) {
1214           q--;
1215           yh -= d;
1216           if (yl < dd)
1217             yh--;
1218           yl = MPW(yl - dd);
1219         }
1220       }
1221
1222       /* --- Remove a chunk from the dividend --- */
1223
1224       {
1225         mpw *svv;
1226         const mpw *dvv;
1227         mpw mc = 0, sc = 0;
1228
1229         /* --- Calculate the size of the chunk --- *
1230          *
1231          * This does the whole job of calculating @r >> scale - qd@.
1232          */
1233
1234         for (svv = rv + scale, dvv = dv;
1235              dvv < dvl && svv < rvl;
1236              svv++, dvv++) {
1237           mpd x = (mpd)*dvv * (mpd)q + mc;
1238           mc = x >> MPW_BITS;
1239           x = (mpd)*svv - MPW(x) - sc;
1240           *svv = MPW(x);
1241           if (x >> MPW_BITS)
1242             sc = 1;
1243           else
1244             sc = 0;
1245         }
1246
1247         if (svv < rvl) {
1248           mpd x = (mpd)*svv - mc - sc;
1249           *svv++ = MPW(x);
1250           if (x >> MPW_BITS)
1251             sc = MPW_MAX;
1252           else
1253             sc = 0;
1254           while (svv < rvl)
1255             *svv++ = sc;
1256         }
1257
1258         /* --- Fix if the quotient was too large --- *
1259          *
1260          * This doesn't seem to happen very often.
1261          */
1262
1263         if (rvl[-1] > MPW_MAX / 2) {
1264           mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1265           q--;
1266         }
1267       }
1268
1269       /* --- Done for another iteration --- */
1270
1271       if (qvl - qv > scale)
1272         qv[scale] = q;
1273       r = rr;
1274       rr = rrr;
1275     }
1276   }
1277
1278   /* --- Now fiddle with unnormalizing and things --- */
1279
1280   mpx_lsr(rv, rvl, rv, rvl, norm);
1281 }
1282
1283 /* --- @mpx_udivn@ --- *
1284  *
1285  * Arguments:   @mpw *qv, *qvl@ = storage for the quotient (may overlap
1286  *                      dividend)
1287  *              @const mpw *rv, *rvl@ = dividend
1288  *              @mpw d@ = single-precision divisor
1289  *
1290  * Returns:     Remainder after divison.
1291  *
1292  * Use:         Performs a single-precision division operation.
1293  */
1294
1295 mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
1296 {
1297   size_t i;
1298   size_t ql = qvl - qv;
1299   mpd r = 0;
1300
1301   i = rvl - rv;
1302   while (i > 0) {
1303     i--;
1304     r = (r << MPW_BITS) | rv[i];
1305     if (i < ql)
1306       qv[i] = r / d;
1307     r %= d;
1308   }
1309   return (MPW(r));
1310 }
1311
1312 /*----- Test rig ----------------------------------------------------------*/
1313
1314 #ifdef TEST_RIG
1315
1316 #include <mLib/alloc.h>
1317 #include <mLib/dstr.h>
1318 #include <mLib/quis.h>
1319 #include <mLib/testrig.h>
1320
1321 #include "mpscan.h"
1322
1323 #define ALLOC(v, vl, sz) do {                                           \
1324   size_t _sz = (sz);                                                    \
1325   mpw *_vv = xmalloc(MPWS(_sz));                                        \
1326   mpw *_vvl = _vv + _sz;                                                \
1327   (v) = _vv;                                                            \
1328   (vl) = _vvl;                                                          \
1329 } while (0)
1330
1331 #define LOAD(v, vl, d) do {                                             \
1332   const dstr *_d = (d);                                                 \
1333   mpw *_v, *_vl;                                                        \
1334   ALLOC(_v, _vl, MPW_RQ(_d->len));                                      \
1335   mpx_loadb(_v, _vl, _d->buf, _d->len);                                 \
1336   (v) = _v;                                                             \
1337   (vl) = _vl;                                                           \
1338 } while (0)
1339
1340 #define MAX(x, y) ((x) > (y) ? (x) : (y))
1341
1342 static void dumpbits(const char *msg, const void *pp, size_t sz)
1343 {
1344   const octet *p = pp;
1345   fputs(msg, stderr);
1346   for (; sz; sz--)
1347     fprintf(stderr, " %02x", *p++);
1348   fputc('\n', stderr);
1349 }
1350
1351 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
1352 {
1353   fputs(msg, stderr);
1354   MPX_SHRINK(v, vl);
1355   while (v < vl)
1356     fprintf(stderr, " %08lx", (unsigned long)*--vl);
1357   fputc('\n', stderr);
1358 }
1359
1360 static int chkscan(const mpw *v, const mpw *vl,
1361                    const void *pp, size_t sz, int step)
1362 {
1363   mpscan mps;
1364   const octet *p = pp;
1365   unsigned bit = 0;
1366   int ok = 1;
1367
1368   mpscan_initx(&mps, v, vl);
1369   while (sz) {
1370     unsigned x = *p;
1371     int i;
1372     p += step;
1373     for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
1374       if (MPSCAN_BIT(&mps) != (x & 1)) {
1375         fprintf(stderr,
1376                 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1377                 step, bit, x & 1, MPSCAN_BIT(&mps));
1378         ok = 0;
1379       }
1380       x >>= 1;
1381       bit++;
1382     }
1383     sz--;
1384   }
1385
1386   return (ok);
1387 }
1388
1389 static int loadstore(dstr *v)
1390 {
1391   dstr d = DSTR_INIT;
1392   size_t sz = MPW_RQ(v->len) * 2, diff;
1393   mpw *m, *ml;
1394   int ok = 1;
1395
1396   dstr_ensure(&d, v->len);
1397   m = xmalloc(MPWS(sz));
1398
1399   for (diff = 0; diff < sz; diff += 5) {
1400     size_t oct;
1401
1402     ml = m + sz - diff;
1403
1404     mpx_loadl(m, ml, v->buf, v->len);
1405     if (!chkscan(m, ml, v->buf, v->len, +1))
1406       ok = 0;
1407     MPX_OCTETS(oct, m, ml);
1408     mpx_storel(m, ml, d.buf, d.sz);
1409     if (memcmp(d.buf, v->buf, oct) != 0) {
1410       dumpbits("\n*** storel failed", d.buf, d.sz);
1411       ok = 0;
1412     }
1413
1414     mpx_loadb(m, ml, v->buf, v->len);
1415     if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1416       ok = 0;
1417     MPX_OCTETS(oct, m, ml);
1418     mpx_storeb(m, ml, d.buf, d.sz);
1419     if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1420       dumpbits("\n*** storeb failed", d.buf, d.sz);
1421       ok = 0;
1422     }
1423   }
1424
1425   if (!ok)
1426     dumpbits("input data", v->buf, v->len);
1427
1428   xfree(m);
1429   dstr_destroy(&d);
1430   return (ok);
1431 }
1432
1433 static int twocl(dstr *v)
1434 {
1435   dstr d = DSTR_INIT;
1436   mpw *m, *ml;
1437   size_t sz;
1438   int ok = 1;
1439
1440   sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1441   dstr_ensure(&d, sz);
1442
1443   sz = MPW_RQ(sz);
1444   m = xmalloc(MPWS(sz));
1445   ml = m + sz;
1446
1447   mpx_loadl(m, ml, v[0].buf, v[0].len);
1448   mpx_storel2cn(m, ml, d.buf, v[1].len);
1449   if (memcmp(d.buf, v[1].buf, v[1].len)) {
1450     dumpbits("\n*** storel2cn failed", d.buf, v[1].len);
1451     ok = 0;
1452   }
1453
1454   mpx_loadl2cn(m, ml, v[1].buf, v[1].len);
1455   mpx_storel(m, ml, d.buf, v[0].len);
1456   if (memcmp(d.buf, v[0].buf, v[0].len)) {
1457     dumpbits("\n*** loadl2cn failed", d.buf, v[0].len);
1458     ok = 0;
1459   }
1460
1461   if (!ok) {
1462     dumpbits("pos", v[0].buf, v[0].len);
1463     dumpbits("neg", v[1].buf, v[1].len);
1464   }
1465
1466   xfree(m);
1467   dstr_destroy(&d);
1468
1469   return (ok);
1470 }
1471
1472 static int twocb(dstr *v)
1473 {
1474   dstr d = DSTR_INIT;
1475   mpw *m, *ml;
1476   size_t sz;
1477   int ok = 1;
1478
1479   sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1480   dstr_ensure(&d, sz);
1481
1482   sz = MPW_RQ(sz);
1483   m = xmalloc(MPWS(sz));
1484   ml = m + sz;
1485
1486   mpx_loadb(m, ml, v[0].buf, v[0].len);
1487   mpx_storeb2cn(m, ml, d.buf, v[1].len);
1488   if (memcmp(d.buf, v[1].buf, v[1].len)) {
1489     dumpbits("\n*** storeb2cn failed", d.buf, v[1].len);
1490     ok = 0;
1491   }
1492
1493   mpx_loadb2cn(m, ml, v[1].buf, v[1].len);
1494   mpx_storeb(m, ml, d.buf, v[0].len);
1495   if (memcmp(d.buf, v[0].buf, v[0].len)) {
1496     dumpbits("\n*** loadb2cn failed", d.buf, v[0].len);
1497     ok = 0;
1498   }
1499
1500   if (!ok) {
1501     dumpbits("pos", v[0].buf, v[0].len);
1502     dumpbits("neg", v[1].buf, v[1].len);
1503   }
1504
1505   xfree(m);
1506   dstr_destroy(&d);
1507
1508   return (ok);
1509 }
1510
1511 static int lsl(dstr *v)
1512 {
1513   mpw *a, *al;
1514   int n = *(int *)v[1].buf;
1515   mpw *c, *cl;
1516   mpw *d, *dl;
1517   int ok = 1;
1518
1519   LOAD(a, al, &v[0]);
1520   LOAD(c, cl, &v[2]);
1521   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1522
1523   mpx_lsl(d, dl, a, al, n);
1524   if (!mpx_ueq(d, dl, c, cl)) {
1525     fprintf(stderr, "\n*** lsl(%i) failed\n", n);
1526     dumpmp("       a", a, al);
1527     dumpmp("expected", c, cl);
1528     dumpmp("  result", d, dl);
1529     ok = 0;
1530   }
1531
1532   xfree(a); xfree(c); xfree(d);
1533   return (ok);
1534 }
1535
1536 static int lslc(dstr *v)
1537 {
1538   mpw *a, *al;
1539   int n = *(int *)v[1].buf;
1540   mpw *c, *cl;
1541   mpw *d, *dl;
1542   int ok = 1;
1543
1544   LOAD(a, al, &v[0]);
1545   LOAD(c, cl, &v[2]);
1546   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1547
1548   mpx_lslc(d, dl, a, al, n);
1549   if (!mpx_ueq(d, dl, c, cl)) {
1550     fprintf(stderr, "\n*** lslc(%i) failed\n", n);
1551     dumpmp("       a", a, al);
1552     dumpmp("expected", c, cl);
1553     dumpmp("  result", d, dl);
1554     ok = 0;
1555   }
1556
1557   xfree(a); xfree(c); xfree(d);
1558   return (ok);
1559 }
1560
1561 static int lsr(dstr *v)
1562 {
1563   mpw *a, *al;
1564   int n = *(int *)v[1].buf;
1565   mpw *c, *cl;
1566   mpw *d, *dl;
1567   int ok = 1;
1568
1569   LOAD(a, al, &v[0]);
1570   LOAD(c, cl, &v[2]);
1571   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1572
1573   mpx_lsr(d, dl, a, al, n);
1574   if (!mpx_ueq(d, dl, c, cl)) {
1575     fprintf(stderr, "\n*** lsr(%i) failed\n", n);
1576     dumpmp("       a", a, al);
1577     dumpmp("expected", c, cl);
1578     dumpmp("  result", d, dl);
1579     ok = 0;
1580   }
1581
1582   xfree(a); xfree(c); xfree(d);
1583   return (ok);
1584 }
1585
1586 static int uadd(dstr *v)
1587 {
1588   mpw *a, *al;
1589   mpw *b, *bl;
1590   mpw *c, *cl;
1591   mpw *d, *dl;
1592   int ok = 1;
1593
1594   LOAD(a, al, &v[0]);
1595   LOAD(b, bl, &v[1]);
1596   LOAD(c, cl, &v[2]);
1597   ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1598
1599   mpx_uadd(d, dl, a, al, b, bl);
1600   if (!mpx_ueq(d, dl, c, cl)) {
1601     fprintf(stderr, "\n*** uadd failed\n");
1602     dumpmp("       a", a, al);
1603     dumpmp("       b", b, bl);
1604     dumpmp("expected", c, cl);
1605     dumpmp("  result", d, dl);
1606     ok = 0;
1607   }
1608
1609   xfree(a); xfree(b); xfree(c); xfree(d);
1610   return (ok);
1611 }
1612
1613 static int usub(dstr *v)
1614 {
1615   mpw *a, *al;
1616   mpw *b, *bl;
1617   mpw *c, *cl;
1618   mpw *d, *dl;
1619   int ok = 1;
1620
1621   LOAD(a, al, &v[0]);
1622   LOAD(b, bl, &v[1]);
1623   LOAD(c, cl, &v[2]);
1624   ALLOC(d, dl, al - a);
1625
1626   mpx_usub(d, dl, a, al, b, bl);
1627   if (!mpx_ueq(d, dl, c, cl)) {
1628     fprintf(stderr, "\n*** usub failed\n");
1629     dumpmp("       a", a, al);
1630     dumpmp("       b", b, bl);
1631     dumpmp("expected", c, cl);
1632     dumpmp("  result", d, dl);
1633     ok = 0;
1634   }
1635
1636   xfree(a); xfree(b); xfree(c); xfree(d);
1637   return (ok);
1638 }
1639
1640 static int umul(dstr *v)
1641 {
1642   mpw *a, *al;
1643   mpw *b, *bl;
1644   mpw *c, *cl;
1645   mpw *d, *dl;
1646   int ok = 1;
1647
1648   LOAD(a, al, &v[0]);
1649   LOAD(b, bl, &v[1]);
1650   LOAD(c, cl, &v[2]);
1651   ALLOC(d, dl, (al - a) + (bl - b));
1652
1653   mpx_umul(d, dl, a, al, b, bl);
1654   if (!mpx_ueq(d, dl, c, cl)) {
1655     fprintf(stderr, "\n*** umul failed\n");
1656     dumpmp("       a", a, al);
1657     dumpmp("       b", b, bl);
1658     dumpmp("expected", c, cl);
1659     dumpmp("  result", d, dl);
1660     ok = 0;
1661   }
1662
1663   xfree(a); xfree(b); xfree(c); xfree(d);
1664   return (ok);
1665 }
1666
1667 static int usqr(dstr *v)
1668 {
1669   mpw *a, *al;
1670   mpw *c, *cl;
1671   mpw *d, *dl;
1672   int ok = 1;
1673
1674   LOAD(a, al, &v[0]);
1675   LOAD(c, cl, &v[1]);
1676   ALLOC(d, dl, 2 * (al - a));
1677
1678   mpx_usqr(d, dl, a, al);
1679   if (!mpx_ueq(d, dl, c, cl)) {
1680     fprintf(stderr, "\n*** usqr failed\n");
1681     dumpmp("       a", a, al);
1682     dumpmp("expected", c, cl);
1683     dumpmp("  result", d, dl);
1684     ok = 0;
1685   }
1686
1687   xfree(a); xfree(c); xfree(d);
1688   return (ok);
1689 }
1690
1691 static int udiv(dstr *v)
1692 {
1693   mpw *a, *al;
1694   mpw *b, *bl;
1695   mpw *q, *ql;
1696   mpw *r, *rl;
1697   mpw *qq, *qql;
1698   mpw *s, *sl;
1699   int ok = 1;
1700
1701   ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1702   LOAD(b, bl, &v[1]);
1703   LOAD(q, ql, &v[2]);
1704   LOAD(r, rl, &v[3]);
1705   ALLOC(qq, qql, al - a);
1706   ALLOC(s, sl, (bl - b) + 1);
1707
1708   mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1709   if (!mpx_ueq(qq, qql, q, ql) ||
1710       !mpx_ueq(a, al, r, rl)) {
1711     fprintf(stderr, "\n*** udiv failed\n");
1712     dumpmp(" divisor", b, bl);
1713     dumpmp("expect r", r, rl);
1714     dumpmp("result r", a, al);
1715     dumpmp("expect q", q, ql);
1716     dumpmp("result q", qq, qql);
1717     ok = 0;
1718   }
1719
1720   xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq);
1721   return (ok);
1722 }
1723
1724 static test_chunk defs[] = {
1725   { "load-store", loadstore, { &type_hex, 0 } },
1726   { "2cl", twocl, { &type_hex, &type_hex, } },
1727   { "2cb", twocb, { &type_hex, &type_hex, } },
1728   { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
1729   { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
1730   { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1731   { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1732   { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1733   { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1734   { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1735   { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1736   { 0, 0, { 0 } }
1737 };
1738
1739 int main(int argc, char *argv[])
1740 {
1741   test_run(argc, argv, defs, SRCDIR"/t/mpx");
1742   return (0);
1743 }
1744
1745 #endif
1746
1747 /*----- That's all, folks -------------------------------------------------*/