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

   1 /* -*-c-*-
   2  *
   3  * $Id: mpx.h,v 1.17 2004/03/27 00:04:46 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.h,v $
  33  * Revision 1.17  2004/03/27 00:04:46  mdw
  34  * Implement efficient reduction for pleasant-looking primes.
  35  *
  36  * Revision 1.16  2003/05/16 09:09:24  mdw
  37  * Fix @mp_lsl2c@.  Turns out to be surprisingly tricky.
  38  *
  39  * Revision 1.15  2002/10/19 17:56:50  mdw
  40  * Fix bit operations.  Test them (a bit) better.
  41  *
  42  * Revision 1.14  2002/10/09 00:36:03  mdw
  43  * Fix bounds on workspace for Karatsuba operations.
  44  *
  45  * Revision 1.13  2002/10/06 22:52:50  mdw
  46  * Pile of changes for supporting two's complement properly.
  47  *
  48  * Revision 1.12  2001/04/03 19:36:05  mdw
  49  * Add some simple bitwise operations so that Perl can use them.
  50  *
  51  * Revision 1.11  2000/10/08 15:48:35  mdw
  52  * Rename Karatsuba constants now that we have @gfx_kmul@ too.
  53  *
  54  * Revision 1.10  2000/10/08 12:06:12  mdw
  55  * Provide @mpx_ueq@ for rapidly testing equality of two integers.
  56  *
  57  * Revision 1.9  1999/12/22 15:49:07  mdw
  58  * New function for division by a small integer.
  59  *
  60  * Revision 1.8  1999/12/11 10:57:43  mdw
  61  * Karatsuba squaring algorithm.
  62  *
  63  * Revision 1.7  1999/12/11 01:51:28  mdw
  64  * Change Karatsuba parameters slightly.
  65  *
  66  * Revision 1.6  1999/12/10 23:23:51  mdw
  67  * Karatsuba-Ofman multiplication algorithm.
  68  *
  69  * Revision 1.5  1999/11/20 22:23:27  mdw
  70  * Add function versions of some low-level macros with wider use.
  71  *
  72  * Revision 1.4  1999/11/17 18:04:43  mdw
  73  * Add two's complement support.  Fix a bug in MPX_UMLAN.
  74  *
  75  * Revision 1.3  1999/11/13 01:51:29  mdw
  76  * Minor interface changes.  Should be stable now.
  77  *
  78  * Revision 1.2  1999/11/11 17:47:55  mdw
  79  * Minor changes for different `mptypes.h' format.
  80  *
  81  * Revision 1.1  1999/09/03 08:41:12  mdw
  82  * Initial import.
  83  *
  84  */
  85
  86 #ifndef CATACOMB_MPX_H
  87 #define CATACOMB_MPX_H
  88
  89 #ifdef __cplusplus
  90   extern "C" {
  91 #endif
  92
  93 /*----- The idea ----------------------------------------------------------*
  94  *
  95  * This file provides functions and macros which work on vectors of words as
  96  * unsigned multiprecision integers.  The interface works in terms of base
  97  * and limit pointers (i.e., a pointer to the start of a vector, and a
  98  * pointer just past its end) rather than base pointer and length, because
  99  * that requires more arithmetic and state to work on.
 100  *
 101  * The interfaces are slightly bizarre in other ways.  Try to use the
 102  * higher-level functions where you can: they're rather better designed to
 103  * actually be friendly and useful.
 104  */
 105
 106 /*----- Header files ------------------------------------------------------*/
 107
 108 #include <string.h>
 109
 110 #ifndef CATACOMB_MPW_H
 111 #  include "mpw.h"
 112 #endif
 113
 114 /*----- General manipulation ----------------------------------------------*/
 115
 116 /* --- @MPX_SHRINK@ --- *
 117  *
 118  * Arguments:   @const mpw *v@ = pointer to vector of words
 119  *              @const mpw *vl@ = (updated) current limit of vector
 120  *
 121  * Use:         Shrinks down the limit of a multiprecision integer vector.
 122  */
 123
 124 #define MPX_SHRINK(v, vl) do {                                          \
 125   const mpw *_vv = (v), *_vvl = (vl);                                   \
 126   while (_vvl > _vv && !_vvl[-1])                                       \
 127     _vvl--;                                                             \
 128   (vl) = (mpw *)_vvl;                                                   \
 129 } while (0)
 130
 131 /* --- @MPX_BITS@ --- *
 132  *
 133  * Arguments:   @unsigned long b@ = result variable
 134  *              @const mpw *v@ = pointer to array of words
 135  *              @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
 136  *
 137  * Use:         Calculates the number of bits in a multiprecision value.
 138  */
 139
 140 #define MPX_BITS(b, v, vl) do {                                         \
 141   const mpw *_v = (v), *_vl = (vl);                                     \
 142   MPX_SHRINK(_v, _vl);                                                  \
 143   if (_v == _vl)                                                        \
 144     (b) = 0;                                                            \
 145   else {                                                                \
 146     unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1;                   \
 147     mpw _w = _vl[-1];                                                   \
 148     unsigned _k = MPW_BITS / 2;                                         \
 149     while (_k) {                                                        \
 150       if (_w >> _k) {                                                   \
 151         _w >>= _k;                                                      \
 152         _b += _k;                                                       \
 153       }                                                                 \
 154       _k >>= 1;                                                         \
 155     }                                                                   \
 156     (b) = _b;                                                           \
 157   }                                                                     \
 158 } while (0)
 159
 160 /* --- @MPX_OCTETS@ --- *
 161  *
 162  * Arguments:   @size_t o@ = result variable
 163  *              @const mpw *v, *vl@ = pointer to array of words
 164  *
 165  * Use:         Calculates the number of octets in a multiprecision value.
 166  */
 167
 168 #define MPX_OCTETS(o, v, vl) do {                                       \
 169   unsigned long _bb;                                                    \
 170   MPX_BITS(_bb, (v), (vl));                                             \
 171   (o) = (_bb + 7) >> 3;                                                 \
 172 } while (0)
 173
 174 /* --- @MPX_OCTETS2C@ --- *
 175  *
 176  * Arguments:   @size_t o@ = result variable
 177  *              @const mpw *v, *vl@ = pointer to array of words
 178  *
 179  * Use:         Calculates the number of octets in a multiprecision value, if
 180  *              you represent it as two's complement.
 181  */
 182
 183 #define MPX_OCTETS2C(o, v, vl) do {                                     \
 184   unsigned long _bb;                                                    \
 185   MPX_BITS(_bb, (v), (vl));                                             \
 186   (o) = (_bb >> 3) + 1;                                                 \
 187 } while (0)
 188
 189 /* --- @MPX_COPY@ --- *
 190  *
 191  * Arguments:   @dv, dvl@ = destination vector base and limit
 192  *              @av, avl@ = source vector base and limit
 193  *
 194  * Use:         Copies a multiprecision integer.
 195  */
 196
 197 #define MPX_COPY(dv, dvl, av, avl) do {                                 \
 198   mpw *_dv = (dv), *_dvl = (dvl);                                       \
 199   size_t _dn = _dvl - _dv;                                              \
 200   const mpw *_av = (av), *_avl = (avl);                                 \
 201   size_t _an = _avl - _av;                                              \
 202   if (_av == _dv) {                                                     \
 203     if (_dvl > _avl)                                                    \
 204       memset(_dv, 0, MPWS(_dn - _an));                                  \
 205   } else if (_an >= _dn)                                                \
 206     memmove(_dv, _av, MPWS(_dn));                                       \
 207   else {                                                                \
 208     memmove(_dv, _av, MPWS(_an));                                       \
 209     memset(_dv + _an, 0, MPWS(_dn - _an));                              \
 210   }                                                                     \
 211 } while (0)
 212
 213 /* --- @MPX_ZERO@ --- *
 214  *
 215  * Arguments:   @v, vl@ = base and limit of vector to clear
 216  *
 217  * Use:         Zeroes the area between the two vector pointers.
 218  */
 219
 220 #define MPX_ZERO(v, vl) do {                                            \
 221   mpw *_v = (v), *_vl = (vl);                                           \
 222   if (_v < _vl)                                                         \
 223     memset(_v, 0, MPWS(_vl - _v));                                      \
 224 } while (0)
 225
 226 /* --- @MPX_ONE@ --- *
 227  *
 228  * Arguments:   @v, vl@ = base and limit of vector to clear
 229  *
 230  * Use:         Fills the area between the two vector pointers with ones.
 231  */
 232
 233 #define MPX_ONE(v, vl) do {                                             \
 234   mpw * _v = (v);                                                       \
 235   const mpw *_vl = (vl);                                                \
 236   while (_v < _vl)                                                      \
 237     *_v++ = MPW_MAX;                                                    \
 238 } while (0)
 239
 240 /*----- Loading and storing -----------------------------------------------*/
 241
 242 /* --- @mpx_storel@ --- *
 243  *
 244  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 245  *              @void *p@ = pointer to octet array
 246  *              @size_t sz@ = size of octet array
 247  *
 248  * Returns:     ---
 249  *
 250  * Use:         Stores an MP in an octet array, least significant octet
 251  *              first.  High-end octets are silently discarded if there
 252  *              isn't enough space for them.
 253  */
 254
 255 extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
 256                        void */*p*/, size_t /*sz*/);
 257
 258 /* --- @mpx_loadl@ --- *
 259  *
 260  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 261  *              @const void *p@ = pointer to octet array
 262  *              @size_t sz@ = size of octet array
 263  *
 264  * Returns:     ---
 265  *
 266  * Use:         Loads an MP in an octet array, least significant octet
 267  *              first.  High-end octets are ignored if there isn't enough
 268  *              space for them.
 269  */
 270
 271 extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
 272                       const void */*p*/, size_t /*sz*/);
 273
 274 /* --- @mpx_storeb@ --- *
 275  *
 276  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 277  *              @void *p@ = pointer to octet array
 278  *              @size_t sz@ = size of octet array
 279  *
 280  * Returns:     ---
 281  *
 282  * Use:         Stores an MP in an octet array, most significant octet
 283  *              first.  High-end octets are silently discarded if there
 284  *              isn't enough space for them.
 285  */
 286
 287 extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
 288                        void */*p*/, size_t /*sz*/);
 289
 290 /* --- @mpx_loadb@ --- *
 291  *
 292  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 293  *              @const void *p@ = pointer to octet array
 294  *              @size_t sz@ = size of octet array
 295  *
 296  * Returns:     ---
 297  *
 298  * Use:         Loads an MP in an octet array, most significant octet
 299  *              first.  High-end octets are ignored if there isn't enough
 300  *              space for them.
 301  */
 302
 303 extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
 304                       const void */*p*/, size_t /*sz*/);
 305
 306 /* --- @mpx_storel2cn@ --- *
 307  *
 308  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 309  *              @void *pp@ = pointer to octet array
 310  *              @size_t sz@ = size of octet array
 311  *
 312  * Returns:     ---
 313  *
 314  * Use:         Stores a negative MP in an octet array, least significant
 315  *              octet first, as two's complement.  High-end octets are
 316  *              silently discarded if there isn't enough space for them.
 317  *              This obviously makes the output bad.
 318  */
 319
 320 extern void mpx_storel2cn(const mpw */*v*/, const mpw */*vl*/,
 321                           void */*p*/, size_t /*sz*/);
 322
 323 /* --- @mpx_loadl2cn@ --- *
 324  *
 325  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 326  *              @const void *pp@ = pointer to octet array
 327  *              @size_t sz@ = size of octet array
 328  *
 329  * Returns:     ---
 330  *
 331  * Use:         Loads a negative MP in an octet array, least significant
 332  *              octet first, as two's complement.  High-end octets are
 333  *              ignored if there isn't enough space for them.  This probably
 334  *              means you made the wrong choice coming here.
 335  */
 336
 337 extern void mpx_loadl2cn(mpw */*v*/, mpw */*vl*/,
 338                          const void */*p*/, size_t /*sz*/);
 339
 340 /* --- @mpx_storeb2cn@ --- *
 341  *
 342  * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 343  *              @void *pp@ = pointer to octet array
 344  *              @size_t sz@ = size of octet array
 345  *
 346  * Returns:     ---
 347  *
 348  * Use:         Stores a negative MP in an octet array, most significant
 349  *              octet first, as two's complement.  High-end octets are
 350  *              silently discarded if there isn't enough space for them,
 351  *              which probably isn't what you meant.
 352  */
 353
 354 extern void mpx_storeb2cn(const mpw */*v*/, const mpw */*vl*/,
 355                           void */*p*/, size_t /*sz*/);
 356
 357 /* --- @mpx_loadb2cn@ --- *
 358  *
 359  * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 360  *              @const void *pp@ = pointer to octet array
 361  *              @size_t sz@ = size of octet array
 362  *
 363  * Returns:     ---
 364  *
 365  * Use:         Loads a negative MP in an octet array, most significant octet
 366  *              first as two's complement.  High-end octets are ignored if
 367  *              there isn't enough space for them.  This probably means you
 368  *              chose this function wrongly.
 369  */
 370
 371 extern void mpx_loadb2cn(mpw */*v*/, mpw */*vl*/,
 372                          const void */*p*/, size_t /*sz*/);
 373
 374
 375 /*----- Logical shifting --------------------------------------------------*/
 376
 377 /* --- @mpx_lsl@ --- *
 378  *
 379  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 380  *              @const mpw *av, *avl@ = source vector base and limit
 381  *              @size_t n@ = number of bit positions to shift by
 382  *
 383  * Returns:     ---
 384  *
 385  * Use:         Performs a logical shift left operation on an integer.
 386  */
 387
 388 extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/,
 389                     const mpw */*av*/, const mpw */*avl*/,
 390                     size_t /*n*/);
 391
 392 /* --- @mpx_lslc@ --- *
 393  *
 394  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 395  *              @const mpw *av, *avl@ = source vector base and limit
 396  *              @size_t n@ = number of bit positions to shift by
 397  *
 398  * Returns:     ---
 399  *
 400  * Use:         Performs a logical shift left operation on an integer, only
 401  *              it fills in the bits with ones instead of zeroes.
 402  */
 403
 404 extern void mpx_lslc(mpw */*dv*/, mpw */*dvl*/,
 405                      const mpw */*av*/, const mpw */*avl*/,
 406                      size_t /*n*/);
 407
 408 /* --- @mpx_lsr@ --- *
 409  *
 410  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 411  *              @const mpw *av, *avl@ = source vector base and limit
 412  *              @size_t n@ = number of bit positions to shift by
 413  *
 414  * Returns:     ---
 415  *
 416  * Use:         Performs a logical shift right operation on an integer.
 417  */
 418
 419 extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
 420                     const mpw */*av*/, const mpw */*avl*/,
 421                     size_t /*n*/);
 422
 423 /*----- Bitwise operations ------------------------------------------------*/
 424
 425 /* --- @mpx_bitop@ --- *
 426  *
 427  * Arguments:   @mpw *dv, *dvl@ = destination vector
 428  *              @const mpw *av, *avl@ = first source vector
 429  *              @const mpw *bv, *bvl@ = second source vector
 430  *
 431  * Returns:     ---
 432  *
 433  * Use:         Provide the dyadic boolean functions.  The functions are
 434  *              named after the truth table they generate:
 435  *
 436  *                      a:      0011
 437  *                      b:      0101
 438  *                      @mpx_bitXXXX@
 439  */
 440
 441 #define MPX_DOBIN(what)                                                 \
 442   what(0000) what(0001) what(0010) what(0011)                           \
 443   what(0100) what(0101) what(0110) what(0111)                           \
 444   what(1000) what(1001) what(1010) what(1011)                           \
 445   what(1100) what(1101) what(1110) what(1111)
 446
 447 #define MPX_BITDECL(string)                                             \
 448   extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/,                \
 449                               const mpw */*av*/, const mpw */*avl*/,    \
 450                               const mpw */*bv*/, const mpw */*bvl*/);
 451 MPX_DOBIN(MPX_BITDECL)
 452
 453 /* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- *
 454  *
 455  * Synonyms for the commonly-used functions above.
 456  */
 457
 458 #define mpx_and  mpx_bit0001
 459 #define mpx_or   mpx_bit0111
 460 #define mpx_nand mpx_bit1110
 461 #define mpx_nor  mpx_bit1000
 462 #define mpx_xor  mpx_bit0110
 463
 464 /* --- @mpx_not@ --- *
 465  *
 466  * Arguments:   @mpw *dv, *dvl@ = destination vector
 467  *              @const mpw *av, *avl@ = first source vector
 468  *
 469  * Returns:     ---
 470  *
 471  * Use;         Bitwise NOT.
 472  */
 473
 474 extern void mpx_not(mpw */*dv*/, mpw */*dvl*/,
 475                     const mpw */*av*/, const mpw */*avl*/);
 476
 477 /*----- Unsigned arithmetic -----------------------------------------------*/
 478
 479 /* --- @mpx_2c@ --- *
 480  *
 481  * Arguments:   @mpw *dv, *dvl@ = destination vector
 482  *              @const mpw *v, *vl@ = source vector
 483  *
 484  * Returns:     ---
 485  *
 486  * Use:         Calculates the two's complement of @v@.
 487  */
 488
 489 extern void mpx_2c(mpw */*dv*/, mpw */*dvl*/,
 490                    const mpw */*v*/, const mpw */*vl*/);
 491
 492 /* --- @mpx_ueq@ --- *
 493  *
 494  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 495  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 496  *
 497  * Returns:     Nonzero if the two vectors are equal.
 498  *
 499  * Use:         Performs an unsigned integer test for equality.
 500  */
 501
 502 extern int mpx_ueq(const mpw */*av*/, const mpw */*avl*/,
 503                    const mpw */*bv*/, const mpw */*bvl*/);
 504
 505 /* --- @mpx_ucmp@ --- *
 506  *
 507  * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 508  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 509  *
 510  * Returns:     Less than, equal to, or greater than zero depending on
 511  *              whether @a@ is less than, equal to or greater than @b@,
 512  *              respectively.
 513  *
 514  * Use:         Performs an unsigned integer comparison.
 515  */
 516
 517 #define MPX_UCMP(av, avl, op, dv, dvl)                                  \
 518   (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
 519
 520 extern int mpx_ucmp(const mpw */*av*/, const mpw */*avl*/,
 521                     const mpw */*bv*/, const mpw */*bvl*/);
 522
 523 /* --- @mpx_uadd@ --- *
 524  *
 525  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 526  *              @const mpw *av, *avl@ = first addend vector base and limit
 527  *              @const mpw *bv, *bvl@ = second addend vector base and limit
 528  *
 529  * Returns:     ---
 530  *
 531  * Use:         Performs unsigned integer addition.  If the result overflows
 532  *              the destination vector, high-order bits are discarded.  This
 533  *              means that two's complement addition happens more or less for
 534  *              free, although that's more a side-effect than anything else.
 535  *              The result vector may be equal to either or both source
 536  *              vectors, but may not otherwise overlap them.
 537  */
 538
 539 extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
 540                      const mpw */*av*/, const mpw */*avl*/,
 541                      const mpw */*bv*/, const mpw */*bvl*/);
 542
 543 /* --- @mpx_uaddn@ --- *
 544  *
 545  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 546  *              @mpw n@ = other addend
 547  *
 548  * Returns:     ---
 549  *
 550  * Use:         Adds a small integer to a multiprecision number.
 551  */
 552
 553 #define MPX_UADDN(dv, dvl, n) do {                                      \
 554   mpw *_ddv = (dv), *_ddvl = (dvl);                                     \
 555   mpw _c = (n);                                                         \
 556                                                                         \
 557   while (_c && _ddv < _ddvl) {                                          \
 558     mpd _x = (mpd)*_ddv + (mpd)_c;                                      \
 559     *_ddv++ = MPW(_x);                                                  \
 560     _c = _x >> MPW_BITS;                                                \
 561   }                                                                     \
 562 } while (0)
 563
 564 extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
 565
 566 /* --- @mpx_uaddnlsl@ --- *
 567  *
 568  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 569  *              @mpw a@ = second argument
 570  *              @unsigned o@ = offset in bits
 571  *
 572  * Returns:     ---
 573  *
 574  * Use:         Computes %$d + 2^o a$%.  If the result overflows then
 575  *              high-order bits are discarded, as usual.  We must have
 576  *              @0 < o < MPW_BITS@.
 577  */
 578
 579 extern void mpx_uaddnlsl(mpw */*dv*/, mpw */*dvl*/,
 580                          mpw /*a*/, unsigned /*o*/);
 581
 582 /* --- @mpx_usub@ --- *
 583  *
 584  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 585  *              @const mpw *av, *avl@ = first argument vector base and limit
 586  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 587  *
 588  * Returns:     ---
 589  *
 590  * Use:         Performs unsigned integer subtraction.  If the result
 591  *              overflows the destination vector, high-order bits are
 592  *              discarded.  This means that two's complement subtraction
 593  *              happens more or less for free, although that's more a side-
 594  *              effect than anything else.  The result vector may be equal to
 595  *              either or both source vectors, but may not otherwise overlap
 596  *              them.
 597  */
 598
 599 extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
 600                      const mpw */*av*/, const mpw */*avl*/,
 601                      const mpw */*bv*/, const mpw */*bvl*/);
 602
 603 /* --- @mpx_usubn@ --- *
 604  *
 605  * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 606  *              @n@ = subtrahend
 607  *
 608  * Returns:     ---
 609  *
 610  * Use:         Subtracts a small integer from a multiprecision number.
 611  */
 612
 613 #define MPX_USUBN(dv, dvl, n) do {                                      \
 614   mpw *_ddv = (dv), *_ddvl = (dvl);                                     \
 615   mpw _c = (n);                                                         \
 616                                                                         \
 617   while (_ddv < _ddvl) {                                                \
 618     mpd _x = (mpd)*_ddv - (mpd)_c;                                      \
 619     *_ddv++ = MPW(_x);                                                  \
 620     if (_x >> MPW_BITS)                                                 \
 621       _c = 1;                                                           \
 622     else                                                                \
 623       break;                                                            \
 624   }                                                                     \
 625 } while (0)
 626
 627 extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
 628
 629 /* --- @mpx_usubnlsl@ --- *
 630  *
 631  * Arguments:   @mpw *dv, *dvl@ = destination and first argument vector
 632  *              @mpw a@ = second argument
 633  *              @unsigned o@ = offset in bits
 634  *
 635  * Returns:     ---
 636  *
 637  * Use:         Computes %$d - 2^o a$%.  If the result overflows then
 638  *              high-order bits are discarded, as usual, so you get two's
 639  *              complement.  Which might be what you wanted...  We must have
 640  *              @0 < o < MPW_BITS@.
 641  */
 642
 643 extern void mpx_usubnlsl(mpw */*dv*/, mpw */*dvl*/,
 644                          mpw /*a*/, unsigned /*o*/);
 645
 646 /* --- @mpx_umul@ --- *
 647  *
 648  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 649  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 650  *              @const mpw *bv, *bvl@ = multiplier vector base and limit
 651  *
 652  * Returns:     ---
 653  *
 654  * Use:         Performs unsigned integer multiplication.  If the result
 655  *              overflows the desination vector, high-order bits are
 656  *              discarded.  The result vector may not overlap the argument
 657  *              vectors in any way.
 658  */
 659
 660 extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
 661                      const mpw */*av*/, const mpw */*avl*/,
 662                      const mpw */*bv*/, const mpw */*bvl*/);
 663
 664 /* --- @mpx_umuln@ --- *
 665  *
 666  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 667  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 668  *              @mpw m@ = multiplier
 669  *
 670  * Returns:     ---
 671  *
 672  * Use:         Multiplies a multiprecision integer by a single-word value.
 673  *              The destination and source may be equal.  The destination
 674  *              is completely cleared after use.
 675  */
 676
 677 #define MPX_UMULN(dv, dvl, av, avl, m) do {                             \
 678   mpw *_dv = (dv), *_dvl = (dvl);                                       \
 679   const mpw *_av = (av), *_avl = (avl);                                 \
 680   mpw _c = 0;                                                           \
 681   mpd _m = (m);                                                         \
 682                                                                         \
 683   while (_av < _avl) {                                                  \
 684     mpd _x;                                                             \
 685     if (_dv >= _dvl)                                                    \
 686       break;                                                            \
 687     _x = (mpd)_m * (mpd)*_av++ + _c;                                    \
 688     *_dv++ = MPW(_x);                                                   \
 689     _c = _x >> MPW_BITS;                                                \
 690   }                                                                     \
 691   if (_dv < _dvl) {                                                     \
 692     *_dv++ = MPW(_c);                                                   \
 693     MPX_ZERO(_dv, _dvl);                                                \
 694   }                                                                     \
 695 } while (0)
 696
 697 extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/,
 698                       const mpw */*av*/, const mpw */*avl*/, mpw m);
 699
 700 /* --- @mpx_umlan@ --- *
 701  *
 702  * Arguments:   @mpw *dv, *dvl@ = destination/accumulator base and limit
 703  *              @const mpw *av, *avl@ = multiplicand vector base and limit
 704  *              @mpw m@ = multiplier
 705  *
 706  * Returns:     ---
 707  *
 708  * Use:         Multiplies a multiprecision integer by a single-word value
 709  *              and adds the result to an accumulator.
 710  */
 711
 712 #define MPX_UMLAN(dv, dvl, av, avl, m) do {                             \
 713   mpw *_dv = (dv), *_dvl = (dvl);                                       \
 714   const mpw *_av = (av), *_avl = (avl);                                 \
 715   mpw _cc = 0;                                                          \
 716   mpd _m = (m);                                                         \
 717                                                                         \
 718   while (_dv < _dvl && _av < _avl) {                                    \
 719     mpd _x;                                                             \
 720     _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc;                       \
 721     *_dv++ = MPW(_x);                                                   \
 722     _cc = _x >> MPW_BITS;                                               \
 723   }                                                                     \
 724   MPX_UADDN(_dv, _dvl, _cc);                                            \
 725 } while (0)
 726
 727 extern void mpx_umlan(mpw */*dv*/, mpw */*dvl*/,
 728                       const mpw */*av*/, const mpw */*avl*/, mpw m);
 729
 730 /* --- @mpx_usqr@ --- *
 731  *
 732  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 733  *              @const mpw *av, *av@ = source vector base and limit
 734  *
 735  * Returns:     ---
 736  *
 737  * Use:         Performs unsigned integer squaring.  The result vector must
 738  *              not overlap the source vector in any way.
 739  */
 740
 741 extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
 742                      const mpw */*av*/, const mpw */*avl*/);
 743
 744 /* --- @mpx_udiv@ --- *
 745  *
 746  * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
 747  *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
 748  *              @const mpw *dv, *dvl@ = divisor vector base and limit
 749  *              @mpw *sv, *svl@ = scratch workspace
 750  *
 751  * Returns:     ---
 752  *
 753  * Use:         Performs unsigned integer division.  If the result overflows
 754  *              the quotient vector, high-order bits are discarded.  (Clearly
 755  *              the remainder vector can't overflow.)  The various vectors
 756  *              may not overlap in any way.  Yes, I know it's a bit odd
 757  *              requiring the dividend to be in the result position but it
 758  *              does make some sense really.  The remainder must have
 759  *              headroom for at least two extra words.  The scratch space
 760  *              must be at least one word larger than the divisor.
 761  */
 762
 763 extern void mpx_udiv(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/,
 764                      const mpw */*dv*/, const mpw */*dvl*/,
 765                      mpw */*sv*/, mpw */*svl*/);
 766
 767 /* --- @mpx_udivn@ --- *
 768  *
 769  * Arguments:   @mpw *qv, *qvl@ = storage for the quotient (may overlap
 770  *                      dividend)
 771  *              @const mpw *rv, *rvl@ = dividend
 772  *              @mpw d@ = single-precision divisor
 773  *
 774  * Returns:     Remainder after divison.
 775  *
 776  * Use:         Performs a single-precision division operation.
 777  */
 778
 779 extern mpw mpx_udivn(mpw */*qv*/, mpw */*qvl*/,
 780                      const mpw */*rv*/, const mpw */*rvl*/, mpw /*d*/);
 781
 782 /*----- Karatsuba multiplication algorithms -------------------------------*/
 783
 784 /* --- @MPK_THRESH@ --- *
 785  *
 786  * This is the limiting length for using Karatsuba algorithms.  It's best to
 787  * use the simpler classical multiplication method on numbers smaller than
 788  * this.  It is unsafe to make this constant less than four (i.e., the
 789  * algorithms will fail).
 790  */
 791
 792 #define MPK_THRESH 16
 793
 794 /* --- @mpx_kmul@ --- *
 795  *
 796  * Arguments:   @mpw *dv, *dvl@ = pointer to destination buffer
 797  *              @const mpw *av, *avl@ = pointer to first argument
 798  *              @const mpw *bv, *bvl@ = pointer to second argument
 799  *              @mpw *sv, *svl@ = pointer to scratch workspace
 800  *
 801  * Returns:     ---
 802  *
 803  * Use:         Multiplies two multiprecision integers using Karatsuba's
 804  *              algorithm.  This is rather faster than traditional long
 805  *              multiplication (e.g., @mpx_umul@) on large numbers, although
 806  *              more expensive on small ones.
 807  *
 808  *              The destination must be three times as large as the larger
 809  *              argument.  The scratch space must be five times as large as
 810  *              the larger argument.
 811  */
 812
 813 extern void mpx_kmul(mpw */*dv*/, mpw */*dvl*/,
 814                      const mpw */*av*/, const mpw */*avl*/,
 815                      const mpw */*bv*/, const mpw */*bvl*/,
 816                      mpw */*sv*/, mpw */*svl*/);
 817
 818 /* --- @mpx_ksqr@ --- *
 819  *
 820  * Arguments:   @mpw *dv, *dvl@ = pointer to destination buffer
 821  *              @const mpw *av, *avl@ = pointer to first argument
 822  *              @mpw *sv, *svl@ = pointer to scratch workspace
 823  *
 824  * Returns:     ---
 825  *
 826  * Use:         Squares a multiprecision integers using something similar to
 827  *              Karatsuba's multiplication algorithm.  This is rather faster
 828  *              than traditional long multiplication (e.g., @mpx_umul@) on
 829  *              large numbers, although more expensive on small ones, and
 830  *              rather simpler than full-blown Karatsuba multiplication.
 831  *
 832  *              The destination must be three times as large as the larger
 833  *              argument.  The scratch space must be five times as large as
 834  *              the larger argument.
 835  */
 836
 837 extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
 838                      const mpw */*av*/, const mpw */*avl*/,
 839                      mpw */*sv*/, mpw */*svl*/);
 840
 841 /*----- That's all, folks -------------------------------------------------*/
 842
 843 #ifdef __cplusplus
 844   }
 845 #endif
 846
 847 #endif