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 #define MPX_BITBINOP(string)                                            \
 665                                                                         \
 666 void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,  \
 667                      const mpw *bv, const mpw *bvl)                     \
 668 {                                                                       \
 669   MPX_SHRINK(av, avl);                                                  \
 670   MPX_SHRINK(bv, bvl);                                                  \
 671                                                                         \
 672   while (dv < dvl) {                                                    \
 673     mpw a, b;                                                           \
 674     a = (av < avl) ? *av++ : 0;                                         \
 675     b = (bv < bvl) ? *bv++ : 0;                                         \
 676     *dv++ = B##string(a, b);                                            \
 677   }                                                                     \
 678 }
 679
 680 MPX_DOBIN(MPX_BITBINOP)
 681
 682 void mpx_not(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
 683 {
 684   MPX_SHRINK(av, avl);
 685
 686   while (dv < dvl) {
 687     mpw a;
 688     a = (av < avl) ? *av++ : 0;
 689     *dv++ = ~a;
 690   }
 691 }
 692
 693 /*----- Unsigned arithmetic -----------------------------------------------*/
 694
 695 /* --- @mpx_2c@ --- *
 696  *
 697  * Arguments:   @mpw *dv, *dvl@ = destination vector
 698  *              @const mpw *v, *vl@ = source vector
 699  *
 700  * Returns:     ---
 701  *
 702  * Use:         Calculates the two's complement of @v@.
 703  */
 704
 705 void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
 706 {
 707   mpw c = 0;
 708   while (dv < dvl && v < vl)
 709     *dv++ = c = MPW(~*v++);
 710   if (dv < dvl) {
 711     if (c > MPW_MAX / 2)
 712       c = MPW(~0);
 713     while (dv < dvl)
 714       *dv++ = c;
 715   }
 716   MPX_UADDN(dv, dvl, 1);
 717 }
 718
 719 /* --- @mpx_ueq@ --- *
 720  *
 721  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 722  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 723  *
 724  * Returns:     Nonzero if the two vectors are equal.
 725  *
 726  * Use:         Performs an unsigned integer test for equality.
 727  */
 728
 729 int mpx_ueq(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
 730 {
 731   MPX_SHRINK(av, avl);
 732   MPX_SHRINK(bv, bvl);
 733   if (avl - av != bvl - bv)
 734     return (0);
 735   while (av < avl) {
 736     if (*av++ != *bv++)
 737       return (0);
 738   }
 739   return (1);
 740 }
 741
 742 /* --- @mpx_ucmp@ --- *
 743  *
 744  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 745  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 746  *
 747  * Returns:     Less than, equal to, or greater than zero depending on
 748  *              whether @a@ is less than, equal to or greater than @b@,
 749  *              respectively.
 750  *
 751  * Use:         Performs an unsigned integer comparison.
 752  */
 753
 754 int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
 755 {
 756   MPX_SHRINK(av, avl);
 757   MPX_SHRINK(bv, bvl);
 758
 759   if (avl - av > bvl - bv)
 760     return (+1);
 761   else if (avl - av < bvl - bv)
 762     return (-1);
 763   else while (avl > av) {
 764     mpw a = *--avl, b = *--bvl;
 765     if (a > b)
 766       return (+1);
 767     else if (a < b)
 768       return (-1);
 769   }
 770   return (0);
 771 }
 772
 773 /* --- @mpx_uadd@ --- *
 774  *
 775  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 776  *              @const mpw *av, *avl@ = first addend vector base and limit
 777  *              @const mpw *bv, *bvl@ = second addend vector base and limit
 778  *
 779  * Returns:     ---
 780  *
 781  * Use:         Performs unsigned integer addition.  If the result overflows
 782  *              the destination vector, high-order bits are discarded.  This
 783  *              means that two's complement addition happens more or less for
 784  *              free, although that's more a side-effect than anything else.
 785  *              The result vector may be equal to either or both source
 786  *              vectors, but may not otherwise overlap them.
 787  */
 788
 789 void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 790               const mpw *bv, const mpw *bvl)
 791 {
 792   mpw c = 0;
 793
 794   while (av < avl || bv < bvl) {
 795     mpw a, b;
 796     mpd x;
 797     if (dv >= dvl)
 798       return;
 799     a = (av < avl) ? *av++ : 0;
 800     b = (bv < bvl) ? *bv++ : 0;
 801     x = (mpd)a + (mpd)b + c;
 802     *dv++ = MPW(x);
 803     c = x >> MPW_BITS;
 804   }
 805   if (dv < dvl) {
 806     *dv++ = c;
 807     MPX_ZERO(dv, dvl);
 808   }
 809 }
 810
 811 /* --- @mpx_uaddn@ --- *
 812  *
 813  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 814  *              @mpw n@ = other addend
 815  *
 816  * Returns:     ---
 817  *
 818  * Use:         Adds a small integer to a multiprecision number.
 819  */
 820
 821 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
 822
 823 /* --- @mpx_uaddnlsl@ --- *
 824  *
 825  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 826  *              @mpw a@ = second argument
 827  *              @unsigned o@ = offset in bits
 828  *
 829  * Returns:     ---
 830  *
 831  * Use:         Computes %$d + 2^o a$%.  If the result overflows then
 832  *              high-order bits are discarded, as usual.  We must have
 833  *              @0 < o < MPW_BITS@.
 834  */
 835
 836 void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
 837 {
 838   mpd x = (mpd)a << o;
 839
 840   while (x && dv < dvl) {
 841     x += *dv;
 842     *dv++ = MPW(x);
 843     x >>= MPW_BITS;
 844   }
 845 }
 846
 847 /* --- @mpx_usub@ --- *
 848  *
 849  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 850  *              @const mpw *av, *avl@ = first argument vector base and limit
 851  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 852  *
 853  * Returns:     ---
 854  *
 855  * Use:         Performs unsigned integer subtraction.  If the result
 856  *              overflows the destination vector, high-order bits are
 857  *              discarded.  This means that two's complement subtraction
 858  *              happens more or less for free, althuogh that's more a side-
 859  *              effect than anything else.  The result vector may be equal to
 860  *              either or both source vectors, but may not otherwise overlap
 861  *              them.
 862  */
 863
 864 void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 865               const mpw *bv, const mpw *bvl)
 866 {
 867   mpw c = 0;
 868
 869   while (av < avl || bv < bvl) {
 870     mpw a, b;
 871     mpd x;
 872     if (dv >= dvl)
 873       return;
 874     a = (av < avl) ? *av++ : 0;
 875     b = (bv < bvl) ? *bv++ : 0;
 876     x = (mpd)a - (mpd)b - c;
 877     *dv++ = MPW(x);
 878     if (x >> MPW_BITS)
 879       c = 1;
 880     else
 881       c = 0;
 882   }
 883   if (c)
 884     c = MPW_MAX;
 885   while (dv < dvl)
 886     *dv++ = c;
 887 }
 888
 889 /* --- @mpx_usubn@ --- *
 890  *
 891  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 892  *              @n@ = subtrahend
 893  *
 894  * Returns:     ---
 895  *
 896  * Use:         Subtracts a small integer from a multiprecision number.
 897  */
 898
 899 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
 900
 901 /* --- @mpx_uaddnlsl@ --- *
 902  *
 903  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 904  *              @mpw a@ = second argument
 905  *              @unsigned o@ = offset in bits
 906  *
 907  * Returns:     ---
 908  *
 909  * Use:         Computes %$d + 2^o a$%.  If the result overflows then
 910  *              high-order bits are discarded, as usual.  We must have
 911  *              @0 < o < MPW_BITS@.
 912  */
 913
 914 void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
 915 {
 916   mpw b = a >> (MPW_BITS - o);
 917   a <<= o;
 918
 919   if (dv < dvl) {
 920     mpd x = (mpd)*dv - MPW(a);
 921     *dv++ = MPW(x);
 922     if (x >> MPW_BITS)
 923       b++;
 924     MPX_USUBN(dv, dvl, b);
 925   }
 926 }
 927
 928 /* --- @mpx_umul@ --- *
 929  *
 930  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 931  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 932  *              @const mpw *bv, *bvl@ = multiplier vector base and limit
 933  *
 934  * Returns:     ---
 935  *
 936  * Use:         Performs unsigned integer multiplication.  If the result
 937  *              overflows the desination vector, high-order bits are
 938  *              discarded.  The result vector may not overlap the argument
 939  *              vectors in any way.
 940  */
 941
 942 void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 943               const mpw *bv, const mpw *bvl)
 944 {
 945   /* --- This is probably worthwhile on a multiply --- */
 946
 947   MPX_SHRINK(av, avl);
 948   MPX_SHRINK(bv, bvl);
 949
 950   /* --- Deal with a multiply by zero --- */
 951
 952   if (bv == bvl) {
 953     MPX_ZERO(dv, dvl);
 954     return;
 955   }
 956
 957   /* --- Do the initial multiply and initialize the accumulator --- */
 958
 959   MPX_UMULN(dv, dvl, av, avl, *bv++);
 960
 961   /* --- Do the remaining multiply/accumulates --- */
 962
 963   while (dv < dvl && bv < bvl) {
 964     mpw m = *bv++;
 965     mpw c = 0;
 966     const mpw *avv = av;
 967     mpw *dvv = ++dv;
 968
 969     while (avv < avl) {
 970       mpd x;
 971       if (dvv >= dvl)
 972         goto next;
 973       x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
 974       *dvv++ = MPW(x);
 975       c = x >> MPW_BITS;
 976     }
 977     MPX_UADDN(dvv, dvl, c);
 978   next:;
 979   }
 980 }
 981
 982 /* --- @mpx_umuln@ --- *
 983  *
 984  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 985  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 986  *              @mpw m@ = multiplier
 987  *
 988  * Returns:     ---
 989  *
 990  * Use:         Multiplies a multiprecision integer by a single-word value.
 991  *              The destination and source may be equal.  The destination
 992  *              is completely cleared after use.
 993  */
 994
 995 void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
 996   { MPX_UMULN(dv, dvl, av, avl, m); }
 997
 998 /* --- @mpx_umlan@ --- *
 999  *
1000  * Arguments:   @mpw *dv, *dvl@ = destination/accumulator base and limit
1001  *              @const mpw *av, *avl@ = multiplicand vector base and limit
1002  *              @mpw m@ = multiplier
1003  *
1004  * Returns:     ---
1005  *
1006  * Use:         Multiplies a multiprecision integer by a single-word value
1007  *              and adds the result to an accumulator.
1008  */
1009
1010 void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
1011   { MPX_UMLAN(dv, dvl, av, avl, m); }
1012
1013 /* --- @mpx_usqr@ --- *
1014  *
1015  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
1016  *              @const mpw *av, *av@ = source vector base and limit
1017  *
1018  * Returns:     ---
1019  *
1020  * Use:         Performs unsigned integer squaring.  The result vector must
1021  *              not overlap the source vector in any way.
1022  */
1023
1024 void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
1025 {
1026   MPX_ZERO(dv, dvl);
1027
1028   /* --- Main loop --- */
1029
1030   while (av < avl) {
1031     const mpw *avv = av;
1032     mpw *dvv = dv;
1033     mpw a = *av;
1034     mpd c;
1035
1036     /* --- Stop if I've run out of destination --- */
1037
1038     if (dvv >= dvl)
1039       break;
1040
1041     /* --- Work out the square at this point in the proceedings --- */
1042
1043     {
1044       mpd x = (mpd)a * (mpd)a + *dvv;
1045       *dvv++ = MPW(x);
1046       c = MPW(x >> MPW_BITS);
1047     }
1048
1049     /* --- Now fix up the rest of the vector upwards --- */
1050
1051     avv++;
1052     while (dvv < dvl && avv < avl) {
1053       mpd x = (mpd)a * (mpd)*avv++;
1054       mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
1055       c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
1056       *dvv++ = MPW(y);
1057     }
1058     while (dvv < dvl && c) {
1059       mpd x = c + *dvv;
1060       *dvv++ = MPW(x);
1061       c = x >> MPW_BITS;
1062     }
1063
1064     /* --- Get ready for the next round --- */
1065
1066     av++;
1067     dv += 2;
1068   }
1069 }
1070
1071 /* --- @mpx_udiv@ --- *
1072  *
1073  * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
1074  *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
1075  *              @const mpw *dv, *dvl@ = divisor vector base and limit
1076  *              @mpw *sv, *svl@ = scratch workspace
1077  *
1078  * Returns:     ---
1079  *
1080  * Use:         Performs unsigned integer division.  If the result overflows
1081  *              the quotient vector, high-order bits are discarded.  (Clearly
1082  *              the remainder vector can't overflow.)  The various vectors
1083  *              may not overlap in any way.  Yes, I know it's a bit odd
1084  *              requiring the dividend to be in the result position but it
1085  *              does make some sense really.  The remainder must have
1086  *              headroom for at least two extra words.  The scratch space
1087  *              must be at least one word larger than the divisor.
1088  */
1089
1090 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
1091               const mpw *dv, const mpw *dvl,
1092               mpw *sv, mpw *svl)
1093 {
1094   unsigned norm = 0;
1095   size_t scale;
1096   mpw d, dd;
1097
1098   /* --- Initialize the quotient --- */
1099
1100   MPX_ZERO(qv, qvl);
1101
1102   /* --- Perform some sanity checks --- */
1103
1104   MPX_SHRINK(dv, dvl);
1105   assert(((void)"division by zero in mpx_udiv", dv < dvl));
1106
1107   /* --- Normalize the divisor --- *
1108    *
1109    * The algorithm requires that the divisor be at least two digits long.
1110    * This is easy to fix.
1111    */
1112
1113   {
1114     unsigned b;
1115
1116     d = dvl[-1];
1117     for (b = MPW_P2; b; b >>= 1) {
1118       if (d <= (MPW_MAX >> b)) {
1119         d <<= b;
1120         norm += b;
1121       }
1122     }
1123     if (dv + 1 == dvl)
1124       norm += MPW_BITS;
1125   }
1126
1127   /* --- Normalize the dividend/remainder to match --- */
1128
1129   if (norm) {
1130     mpx_lsl(rv, rvl, rv, rvl, norm);
1131     mpx_lsl(sv, svl, dv, dvl, norm);
1132     dv = sv;
1133     dvl = svl;
1134     MPX_SHRINK(dv, dvl);
1135   }
1136
1137   MPX_SHRINK(rv, rvl);
1138   d = dvl[-1];
1139   dd = dvl[-2];
1140
1141   /* --- Work out the relative scales --- */
1142
1143   {
1144     size_t rvn = rvl - rv;
1145     size_t dvn = dvl - dv;
1146
1147     /* --- If the divisor is clearly larger, notice this --- */
1148
1149     if (dvn > rvn) {
1150       mpx_lsr(rv, rvl, rv, rvl, norm);
1151       return;
1152     }
1153
1154     scale = rvn - dvn;
1155   }
1156
1157   /* --- Calculate the most significant quotient digit --- *
1158    *
1159    * Because the divisor has its top bit set, this can only happen once.  The
1160    * pointer arithmetic is a little contorted, to make sure that the
1161    * behaviour is defined.
1162    */
1163
1164   if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
1165     mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1166     if (qvl - qv > scale)
1167       qv[scale] = 1;
1168   }
1169
1170   /* --- Now for the main loop --- */
1171
1172   {
1173     mpw *rvv = rvl - 2;
1174
1175     while (scale) {
1176       mpw q;
1177       mpd rh;
1178
1179       /* --- Get an estimate for the next quotient digit --- */
1180
1181       mpw r = rvv[1];
1182       mpw rr = rvv[0];
1183       mpw rrr = *--rvv;
1184
1185       scale--;
1186       rh = ((mpd)r << MPW_BITS) | rr;
1187       if (r == d)
1188         q = MPW_MAX;
1189       else
1190         q = MPW(rh / d);
1191
1192       /* --- Refine the estimate --- */
1193
1194       {
1195         mpd yh = (mpd)d * q;
1196         mpd yy = (mpd)dd * q;
1197         mpw yl;
1198
1199         if (yy > MPW_MAX)
1200           yh += yy >> MPW_BITS;
1201         yl = MPW(yy);
1202
1203         while (yh > rh || (yh == rh && yl > rrr)) {
1204           q--;
1205           yh -= d;
1206           if (yl < dd)
1207             yh--;
1208           yl = MPW(yl - dd);
1209         }
1210       }
1211
1212       /* --- Remove a chunk from the dividend --- */
1213
1214       {
1215         mpw *svv;
1216         const mpw *dvv;
1217         mpw mc = 0, sc = 0;
1218
1219         /* --- Calculate the size of the chunk --- *
1220          *
1221          * This does the whole job of calculating @r >> scale - qd@.
1222          */
1223
1224         for (svv = rv + scale, dvv = dv;
1225              dvv < dvl && svv < rvl;
1226              svv++, dvv++) {
1227           mpd x = (mpd)*dvv * (mpd)q + mc;
1228           mc = x >> MPW_BITS;
1229           x = (mpd)*svv - MPW(x) - sc;
1230           *svv = MPW(x);
1231           if (x >> MPW_BITS)
1232             sc = 1;
1233           else
1234             sc = 0;
1235         }
1236
1237         if (svv < rvl) {
1238           mpd x = (mpd)*svv - mc - sc;
1239           *svv++ = MPW(x);
1240           if (x >> MPW_BITS)
1241             sc = MPW_MAX;
1242           else
1243             sc = 0;
1244           while (svv < rvl)
1245             *svv++ = sc;
1246         }
1247
1248         /* --- Fix if the quotient was too large --- *
1249          *
1250          * This doesn't seem to happen very often.
1251          */
1252
1253         if (rvl[-1] > MPW_MAX / 2) {
1254           mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1255           q--;
1256         }
1257       }
1258
1259       /* --- Done for another iteration --- */
1260
1261       if (qvl - qv > scale)
1262         qv[scale] = q;
1263       r = rr;
1264       rr = rrr;
1265     }
1266   }
1267
1268   /* --- Now fiddle with unnormalizing and things --- */
1269
1270   mpx_lsr(rv, rvl, rv, rvl, norm);
1271 }
1272
1273 /* --- @mpx_udivn@ --- *
1274  *
1275  * Arguments:   @mpw *qv, *qvl@ = storage for the quotient (may overlap
1276  *                      dividend)
1277  *              @const mpw *rv, *rvl@ = dividend
1278  *              @mpw d@ = single-precision divisor
1279  *
1280  * Returns:     Remainder after divison.
1281  *
1282  * Use:         Performs a single-precision division operation.
1283  */
1284
1285 mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
1286 {
1287   size_t i;
1288   size_t ql = qvl - qv;
1289   mpd r = 0;
1290
1291   i = rvl - rv;
1292   while (i > 0) {
1293     i--;
1294     r = (r << MPW_BITS) | rv[i];
1295     if (i < ql)
1296       qv[i] = r / d;
1297     r %= d;
1298   }
1299   return (MPW(r));
1300 }
1301
1302 /*----- Test rig ----------------------------------------------------------*/
1303
1304 #ifdef TEST_RIG
1305
1306 #include <mLib/alloc.h>
1307 #include <mLib/dstr.h>
1308 #include <mLib/quis.h>
1309 #include <mLib/testrig.h>
1310
1311 #include "mpscan.h"
1312
1313 #define ALLOC(v, vl, sz) do {                                           \
1314   size_t _sz = (sz);                                                    \
1315   mpw *_vv = xmalloc(MPWS(_sz));                                        \
1316   mpw *_vvl = _vv + _sz;                                                \
1317   (v) = _vv;                                                            \
1318   (vl) = _vvl;                                                          \
1319 } while (0)
1320
1321 #define LOAD(v, vl, d) do {                                             \
1322   const dstr *_d = (d);                                                 \
1323   mpw *_v, *_vl;                                                        \
1324   ALLOC(_v, _vl, MPW_RQ(_d->len));                                      \
1325   mpx_loadb(_v, _vl, _d->buf, _d->len);                                 \
1326   (v) = _v;                                                             \
1327   (vl) = _vl;                                                           \
1328 } while (0)
1329
1330 #define MAX(x, y) ((x) > (y) ? (x) : (y))
1331
1332 static void dumpbits(const char *msg, const void *pp, size_t sz)
1333 {
1334   const octet *p = pp;
1335   fputs(msg, stderr);
1336   for (; sz; sz--)
1337     fprintf(stderr, " %02x", *p++);
1338   fputc('\n', stderr);
1339 }
1340
1341 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
1342 {
1343   fputs(msg, stderr);
1344   MPX_SHRINK(v, vl);
1345   while (v < vl)
1346     fprintf(stderr, " %08lx", (unsigned long)*--vl);
1347   fputc('\n', stderr);
1348 }
1349
1350 static int chkscan(const mpw *v, const mpw *vl,
1351                    const void *pp, size_t sz, int step)
1352 {
1353   mpscan mps;
1354   const octet *p = pp;
1355   unsigned bit = 0;
1356   int ok = 1;
1357
1358   mpscan_initx(&mps, v, vl);
1359   while (sz) {
1360     unsigned x = *p;
1361     int i;
1362     p += step;
1363     for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
1364       if (MPSCAN_BIT(&mps) != (x & 1)) {
1365         fprintf(stderr,
1366                 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1367                 step, bit, x & 1, MPSCAN_BIT(&mps));
1368         ok = 0;
1369       }
1370       x >>= 1;
1371       bit++;
1372     }
1373     sz--;
1374   }
1375
1376   return (ok);
1377 }
1378
1379 static int loadstore(dstr *v)
1380 {
1381   dstr d = DSTR_INIT;
1382   size_t sz = MPW_RQ(v->len) * 2, diff;
1383   mpw *m, *ml;
1384   int ok = 1;
1385
1386   dstr_ensure(&d, v->len);
1387   m = xmalloc(MPWS(sz));
1388
1389   for (diff = 0; diff < sz; diff += 5) {
1390     size_t oct;
1391
1392     ml = m + sz - diff;
1393
1394     mpx_loadl(m, ml, v->buf, v->len);
1395     if (!chkscan(m, ml, v->buf, v->len, +1))
1396       ok = 0;
1397     MPX_OCTETS(oct, m, ml);
1398     mpx_storel(m, ml, d.buf, d.sz);
1399     if (memcmp(d.buf, v->buf, oct) != 0) {
1400       dumpbits("\n*** storel failed", d.buf, d.sz);
1401       ok = 0;
1402     }
1403
1404     mpx_loadb(m, ml, v->buf, v->len);
1405     if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1406       ok = 0;
1407     MPX_OCTETS(oct, m, ml);
1408     mpx_storeb(m, ml, d.buf, d.sz);
1409     if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1410       dumpbits("\n*** storeb failed", d.buf, d.sz);
1411       ok = 0;
1412     }
1413   }
1414
1415   if (!ok)
1416     dumpbits("input data", v->buf, v->len);
1417
1418   xfree(m);
1419   dstr_destroy(&d);
1420   return (ok);
1421 }
1422
1423 static int twocl(dstr *v)
1424 {
1425   dstr d = DSTR_INIT;
1426   mpw *m, *ml;
1427   size_t sz;
1428   int ok = 1;
1429
1430   sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1431   dstr_ensure(&d, sz);
1432
1433   sz = MPW_RQ(sz);
1434   m = xmalloc(MPWS(sz));
1435   ml = m + sz;
1436
1437   mpx_loadl(m, ml, v[0].buf, v[0].len);
1438   mpx_storel2cn(m, ml, d.buf, v[1].len);
1439   if (memcmp(d.buf, v[1].buf, v[1].len)) {
1440     dumpbits("\n*** storel2cn failed", d.buf, v[1].len);
1441     ok = 0;
1442   }
1443
1444   mpx_loadl2cn(m, ml, v[1].buf, v[1].len);
1445   mpx_storel(m, ml, d.buf, v[0].len);
1446   if (memcmp(d.buf, v[0].buf, v[0].len)) {
1447     dumpbits("\n*** loadl2cn failed", d.buf, v[0].len);
1448     ok = 0;
1449   }
1450
1451   if (!ok) {
1452     dumpbits("pos", v[0].buf, v[0].len);
1453     dumpbits("neg", v[1].buf, v[1].len);
1454   }
1455
1456   xfree(m);
1457   dstr_destroy(&d);
1458
1459   return (ok);
1460 }
1461
1462 static int twocb(dstr *v)
1463 {
1464   dstr d = DSTR_INIT;
1465   mpw *m, *ml;
1466   size_t sz;
1467   int ok = 1;
1468
1469   sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1470   dstr_ensure(&d, sz);
1471
1472   sz = MPW_RQ(sz);
1473   m = xmalloc(MPWS(sz));
1474   ml = m + sz;
1475
1476   mpx_loadb(m, ml, v[0].buf, v[0].len);
1477   mpx_storeb2cn(m, ml, d.buf, v[1].len);
1478   if (memcmp(d.buf, v[1].buf, v[1].len)) {
1479     dumpbits("\n*** storeb2cn failed", d.buf, v[1].len);
1480     ok = 0;
1481   }
1482
1483   mpx_loadb2cn(m, ml, v[1].buf, v[1].len);
1484   mpx_storeb(m, ml, d.buf, v[0].len);
1485   if (memcmp(d.buf, v[0].buf, v[0].len)) {
1486     dumpbits("\n*** loadb2cn failed", d.buf, v[0].len);
1487     ok = 0;
1488   }
1489
1490   if (!ok) {
1491     dumpbits("pos", v[0].buf, v[0].len);
1492     dumpbits("neg", v[1].buf, v[1].len);
1493   }
1494
1495   xfree(m);
1496   dstr_destroy(&d);
1497
1498   return (ok);
1499 }
1500
1501 static int lsl(dstr *v)
1502 {
1503   mpw *a, *al;
1504   int n = *(int *)v[1].buf;
1505   mpw *c, *cl;
1506   mpw *d, *dl;
1507   int ok = 1;
1508
1509   LOAD(a, al, &v[0]);
1510   LOAD(c, cl, &v[2]);
1511   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1512
1513   mpx_lsl(d, dl, a, al, n);
1514   if (!mpx_ueq(d, dl, c, cl)) {
1515     fprintf(stderr, "\n*** lsl(%i) failed\n", n);
1516     dumpmp("       a", a, al);
1517     dumpmp("expected", c, cl);
1518     dumpmp("  result", d, dl);
1519     ok = 0;
1520   }
1521
1522   xfree(a); xfree(c); xfree(d);
1523   return (ok);
1524 }
1525
1526 static int lslc(dstr *v)
1527 {
1528   mpw *a, *al;
1529   int n = *(int *)v[1].buf;
1530   mpw *c, *cl;
1531   mpw *d, *dl;
1532   int ok = 1;
1533
1534   LOAD(a, al, &v[0]);
1535   LOAD(c, cl, &v[2]);
1536   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1537
1538   mpx_lslc(d, dl, a, al, n);
1539   if (!mpx_ueq(d, dl, c, cl)) {
1540     fprintf(stderr, "\n*** lslc(%i) failed\n", n);
1541     dumpmp("       a", a, al);
1542     dumpmp("expected", c, cl);
1543     dumpmp("  result", d, dl);
1544     ok = 0;
1545   }
1546
1547   xfree(a); xfree(c); xfree(d);
1548   return (ok);
1549 }
1550
1551 static int lsr(dstr *v)
1552 {
1553   mpw *a, *al;
1554   int n = *(int *)v[1].buf;
1555   mpw *c, *cl;
1556   mpw *d, *dl;
1557   int ok = 1;
1558
1559   LOAD(a, al, &v[0]);
1560   LOAD(c, cl, &v[2]);
1561   ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1562
1563   mpx_lsr(d, dl, a, al, n);
1564   if (!mpx_ueq(d, dl, c, cl)) {
1565     fprintf(stderr, "\n*** lsr(%i) failed\n", n);
1566     dumpmp("       a", a, al);
1567     dumpmp("expected", c, cl);
1568     dumpmp("  result", d, dl);
1569     ok = 0;
1570   }
1571
1572   xfree(a); xfree(c); xfree(d);
1573   return (ok);
1574 }
1575
1576 static int uadd(dstr *v)
1577 {
1578   mpw *a, *al;
1579   mpw *b, *bl;
1580   mpw *c, *cl;
1581   mpw *d, *dl;
1582   int ok = 1;
1583
1584   LOAD(a, al, &v[0]);
1585   LOAD(b, bl, &v[1]);
1586   LOAD(c, cl, &v[2]);
1587   ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1588
1589   mpx_uadd(d, dl, a, al, b, bl);
1590   if (!mpx_ueq(d, dl, c, cl)) {
1591     fprintf(stderr, "\n*** uadd failed\n");
1592     dumpmp("       a", a, al);
1593     dumpmp("       b", b, bl);
1594     dumpmp("expected", c, cl);
1595     dumpmp("  result", d, dl);
1596     ok = 0;
1597   }
1598
1599   xfree(a); xfree(b); xfree(c); xfree(d);
1600   return (ok);
1601 }
1602
1603 static int usub(dstr *v)
1604 {
1605   mpw *a, *al;
1606   mpw *b, *bl;
1607   mpw *c, *cl;
1608   mpw *d, *dl;
1609   int ok = 1;
1610
1611   LOAD(a, al, &v[0]);
1612   LOAD(b, bl, &v[1]);
1613   LOAD(c, cl, &v[2]);
1614   ALLOC(d, dl, al - a);
1615
1616   mpx_usub(d, dl, a, al, b, bl);
1617   if (!mpx_ueq(d, dl, c, cl)) {
1618     fprintf(stderr, "\n*** usub failed\n");
1619     dumpmp("       a", a, al);
1620     dumpmp("       b", b, bl);
1621     dumpmp("expected", c, cl);
1622     dumpmp("  result", d, dl);
1623     ok = 0;
1624   }
1625
1626   xfree(a); xfree(b); xfree(c); xfree(d);
1627   return (ok);
1628 }
1629
1630 static int umul(dstr *v)
1631 {
1632   mpw *a, *al;
1633   mpw *b, *bl;
1634   mpw *c, *cl;
1635   mpw *d, *dl;
1636   int ok = 1;
1637
1638   LOAD(a, al, &v[0]);
1639   LOAD(b, bl, &v[1]);
1640   LOAD(c, cl, &v[2]);
1641   ALLOC(d, dl, (al - a) + (bl - b));
1642
1643   mpx_umul(d, dl, a, al, b, bl);
1644   if (!mpx_ueq(d, dl, c, cl)) {
1645     fprintf(stderr, "\n*** umul failed\n");
1646     dumpmp("       a", a, al);
1647     dumpmp("       b", b, bl);
1648     dumpmp("expected", c, cl);
1649     dumpmp("  result", d, dl);
1650     ok = 0;
1651   }
1652
1653   xfree(a); xfree(b); xfree(c); xfree(d);
1654   return (ok);
1655 }
1656
1657 static int usqr(dstr *v)
1658 {
1659   mpw *a, *al;
1660   mpw *c, *cl;
1661   mpw *d, *dl;
1662   int ok = 1;
1663
1664   LOAD(a, al, &v[0]);
1665   LOAD(c, cl, &v[1]);
1666   ALLOC(d, dl, 2 * (al - a));
1667
1668   mpx_usqr(d, dl, a, al);
1669   if (!mpx_ueq(d, dl, c, cl)) {
1670     fprintf(stderr, "\n*** usqr failed\n");
1671     dumpmp("       a", a, al);
1672     dumpmp("expected", c, cl);
1673     dumpmp("  result", d, dl);
1674     ok = 0;
1675   }
1676
1677   xfree(a); xfree(c); xfree(d);
1678   return (ok);
1679 }
1680
1681 static int udiv(dstr *v)
1682 {
1683   mpw *a, *al;
1684   mpw *b, *bl;
1685   mpw *q, *ql;
1686   mpw *r, *rl;
1687   mpw *qq, *qql;
1688   mpw *s, *sl;
1689   int ok = 1;
1690
1691   ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1692   LOAD(b, bl, &v[1]);
1693   LOAD(q, ql, &v[2]);
1694   LOAD(r, rl, &v[3]);
1695   ALLOC(qq, qql, al - a);
1696   ALLOC(s, sl, (bl - b) + 1);
1697
1698   mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1699   if (!mpx_ueq(qq, qql, q, ql) ||
1700       !mpx_ueq(a, al, r, rl)) {
1701     fprintf(stderr, "\n*** udiv failed\n");
1702     dumpmp(" divisor", b, bl);
1703     dumpmp("expect r", r, rl);
1704     dumpmp("result r", a, al);
1705     dumpmp("expect q", q, ql);
1706     dumpmp("result q", qq, qql);
1707     ok = 0;
1708   }
1709
1710   xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq);
1711   return (ok);
1712 }
1713
1714 static test_chunk defs[] = {
1715   { "load-store", loadstore, { &type_hex, 0 } },
1716   { "2cl", twocl, { &type_hex, &type_hex, } },
1717   { "2cb", twocb, { &type_hex, &type_hex, } },
1718   { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
1719   { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
1720   { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1721   { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1722   { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1723   { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1724   { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1725   { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1726   { 0, 0, { 0 } }
1727 };
1728
1729 int main(int argc, char *argv[])
1730 {
1731   test_run(argc, argv, defs, SRCDIR"/t/mpx");
1732   return (0);
1733 }
1734
1735 #endif
1736
1737 /*----- That's all, folks -------------------------------------------------*/