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

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