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

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