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

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