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

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