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

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