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

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