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

   1 /* -*-c-*-
   2  *
   3  * Simple 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_MP_H
  29 #define CATACOMB_MP_H
  30
  31 #ifdef __cplusplus
  32   extern "C" {
  33 #endif
  34
  35 /*----- Header files ------------------------------------------------------*/
  36
  37 #include <assert.h>
  38 #include <string.h>
  39
  40 #include <mLib/macros.h>
  41 #include <mLib/sub.h>
  42
  43 #ifndef CATACOMB_MPW_H
  44 #  include "mpw.h"
  45 #endif
  46
  47 #ifndef CATACOMB_ARENA_H
  48 #  include "arena.h"
  49 #endif
  50
  51 #ifndef CATACOMB_MPARENA_H
  52 #  include "mparena.h"
  53 #endif
  54
  55 #ifndef CATACOMB_MPX_H
  56 #  include "mpx.h"
  57 #endif
  58
  59 /*----- Data structures ---------------------------------------------------*/
  60
  61 /* --- A multiprecision integer --- */
  62
  63 typedef struct mp {
  64   mpw *v, *vl;                          /* Vector of digits, current limit */
  65   size_t sz;                            /* Size of digit buffer in words */
  66   mparena *a;                           /* Arena for buffer allocation */
  67   unsigned f;                           /* Flags (see below) */
  68   unsigned ref;                         /* Reference counter */
  69 } mp;
  70
  71 #define MP_NEG 1u                       /* Negative (signed magnitude) */
  72 #define MP_BURN 2u                      /* Secret (viral flag) */
  73 #define MP_CONST 4u                     /* Uses strange memory allocation */
  74 #define MP_UNDEF 8u                     /* Contains nothing interesting */
  75 #define MP_DESTROYED 16u                /* Has been destroyed */
  76
  77 /* --- A factor for simultaneous exponentation --- *
  78  *
  79  * Used by the Montgomery and Barrett exponentiators.
  80  */
  81
  82 typedef struct mp_expfactor {
  83   mp *base;
  84   mp *exp;
  85 } mp_expfactor;
  86
  87 /*----- Useful constants --------------------------------------------------*/
  88
  89 extern mp mp_const[];
  90
  91 #define MP_ZERO  (&mp_const[0])
  92 #define MP_ONE   (&mp_const[1])
  93 #define MP_TWO   (&mp_const[2])
  94 #define MP_THREE (&mp_const[3])
  95 #define MP_FOUR  (&mp_const[4])
  96 #define MP_FIVE  (&mp_const[5])
  97 #define MP_TEN   (&mp_const[6])
  98 #define MP_256   (&mp_const[7])
  99 #define MP_MONE  (&mp_const[8])
 100
 101 #define MP_NEW ((mp *)0)
 102 #define MP_NEWSEC (&mp_const[9])
 103
 104 /*----- Trivial macros ----------------------------------------------------*/
 105
 106 /* --- @MP_LEN@ --- *
 107  *
 108  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 109  *
 110  * Returns:     Length of the integer, in words.
 111  */
 112
 113 #define MP_LEN(m) ((m)->vl - ((m)->v))
 114
 115 /*----- Memory management and reference counting --------------------------*/
 116
 117 /* --- @mp_new@ --- *
 118  *
 119  * Arguments:   @size_t sz@ = size of vector required
 120  *              @unsigned f@ = flags to set
 121  *
 122  * Returns:     Pointer to a new MP structure.
 123  *
 124  * Use:         Allocates a new multiprecision integer.  The data space is
 125  *              allocated from either the standard global or secret arena,
 126  *              depending on the initial flags requested.
 127  */
 128
 129 extern mp *mp_new(size_t /*sz*/, unsigned /*f*/);
 130
 131 /* --- @mp_create@ --- *
 132  *
 133  * Arguments:   @size_t sz@ = size of vector required
 134  *
 135  * Returns:     Pointer to pristine new MP structure with enough memory
 136  *              bolted onto it.
 137  *
 138  * Use:         Creates a new multiprecision integer with indeterminate
 139  *              contents.  The integer has a single reference.
 140  */
 141
 142 extern mp *mp_create(size_t /*sz*/);
 143
 144 /* --- @mp_createsecure@ --- *
 145  *
 146  * Arguments:   @size_t sz@ = size of vector required
 147  *
 148  * Returns:     Pointer to pristine new MP structure with enough memory
 149  *              bolted onto it.
 150  *
 151  * Use:         Creates a new multiprecision integer with indeterminate
 152  *              contents.  The integer has a single reference.  The integer's
 153  *              data space is allocated from the secure arena.  Its burn flag
 154  *              is set.
 155  */
 156
 157 extern mp *mp_createsecure(size_t /*sz*/);
 158
 159 /* --- @mp_build@ --- *
 160  *
 161  * Arguments:   @mp *m@ = pointer to an MP block to fill in
 162  *              @mpw *v@ = pointer to a word array
 163  *              @mpw *vl@ = pointer just past end of array
 164  *
 165  * Returns:     ---
 166  *
 167  * Use:         Creates a multiprecision integer representing some smallish
 168  *              number.  You must provide storage for the number and dispose
 169  *              of it when you've finished with it.  The number is marked as
 170  *              constant while it exists.
 171  */
 172
 173 extern void mp_build(mp */*m*/, mpw */*v*/, mpw */*vl*/);
 174
 175 /* --- @mp_destroy@ --- *
 176  *
 177  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 178  *
 179  * Returns:     ---
 180  *
 181  * Use:         Destroys a multiprecision integer. The reference count isn't
 182  *              checked.  Don't use this function if you don't know what
 183  *              you're doing: use @mp_drop@ instead.
 184  */
 185
 186 extern void mp_destroy(mp */*m*/);
 187
 188 /* --- @mp_copy@ --- *
 189  *
 190  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 191  *
 192  * Returns:     A copy of the given multiprecision integer.
 193  *
 194  * Use:         Copies the given integer.  In fact you just get another
 195  *              reference to the same old one again.
 196  */
 197
 198 extern mp *mp_copy(mp */*m*/);
 199
 200 #define MP_COPY(m) MUFFLE_WARNINGS_EXPR(GCC_WARNING("-Wunused-value"),  \
 201                                         ((m)->ref++, (m)))
 202
 203 /* --- @mp_drop@ --- *
 204  *
 205  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 206  *
 207  * Returns:     ---
 208  *
 209  * Use:         Drops a reference to an integer which isn't wanted any more.
 210  *              If there are no more references, the integer is destroyed.
 211  */
 212
 213 extern void mp_drop(mp */*m*/);
 214
 215 #define MP_DROP(m) do {                                                 \
 216   mp *_mm = (m);                                                        \
 217   _mm->ref--;                                                           \
 218   if (_mm->ref == 0 && !(_mm->f & MP_CONST))                            \
 219     mp_destroy(_mm);                                                    \
 220 } while (0)
 221
 222 /* --- @mp_split@ --- *
 223  *
 224  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 225  *
 226  * Returns:     A reference to the same integer, possibly with a different
 227  *              address.
 228  *
 229  * Use:         Splits off a modifiable version of the integer referred to.
 230  */
 231
 232 extern mp *mp_split(mp */*m*/);
 233
 234 #define MP_SPLIT(m) do {                                                \
 235   mp *_m = (m);                                                         \
 236   if ((_m->f & MP_CONST) || _m->ref > 1) {                              \
 237     size_t _len = MP_LEN(_m);                                           \
 238     mp *_mm = mp_new(_len, _m->f);                                      \
 239     if (!(_m->f & MP_UNDEF))                                            \
 240       memcpy(_mm->v, _m->v, MPWS(_len));                                \
 241     _m->ref--;                                                          \
 242     _m = _mm;                                                           \
 243   }                                                                     \
 244   (m) = _m;                                                             \
 245 } while (0)
 246
 247 /* --- @mp_resize@ --- *
 248  *
 249  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 250  *              @size_t sz@ = new size
 251  *
 252  * Returns:     ---
 253  *
 254  * Use:         Resizes the vector containing the integer's digits.  The new
 255  *              size must be at least as large as the current integer's
 256  *              length.  This isn't really intended for client use.
 257  */
 258
 259 extern void mp_resize(mp */*m*/, size_t /*sz*/);
 260
 261 #define MP_RESIZE(m, ssz) do {                                          \
 262   mp *_m = (m);                                                         \
 263   size_t _sz = (ssz);                                                   \
 264   mparena *_a = (_m->f & MP_BURN) ? MPARENA_SECURE : MPARENA_GLOBAL;    \
 265   mpw *_v;                                                              \
 266   size_t _len = MP_LEN(_m);                                             \
 267   assert(((void)"can't make size less than length", _sz >= _len));      \
 268   _v = mpalloc(_a, _sz);                                                \
 269   if (!(_m->f & MP_UNDEF))                                              \
 270     memcpy(_v, _m->v, MPWS(_len));                                      \
 271   if (_m->f & MP_BURN)                                                  \
 272     memset(_m->v, 0, MPWS(_m->sz));                                     \
 273   mpfree(_m->a, _m->v);                                                 \
 274   _m->a = _a;                                                           \
 275   _m->v = _v;                                                           \
 276   _m->vl = _v + _len;                                                   \
 277 } while (0)
 278
 279 /* --- @mp_ensure@ --- *
 280  *
 281  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 282  *              @size_t sz@ = required size
 283  *
 284  * Returns:     ---
 285  *
 286  * Use:         Ensures that the integer has enough space for @sz@ digits.
 287  *              The value is not changed.
 288  */
 289
 290 extern void mp_ensure(mp */*m*/, size_t /*sz*/);
 291
 292 #define MP_ENSURE(m, ssz) do {                                          \
 293   mp *_m = (m);                                                         \
 294   size_t _ssz = (ssz);                                                  \
 295   size_t _len = MP_LEN(_m);                                             \
 296   if (_ssz >= _len) {                                                   \
 297     if (_ssz > _m->sz)                                                  \
 298       mp_resize(_m, _ssz);                                              \
 299     if (!(_m->f & MP_UNDEF) && _ssz > _len)                             \
 300       memset(_m->vl, 0, MPWS(_ssz - _len));                             \
 301     _m->vl = _m->v + _ssz;                                              \
 302   }                                                                     \
 303 } while (0)
 304
 305 /* --- @mp_dest@ --- *
 306  *
 307  * Arguments:   @mp *m@ = a suggested destination integer
 308  *              @size_t sz@ = size required for result, in digits
 309  *              @unsigned f@ = various flags
 310  *
 311  * Returns:     A pointer to an appropriate destination.
 312  *
 313  * Use:         Converts a suggested destination into a real destination with
 314  *              the required properties.  If the real destination is @d@,
 315  *              then the following properties will hold:
 316  *
 317  *                * @d@ will have exactly one reference.
 318  *
 319  *                * If @m@ is not @MP_NEW@, then the contents of @m@ will not
 320  *                  change, unless @f@ has the @MP_UNDEF@ flag set.
 321  *
 322  *                * If @m@ is not @MP_NEW@, then he reference count of @m@ on
 323  *                  entry is equal to the sum of the counts of @d@ and @m@ on
 324  *                  exit.
 325  *
 326  *                * The size of @d@ will be at least @sz@.
 327  *
 328  *                * If @f@ has the @MP_BURN@ flag set, then @d@ will be
 329  *                  allocated from @MPARENA_SECURE@.
 330  *
 331  *              Understanding this function is crucial to using Catacomb's
 332  *              multiprecision integer library effectively.
 333  */
 334
 335 extern mp *mp_dest(mp */*m*/, size_t /*sz*/, unsigned /*f*/);
 336
 337 #define MP_DEST(m, ssz, f) do {                                         \
 338   mp *_m = (m);                                                         \
 339   size_t _ssz = (ssz);                                                  \
 340   unsigned _f = (f);                                                    \
 341   _m = mp_dest(_m, _ssz, _f);                                           \
 342   (m) = _m;                                                             \
 343 } while (0)
 344
 345 /*----- Size manipulation -------------------------------------------------*/
 346
 347 /* --- @mp_shrink@ --- *
 348  *
 349  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 350  *
 351  * Returns:     ---
 352  *
 353  * Use:         Reduces the recorded length of an integer.  This doesn't
 354  *              reduce the amount of memory used, although it can improve
 355  *              performance a bit.  To reduce memory, use @mp_minimize@
 356  *              instead.  This can't change the value of an integer, and is
 357  *              therefore safe to use even when there are multiple
 358  *              references.
 359  */
 360
 361 extern void mp_shrink(mp */*m*/);
 362
 363 #define MP_SHRINK(m) do {                                               \
 364   mp *_mm = (m);                                                        \
 365   MPX_SHRINK(_mm->v, _mm->vl);                                          \
 366   if (MP_ZEROP(_mm))                                                    \
 367     _mm->f &= ~MP_NEG;                                                  \
 368 } while (0)
 369
 370 /* --- @mp_minimize@ --- *
 371  *
 372  * Arguments:   @mp *m@ = pointer to a multiprecision integer
 373  *
 374  * Returns:     ---
 375  *
 376  * Use:         Reduces the amount of memory an integer uses.  It's best to
 377  *              do this to numbers which aren't going to change in the
 378  *              future.
 379  */
 380
 381 extern void mp_minimize(mp */*m*/);
 382
 383 /*----- Bit scanning ------------------------------------------------------*/
 384
 385 #ifndef CATACOMB_MPSCAN_H
 386 #  include "mpscan.h"
 387 #endif
 388
 389 /* --- @mp_scan@ --- *
 390  *
 391  * Arguments:   @mpscan *sc@ = pointer to bitscanner block
 392  *              @const mp *m@ = pointer to a multiprecision integer
 393  *
 394  * Returns:     ---
 395  *
 396  * Use:         Initializes a bitscanner on a multiprecision integer.
 397  */
 398
 399 extern void mp_scan(mpscan */*sc*/, const mp */*m*/);
 400
 401 #define MP_SCAN(sc, m) do {                                             \
 402   const mp *_mm = (m);                                                  \
 403   mpscan *_sc = (sc);                                                   \
 404   MPSCAN_INITX(_sc, _mm->v, _mm->vl);                                   \
 405 } while (0)
 406
 407 /* --- @mp_rscan@ --- *
 408  *
 409  * Arguments:   @mpscan *sc@ = pointer to bitscanner block
 410  *              @const mp *m@ = pointer to a multiprecision integer
 411  *
 412  * Returns:     ---
 413  *
 414  * Use:         Initializes a reverse bitscanner on a multiprecision
 415  *              integer.
 416  */
 417
 418 extern void mp_rscan(mpscan */*sc*/, const mp */*m*/);
 419
 420 #define MP_RSCAN(sc, m) do {                                            \
 421   const mp *_mm = (m);                                                  \
 422   mpscan *_sc = (sc);                                                   \
 423   MPSCAN_RINITX(_sc, _mm->v, _mm->vl);                                  \
 424 } while (0)
 425
 426 /* --- Other bitscanning aliases --- */
 427
 428 #define mp_step mpscan_step
 429 #define mp_bit mpscan_bit
 430 #define mp_rstep mpscan_rstep
 431 #define mp_rbit mpscan_rbit
 432
 433 #define MP_STEP MPSCAN_STEP
 434 #define MP_BIT MPSCAN_BIT
 435 #define MP_RSTEP MPSCAN_RSTEP
 436 #define MP_RBIT MPSCAN_RBIT
 437
 438 /*----- Loading and storing -----------------------------------------------*/
 439
 440 /* --- @mp_octets@ --- *
 441  *
 442  * Arguments:   @const mp *m@ = a multiprecision integer
 443  *
 444  * Returns:     The number of octets required to represent @m@.
 445  *
 446  * Use:         Calculates the external storage required for a multiprecision
 447  *              integer.
 448  */
 449
 450 extern size_t mp_octets(const mp */*m*/);
 451
 452 /* --- @mp_octets2c@ --- *
 453  *
 454  * Arguments:   @const mp *m@ = a multiprecision integer
 455  *
 456  * Returns:     The number of octets required to represent @m@.
 457  *
 458  * Use:         Calculates the external storage required for a multiprecision
 459  *              integer represented as two's complement.
 460  */
 461
 462 extern size_t mp_octets2c(const mp */*m*/);
 463
 464 /* --- @mp_bits@ --- *
 465  *
 466  * Arguments:   @const mp *m@ = a multiprecision integer
 467  *
 468  * Returns:     The number of bits required to represent @m@.
 469  *
 470  * Use:         Calculates the external storage required for a multiprecision
 471  *              integer.
 472  */
 473
 474 extern unsigned long mp_bits(const mp */*m*/);
 475
 476 /* --- @mp_loadl@ --- *
 477  *
 478  * Arguments:   @mp *d@ = destination
 479  *              @const void *pv@ = pointer to source data
 480  *              @size_t sz@ = size of the source data
 481  *
 482  * Returns:     Resulting multiprecision number.
 483  *
 484  * Use:         Loads a multiprecision number from an array of octets.  The
 485  *              first byte in the array is the least significant.  More
 486  *              formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$%
 487  *              then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
 488  */
 489
 490 extern mp *mp_loadl(mp */*d*/, const void */*pv*/, size_t /*sz*/);
 491
 492 /* --- @mp_storel@ --- *
 493  *
 494  * Arguments:   @const mp *m@ = source
 495  *              @void *pv@ = pointer to output array
 496  *              @size_t sz@ = size of the output array
 497  *
 498  * Returns:     ---
 499  *
 500  * Use:         Stores a multiprecision number in an array of octets.  The
 501  *              first byte in the array is the least significant.  If the
 502  *              array is too small to represent the number, high-order bits
 503  *              are truncated; if the array is too large, high order bytes
 504  *              are filled with zeros.  More formally, if the number is
 505  *              %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
 506  *              then the array is %$b_0, b_1, \ldots, b_{n-1}$%.
 507  */
 508
 509 extern void mp_storel(const mp */*m*/, void */*pv*/, size_t /*sz*/);
 510
 511 /* --- @mp_loadb@ --- *
 512  *
 513  * Arguments:   @mp *d@ = destination
 514  *              @const void *pv@ = pointer to source data
 515  *              @size_t sz@ = size of the source data
 516  *
 517  * Returns:     Resulting multiprecision number.
 518  *
 519  * Use:         Loads a multiprecision number from an array of octets.  The
 520  *              last byte in the array is the least significant.  More
 521  *              formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$%
 522  *              then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%.
 523  */
 524
 525 extern mp *mp_loadb(mp */*d*/, const void */*pv*/, size_t /*sz*/);
 526
 527 /* --- @mp_storeb@ --- *
 528  *
 529  * Arguments:   @const mp *m@ = source
 530  *              @void *pv@ = pointer to output array
 531  *              @size_t sz@ = size of the output array
 532  *
 533  * Returns:     ---
 534  *
 535  * Use:         Stores a multiprecision number in an array of octets.  The
 536  *              last byte in the array is the least significant.  If the
 537  *              array is too small to represent the number, high-order bits
 538  *              are truncated; if the array is too large, high order bytes
 539  *              are filled with zeros.  More formally, if the number is
 540  *              %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%,
 541  *              then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%.
 542  */
 543
 544 extern void mp_storeb(const mp */*m*/, void */*pv*/, size_t /*sz*/);
 545
 546 /* --- @mp_loadl2c@ --- *
 547  *
 548  * Arguments:   @mp *d@ = destination
 549  *              @const void *pv@ = pointer to source data
 550  *              @size_t sz@ = size of the source data
 551  *
 552  * Returns:     Resulting multiprecision number.
 553  *
 554  * Use:         Loads a multiprecision number from an array of octets as
 555  *              two's complement.  The first byte in the array is the least
 556  *              significant.
 557  */
 558
 559 extern mp *mp_loadl2c(mp */*d*/, const void */*pv*/, size_t /*sz*/);
 560
 561 /* --- @mp_storel2c@ --- *
 562  *
 563  * Arguments:   @const mp *m@ = source
 564  *              @void *pv@ = pointer to output array
 565  *              @size_t sz@ = size of the output array
 566  *
 567  * Returns:     ---
 568  *
 569  * Use:         Stores a multiprecision number in an array of octets as two's
 570  *              complement.  The first byte in the array is the least
 571  *              significant.  If the array is too small to represent the
 572  *              number, high-order bits are truncated; if the array is too
 573  *              large, high order bytes are sign-extended.
 574  */
 575
 576 extern void mp_storel2c(const mp */*m*/, void */*pv*/, size_t /*sz*/);
 577
 578 /* --- @mp_loadb2c@ --- *
 579  *
 580  * Arguments:   @mp *d@ = destination
 581  *              @const void *pv@ = pointer to source data
 582  *              @size_t sz@ = size of the source data
 583  *
 584  * Returns:     Resulting multiprecision number.
 585  *
 586  * Use:         Loads a multiprecision number from an array of octets as
 587  *              two's complement.  The last byte in the array is the least
 588  *              significant.
 589  */
 590
 591 extern mp *mp_loadb2c(mp */*d*/, const void */*pv*/, size_t /*sz*/);
 592
 593 /* --- @mp_storeb2c@ --- *
 594  *
 595  * Arguments:   @const mp *m@ = source
 596  *              @void *pv@ = pointer to output array
 597  *              @size_t sz@ = size of the output array
 598  *
 599  * Returns:     ---
 600  *
 601  * Use:         Stores a multiprecision number in an array of octets, as
 602  *              two's complement.  The last byte in the array is the least
 603  *              significant.  If the array is too small to represent the
 604  *              number, high-order bits are truncated; if the array is too
 605  *              large, high order bytes are sign-extended.
 606  */
 607
 608 extern void mp_storeb2c(const mp */*m*/, void */*pv*/, size_t /*sz*/);
 609
 610 /*----- Bit operations ----------------------------------------------------*/
 611
 612 /* --- @mp_not@ --- *
 613  *
 614  * Arguments:   @mp *d@ = destination
 615  *              @mp *a@ = source
 616  *
 617  * Returns:     The bitwise complement of the source.
 618  */
 619
 620 extern mp *mp_not(mp */*d*/, mp */*a*/);
 621
 622 /* --- @mp_bitop@ --- *
 623  *
 624  * Arguments:   @mp *d@ = destination
 625  *              @mp *a, *b@ = sources
 626  *
 627  * Returns:     The result of the given bitwise operation.  These functions
 628  *              don't handle negative numbers at all sensibly.  For that, use
 629  *              the @...2c@ variants.  The functions are named after the
 630  *              truth tables they generate:
 631  *
 632  *                      a:      0011
 633  *                      b:      0101
 634  *                      @mpx_bitXXXX@
 635  */
 636
 637 #define MP_BITDECL(string)                                              \
 638   extern mp *mp_bit##string(mp */*d*/, mp */*a*/, mp */*b*/);
 639 MPX_DOBIN(MP_BITDECL)
 640
 641 /* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
 642  *
 643  * Synonyms for the commonly-used functions.
 644  */
 645
 646 #define mp_and  mp_bit0001
 647 #define mp_or   mp_bit0111
 648 #define mp_nand mp_bit1110
 649 #define mp_nor  mp_bit1000
 650 #define mp_xor  mp_bit0110
 651
 652 /* --- @mp_testbit@ --- *
 653  *
 654  * Arguments:   @mp *x@ = a large integer
 655  *              @unsigned long n@ = which bit to test
 656  *
 657  * Returns:     Nonzero if the bit is set, zero if not.
 658  */
 659
 660 extern int mp_testbit(mp */*x*/, unsigned long /*n*/);
 661
 662 /* --- @mp_setbit@, @mp_clearbit@ --- *
 663  *
 664  * Arguments:   @mp *d@ = a destination
 665  *              @mp *x@ = a large integer
 666  *              @unsigned long n@ = which bit to modify
 667  *
 668  * Returns:     The argument @x@, with the appropriate bit set or cleared.
 669  */
 670
 671 extern mp *mp_setbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);
 672 extern mp *mp_clearbit(mp */*d*/, mp */*x*/, unsigned long /*n*/);
 673
 674 /* --- @mp_lsl@, @mp_lslc@, @mp_lsr@ --- *
 675  *
 676  * Arguments:   @mp *d@ = destination
 677  *              @mp *a@ = source
 678  *              @size_t n@ = number of bits to move
 679  *
 680  * Returns:     Result, @a@ shifted left or right by @n@.
 681  *
 682  * Use:         Bitwise shift operators.  @mp_lslc@ fills the bits introduced
 683  *              on the right with ones instead of zeroes: it's used
 684  *              internally by @mp_lsl2c@, though it may be useful on its
 685  *              own.
 686  */
 687
 688 extern mp *mp_lsl(mp */*d*/, mp */*a*/, size_t /*n*/);
 689 extern mp *mp_lslc(mp */*d*/, mp */*a*/, size_t /*n*/);
 690 extern mp *mp_lsr(mp */*d*/, mp */*a*/, size_t /*n*/);
 691
 692 /* --- @mp_not2c@ --- *
 693  *
 694  * Arguments:   @mp *d@ = destination
 695  *              @mp *a@ = source
 696  *
 697  * Returns:     The sign-extended complement of the argument.
 698  */
 699
 700 extern mp *mp_not2c(mp */*d*/, mp */*a*/);
 701
 702 /* --- @mp_bitop2c@ --- *
 703  *
 704  * Arguments:   @mp *d@ = destination
 705  *              @mp *a, *b@ = sources
 706  *
 707  * Returns:     The result of the given bitwise operation.  Negative numbers
 708  *              are treated as two's complement, sign-extended infinitely to
 709  *              the left.  The functions are named after the truth tables
 710  *              they generate:
 711  *
 712  *                      a:      0011
 713  *                      b:      0101
 714  *                      @mpx_bitXXXX@
 715  */
 716
 717 #define MP_BIT2CDECL(string)                                            \
 718   extern mp *mp_bit##string##2c(mp */*d*/, mp */*a*/, mp */*b*/);
 719 MPX_DOBIN(MP_BIT2CDECL)
 720
 721 /* --- @mp_[n]and@, @mp_[n]or@, @mp_[n]xor@, @mp_not@ --- *
 722  *
 723  * Synonyms for the commonly-used functions.
 724  */
 725
 726 #define mp_and2c  mp_bit00012c
 727 #define mp_or2c   mp_bit01112c
 728 #define mp_nand2c mp_bit11102c
 729 #define mp_nor2c  mp_bit10002c
 730 #define mp_xor2c  mp_bit01102c
 731
 732 /* --- @mp_lsl2c@, @mp_lsr2c@ --- *
 733  *
 734  * Arguments:   @mp *d@ = destination
 735  *              @mp *a@ = source
 736  *              @size_t n@ = number of bits to move
 737  *
 738  * Returns:     Result, @a@ shifted left or right by @n@.  Handles the
 739  *              pretence of sign-extension for negative numbers.
 740  */
 741
 742 extern mp *mp_lsl2c(mp */*d*/, mp */*a*/, size_t /*n*/);
 743 extern mp *mp_lsr2c(mp */*d*/, mp */*a*/, size_t /*n*/);
 744
 745 /* --- @mp_testbit2c@ --- *
 746  *
 747  * Arguments:   @mp *x@ = a large integer
 748  *              @unsigned long n@ = which bit to test
 749  *
 750  * Returns:     Nonzero if the bit is set, zero if not.  Fakes up two's
 751  *              complement representation.
 752  */
 753
 754 extern int mp_testbit2c(mp */*x*/, unsigned long /*n*/);
 755
 756 /* --- @mp_setbit2c@, @mp_clearbit2c@ --- *
 757  *
 758  * Arguments:   @mp *d@ = a destination
 759  *              @mp *x@ = a large integer
 760  *              @unsigned long n@ = which bit to modify
 761  *
 762  * Returns:     The argument @x@, with the appropriate bit set or cleared.
 763  *              Fakes up two's complement representation.
 764  */
 765
 766 extern mp *mp_setbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);
 767 extern mp *mp_clearbit2c(mp */*d*/, mp */*x*/, unsigned long /*n*/);
 768
 769 /*----- Comparisons -------------------------------------------------------*/
 770
 771 /* --- @mp_eq@ --- *
 772  *
 773  * Arguments:   @const mp *a, *b@ = two numbers
 774  *
 775  * Returns:     Nonzero if the numbers are equal.
 776  */
 777
 778 extern int mp_eq(const mp */*a*/, const mp */*b*/);
 779
 780 #define MP_EQ(a, b)                                                     \
 781   ((((a)->f ^ (b)->f) & MP_NEG) == 0 &&                                 \
 782    mpx_ueq((a)->v, (a)->vl, (b)->v, (b)->vl))
 783
 784 /* --- @mp_cmp@ --- *
 785  *
 786  * Arguments:   @const mp *a, *b@ = two numbers
 787  *
 788  * Returns:     Less than, equal to or greater than zero, according to
 789  *              whether @a@ is less than, equal to or greater than @b@.
 790  */
 791
 792 extern int mp_cmp(const mp */*a*/, const mp */*b*/);
 793
 794 #define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0)
 795
 796 /* --- Other handy macros --- */
 797
 798 #define MP_NEGP(x) ((x)->f & MP_NEG)
 799 #define MP_ZEROP(x) (!MP_LEN(x))
 800 #define MP_POSP(x) (!MP_NEGP(x) && !MP_ZEROP(x))
 801 #define MP_ODDP(x) (!MP_ZEROP(x) && ((x)->v[0] & 1u))
 802 #define MP_EVENP(x) (!MP_ODDP(x))
 803
 804 /*----- Arithmetic operations ---------------------------------------------*/
 805
 806 /* --- @mp_neg@ --- *
 807  *
 808  * Arguments:   @mp *d@ = destination
 809  *              @mp *a@ = argument
 810  *
 811  * Returns:     The negation of the argument.
 812  *
 813  * Use:         Negates its argument.
 814  */
 815
 816 extern mp *mp_neg(mp */*d*/, mp */*a*/);
 817
 818 /* --- @mp_add@ --- *
 819  *
 820  * Arguments:   @mp *d@ = destination
 821  *              @mp *a, *b@ = sources
 822  *
 823  * Returns:     Result, @a@ added to @b@.
 824  */
 825
 826 extern mp *mp_add(mp */*d*/, mp */*a*/, mp */*b*/);
 827
 828 /* --- @mp_sub@ --- *
 829  *
 830  * Arguments:   @mp *d@ = destination
 831  *              @mp *a, *b@ = sources
 832  *
 833  * Returns:     Result, @b@ subtracted from @a@.
 834  */
 835
 836 extern mp *mp_sub(mp */*d*/, mp */*a*/, mp */*b*/);
 837
 838 /* --- @mp_mul@ --- *
 839  *
 840  * Arguments:   @mp *d@ = destination
 841  *              @mp *a, *b@ = sources
 842  *
 843  * Returns:     Result, @a@ multiplied by @b@.
 844  */
 845
 846 extern mp *mp_mul(mp */*d*/, mp */*a*/, mp */*b*/);
 847
 848 /* --- @mp_sqr@ --- *
 849  *
 850  * Arguments:   @mp *d@ = destination
 851  *              @mp *a@ = source
 852  *
 853  * Returns:     Result, @a@ squared.
 854  */
 855
 856 extern mp *mp_sqr(mp */*d*/, mp */*a*/);
 857
 858 /* --- @mp_div@ --- *
 859  *
 860  * Arguments:   @mp **qq, **rr@ = destination, quotient and remainder
 861  *              @mp *a, *b@ = sources
 862  *
 863  * Use:         Calculates the quotient and remainder when @a@ is divided by
 864  *              @b@.
 865  */
 866
 867 extern void mp_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/);
 868
 869 /* --- @mp_exp@ --- *
 870  *
 871  * Arguments:   @mp *d@ = fake destination
 872  *              @mp *a@ = base
 873  *              @mp *e@ = exponent
 874  *
 875  * Returns:     Result, %$a^e$%.
 876  */
 877
 878 extern mp *mp_exp(mp */*d*/, mp */*a*/, mp */*e*/);
 879
 880 /* --- @mp_odd@ --- *
 881  *
 882  * Arguments:   @mp *d@ = pointer to destination integer
 883  *              @mp *m@ = pointer to source integer
 884  *              @size_t *s@ = where to store the power of 2
 885  *
 886  * Returns:     An odd integer integer %$t$% such that %$m = 2^s t$%.
 887  *
 888  * Use:         Computes a power of two and an odd integer which, when
 889  *              multiplied, give a specified result.  This sort of thing is
 890  *              useful in number theory quite often.
 891  */
 892
 893 extern mp *mp_odd(mp */*d*/, mp */*m*/, size_t */*s*/);
 894
 895 /* --- @mp_leastcongruent@ --- *
 896  *
 897  * Arguments:   @mp *d@ = pointer to destination
 898  *              @mp *b@ = lower bound
 899  *              @mp *r@ = representative
 900  *              @mp *m@ = modulus
 901  *
 902  * Returns:     The smallest integer %$x \equiv r \pmod{m}$% such that
 903  *              %$x \ge b$%.
 904  */
 905
 906 extern mp *mp_leastcongruent(mp */*d*/, mp */*b*/, mp */*r*/, mp */*m*/);
 907
 908 /*----- More advanced algorithms ------------------------------------------*/
 909
 910 /* --- @mp_sqrt@ --- *
 911  *
 912  * Arguments:   @mp *d@ = pointer to destination integer
 913  *              @mp *a@ = (nonnegative) integer to take square root of
 914  *
 915  * Returns:     The largest integer %$x$% such that %$x^2 \le a$%.
 916  *
 917  * Use:         Computes integer square roots.
 918  *
 919  *              The current implementation isn't very good: it uses the
 920  *              Newton-Raphson method to find an approximation to %$a$%.  If
 921  *              there's any demand for a better version, I'll write one.
 922  */
 923
 924 extern mp *mp_sqrt(mp */*d*/, mp */*a*/);
 925
 926 /* --- @mp_nthrt@ --- *
 927  *
 928  * Arguments:   @mp *d@ = fake destination
 929  *              @mp *a@ = an integer
 930  *              @mp *n@ = a strictly positive integer
 931  *              @int *exectp_out@ = set nonzero if an exact solution is found
 932  *
 933  * Returns:     The integer %$\bigl\lfloor \sqrt[n]{a} \bigr\rfloor$%.
 934  *
 935  * Use:         Return an approximation to the %$n$%th root of %$a$%.
 936  *              Specifically, it returns the largest integer %$x$% such that
 937  *              %$x^n \le a$%.  If %$x^n = a$% then @*exactp_out@ is set
 938  *              nonzero; otherwise it is set zero.  (If @exactp_out@ is null
 939  *              then this information is discarded.)
 940  *
 941  *              The exponent %$n$% must be strictly positive: it's not clear
 942  *              to me what the right answer is for %$n \le 0$%.  If %$a$% is
 943  *              negative then %$n$% must be odd; otherwise there is no real
 944  *              solution.
 945  */
 946
 947 extern mp *mp_nthrt(mp */*d*/, mp */*a*/, mp */*n*/, int */*exactp_out*/);
 948
 949 /* --- @mp_perfect_power_p@ --- *
 950  *
 951  * Arguments:   @mp **x@ = where to write the base
 952  *              @mp **n@ = where to write the exponent
 953  *              @mp *a@ = an integer
 954  *
 955  * Returns:     Nonzero if %$a$% is a perfect power.
 956  *
 957  * Use:         Returns whether an integer %$a$% is a perfect power, i.e.,
 958  *              whether it can be written in the form %$a = x^n$% where
 959  *              %$|x| > 1$% and %$n > 1$% are integers.  If this is possible,
 960  *              then (a) store %$x$% and the largest such %$n$% in @*x@ and
 961  *              @*n@, and return nonzero; otherwise, store %$x = a$% and
 962  *              %$n = 1$% and return zero.  (Either @x@ or @n@, or both, may
 963  *              be null to discard these outputs.)
 964  *
 965  *              Note that %$-1$%, %$0$% and %$1$% are not considered perfect
 966  *              powers by this definition.  (The exponent is not well-defined
 967  *              in these cases, but it seemed better to implement a function
 968  *              which worked for all integers.)  Note also that %$-4$% is not
 969  *              a perfect power since it has no real square root.
 970  */
 971
 972 extern int mp_perfect_power_p(mp **/*x*/, mp **/*n*/, mp */*a*/);
 973
 974 /* --- @mp_gcd@ --- *
 975  *
 976  * Arguments:   @mp **gcd, **xx, **yy@ = where to write the results
 977  *              @mp *a, *b@ = sources (must be nonzero)
 978  *
 979  * Returns:     ---
 980  *
 981  * Use:         Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
 982  *              @ax + by = gcd(a, b)@.  This is useful for computing modular
 983  *              inverses.  Neither @a@ nor @b@ may be zero.
 984  */
 985
 986 extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
 987                    mp */*a*/, mp */*b*/);
 988
 989 /* -- @mp_modinv@ --- *
 990  *
 991  * Arguments:   @mp *d@ = destination
 992  *              @mp *x@ = argument
 993  *              @mp *p@ = modulus
 994  *
 995  * Returns:     The inverse %$x^{-1} \bmod p$%.
 996  *
 997  * Use:         Computes a modular inverse.    An assertion fails if %$p$%
 998  *              has no inverse.
 999  */
1000
1001 extern mp *mp_modinv(mp */*d*/, mp */*x*/, mp */*p*/);
1002
1003 /* --- @mp_jacobi@ --- *
1004  *
1005  * Arguments:   @mp *a@ = an integer
1006  *              @mp *n@ = another integer
1007  *
1008  * Returns:     @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
1009  *
1010  * Use:         Computes the Kronecker symbol %$\jacobi{a}{n}$%.  If @n@ is
1011  *              prime, this is the Legendre symbol and is equal to 1 if and
1012  *              only if @a@ is a quadratic residue mod @n@.  The result is
1013  *              zero if and only if @a@ and @n@ have a common factor greater
1014  *              than one.
1015  *
1016  *              If @n@ is composite, then this computes the Kronecker symbol
1017  *
1018  *                %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$%
1019  *
1020  *              where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the
1021  *              prime factorization of %$n$%.  The missing bits are:
1022  *
1023  *                * %$\jacobi{a}{1} = 1$%;
1024  *                * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if
1025  *                  positive;
1026  *                * %$\jacobi{a}{0} = 0$%;
1027  *                * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is
1028  *                  congruent to 1 or 7 (mod 8), or %$-1$% otherwise.
1029  *
1030  *              If %$n$% is positive and odd, then this is the Jacobi
1031  *              symbol.  (The Kronecker symbol is a consistent domain
1032  *              extension; the Jacobi symbol was implemented first, and the
1033  *              name stuck.)
1034  */
1035
1036 extern int mp_jacobi(mp */*a*/, mp */*n*/);
1037
1038 /* --- @mp_modsqrt@ --- *
1039  *
1040  * Arguments:   @mp *d@ = destination integer
1041  *              @mp *a@ = source integer
1042  *              @mp *p@ = modulus (must be prime)
1043  *
1044  * Returns:     If %$a$% is a quadratic residue, a square root of %$a$%; else
1045  *              a null pointer.
1046  *
1047  * Use:         Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
1048  *              if one exists; else a null pointer.  This function will not
1049  *              work if %$p$% is composite: you must factor the modulus, take
1050  *              a square root mod each factor, and recombine the results
1051  *              using the Chinese Remainder Theorem.
1052  *
1053  *              We guarantee that the square root returned is the smallest
1054  *              one (i.e., the `positive' square root).
1055  */
1056
1057 extern mp *mp_modsqrt(mp */*d*/, mp */*a*/, mp */*p*/);
1058
1059 /* --- @mp_modexp@ --- *
1060  *
1061  * Arguments:   @mp *d@ = fake destination
1062  *              @mp *x@ = base of exponentiation
1063  *              @mp *e@ = exponent
1064  *              @mp *n@ = modulus (must be positive)
1065  *
1066  * Returns:     The value %$x^e \bmod n$%.
1067  */
1068
1069 extern mp *mp_modexp(mp */*d*/, mp */*x*/, mp */*e*/, mp */*n*/);
1070
1071 /*----- Test harness support ----------------------------------------------*/
1072
1073 #include <mLib/testrig.h>
1074
1075 #ifndef CATACOMB_MPTEXT_H
1076 #  include "mptext.h"
1077 #endif
1078
1079 extern const test_type type_mp;
1080
1081 /*----- That's all, folks -------------------------------------------------*/
1082
1083 #ifdef __cplusplus
1084   }
1085 #endif
1086
1087 #endif