stdlib/gmp-impl.h

   1 /* Include file for internal GNU MP types and definitions.
   2
   3 Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Library General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Library General Public License
  18 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
  20
  21 #if ! defined (alloca)
  22 #if defined (__GNUC__)
  23 #define alloca __builtin_alloca
  24 #endif
  25 #endif
  26
  27 #if ! defined (alloca)
  28 #if defined (__sparc__) || defined (sparc) || defined (__sgi)
  29 #include <alloca.h>
  30 #endif
  31 #endif
  32
  33 #ifndef NULL
  34 #define NULL 0L
  35 #endif
  36
  37 #if ! defined (__GNUC__)
  38 #define inline                  /* Empty */
  39 void *alloca();
  40 #endif
  41
  42 #define ABS(x) (x >= 0 ? x : -x)
  43 #define MIN(l,o) ((l) < (o) ? (l) : (o))
  44 #define MAX(h,i) ((h) > (i) ? (h) : (i))
  45
  46 #include "gmp-mparam.h"
  47 /* #include "longlong.h" */
  48
  49 #ifdef __STDC__
  50 void *malloc (size_t);
  51 void *realloc (void *, size_t);
  52 void free (void *);
  53
  54 extern void *   (*_mp_allocate_func) (size_t);
  55 extern void *   (*_mp_reallocate_func) (void *, size_t, size_t);
  56 extern void     (*_mp_free_func) (void *, size_t);
  57
  58 void *_mp_default_allocate (size_t);
  59 void *_mp_default_reallocate (void *, size_t, size_t);
  60 void _mp_default_free (void *, size_t);
  61
  62 #else
  63
  64 #define const                   /* Empty */
  65 #define signed                  /* Empty */
  66
  67 void *malloc ();
  68 void *realloc ();
  69 void free ();
  70
  71 extern void *   (*_mp_allocate_func) ();
  72 extern void *   (*_mp_reallocate_func) ();
  73 extern void     (*_mp_free_func) ();
  74
  75 void *_mp_default_allocate ();
  76 void *_mp_default_reallocate ();
  77 void _mp_default_free ();
  78 #endif
  79
  80 /* Copy NLIMBS *limbs* from SRC to DST.  */
  81 #define MPN_COPY_INCR(DST, SRC, NLIMBS) \
  82   do {                                                                  \
  83     mp_size_t __i;                                                      \
  84     for (__i = 0; __i < (NLIMBS); __i++)                                \
  85       (DST)[__i] = (SRC)[__i];                                          \
  86   } while (0)
  87 #define MPN_COPY_DECR(DST, SRC, NLIMBS) \
  88   do {                                                                  \
  89     mp_size_t __i;                                                      \
  90     for (__i = (NLIMBS) - 1; __i >= 0; __i--)                           \
  91       (DST)[__i] = (SRC)[__i];                                          \
  92   } while (0)
  93 #define MPN_COPY MPN_COPY_INCR
  94
  95 /* Zero NLIMBS *limbs* AT DST.  */
  96 #define MPN_ZERO(DST, NLIMBS) \
  97   do {                                                                  \
  98     mp_size_t __i;                                                      \
  99     for (__i = 0; __i < (NLIMBS); __i++)                                \
 100       (DST)[__i] = 0;                                                   \
 101   } while (0)
 102
 103 #define MPN_NORMALIZE(DST, NLIMBS) \
 104   do {                                                                  \
 105     while (NLIMBS > 0)                                                  \
 106       {                                                                 \
 107         if ((DST)[(NLIMBS) - 1] != 0)                                   \
 108           break;                                                        \
 109         NLIMBS--;                                                       \
 110       }                                                                 \
 111   } while (0)
 112 #define MPN_NORMALIZE_NOT_ZERO(DST, NLIMBS) \
 113   do {                                                                  \
 114     while (1)                                                           \
 115       {                                                                 \
 116         if ((DST)[(NLIMBS) - 1] != 0)                                   \
 117           break;                                                        \
 118         NLIMBS--;                                                       \
 119       }                                                                 \
 120   } while (0)
 121
 122 /*  Swap (mp_ptr, mp_size_t) (U, UL) with (V, VL)  */
 123 #define MPN_SWAP(u, l, v, m) \
 124   do {                                                                  \
 125     { mp_ptr _; _ = (u), (u) = (v), (v) = _;}                           \
 126     { mp_size_t _; _ = (l), (l) = (m), (m) = _;}                        \
 127   } while (0)
 128
 129 /*  Return true iff the limb X has less bits than the limb Y.  */
 130 #define MPN_LESS_BITS_LIMB(x,y) ((x) < (y) && (x) < ((x) ^ (y)))
 131
 132 /*  Return true iff (mp_ptr, mp_size_t) (U, UL) has less bits than (V, VL).  */
 133 #define MPN_LESS_BITS(u, l, v, m) \
 134   ((l) < (m)                                                            \
 135    || ((l) == (m) && (l) != 0 && MPN_LESS_BITS_LIMB ((u)[(l - 1)], (v)[(l) - 1])))
 136
 137 /*  Return true iff (mp_ptr, mp_size_t) (U, UL) has more bits than (V, VL).  */
 138 #define MPN_MORE_BITS(u, l, v, m) MPN_LESS_BITS (v, m, u, l)
 139
 140 /*  Perform twos complement on (mp_ptr, mp_size_t) (U, UL),
 141     putting result at (v, VL).  Precondition: U[0] != 0.  */
 142 #define MPN_COMPL_INCR(u, v, l) \
 143   do {                                                                  \
 144     mp_size_t _ = 0;                                                    \
 145     (u)[0] = -(v)[_];                                                   \
 146     while (_++ < (l))                                                   \
 147       (u)[_] = ~(v)[_];                                                 \
 148   } while (0)
 149 #define MPN_COMPL MPN_COMPL_INCR
 150
 151 /* Initialize the MP_INT X with space for NLIMBS limbs.
 152    X should be a temporary variable, and it will be automatically
 153    cleared out when the running function returns.
 154    We use __x here to make it possible to accept both mpz_ptr and mpz_t
 155    arguments.  */
 156 #define MPZ_TMP_INIT(X, NLIMBS) \
 157   do {                                                                  \
 158     mpz_ptr __x = (X);                                                  \
 159     __x->alloc = (NLIMBS);                                              \
 160     __x->d = (mp_ptr) alloca ((NLIMBS) * BYTES_PER_MP_LIMB);            \
 161   } while (0)
 162
 163 #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
 164   do {                                                                  \
 165     if ((size) < KARATSUBA_THRESHOLD)                                   \
 166       ____mpn_mul_n_basecase (prodp, up, vp, size);                     \
 167     else                                                                \
 168       ____mpn_mul_n (prodp, up, vp, size, tspace);                      \
 169   } while (0);
 170 #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
 171   do {                                                                  \
 172     if ((size) < KARATSUBA_THRESHOLD)                                   \
 173       ____mpn_sqr_n_basecase (prodp, up, size);                         \
 174     else                                                                \
 175       ____mpn_sqr_n (prodp, up, size, tspace);                          \
 176   } while (0);
 177
 178 /* Structure for conversion between internal binary format and
 179    strings in base 2..36.  */
 180 struct bases
 181 {
 182   /* Number of digits in the conversion base that always fits in
 183      an mp_limb.  For example, for base 10 this is 10, since
 184      2**32 = 4294967296 has ten digits.  */
 185   int chars_per_limb;
 186
 187   /* log(2)/log(conversion_base) */
 188   float chars_per_bit_exactly;
 189
 190   /* big_base is conversion_base**chars_per_limb, i.e. the biggest
 191      number that fits a word, built by factors of conversion_base.
 192      Exception: For 2, 4, 8, etc, big_base is log2(base), i.e. the
 193      number of bits used to represent each digit in the base.  */
 194   mp_limb big_base;
 195
 196   /* big_base_inverted is a BITS_PER_MP_LIMB bit approximation to
 197      1/big_base, represented as a fixed-point number.  Instead of
 198      dividing by big_base an application can choose to multiply
 199      by big_base_inverted.  */
 200   mp_limb big_base_inverted;
 201 };
 202
 203 extern const struct bases __mp_bases[];
 204 extern mp_size_t __gmp_default_fp_limb_precision;
 205
 206 /* Divide the two-limb number in (NH,,NL) by D, with DI being the largest
 207    limb not larger than (2**(2*BITS_PER_MP_LIMB))/D - (2**BITS_PER_MP_LIMB).
 208    If this would yield overflow, DI should be the largest possible number
 209    (i.e., only ones).  For correct operation, the most significant bit of D
 210    has to be set.  Put the quotient in Q and the remainder in R.  */
 211 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
 212   do {                                                                  \
 213     mp_limb _q, _ql, _r;                                                \
 214     mp_limb _xh, _xl;                                                   \
 215     umul_ppmm (_q, _ql, (nh), (di));                                    \
 216     _q += (nh);                 /* DI is 2**BITS_PER_MP_LIMB too small */\
 217     umul_ppmm (_xh, _xl, _q, (d));                                      \
 218     sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl);                         \
 219     if (_xh != 0)                                                       \
 220       {                                                                 \
 221         sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                          \
 222         _q += 1;                                                        \
 223         if (_xh != 0)                                                   \
 224           {                                                             \
 225             sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                      \
 226             _q += 1;                                                    \
 227           }                                                             \
 228       }                                                                 \
 229     if (_r >= (d))                                                      \
 230       {                                                                 \
 231         _r -= (d);                                                      \
 232         _q += 1;                                                        \
 233       }                                                                 \
 234     (r) = _r;                                                           \
 235     (q) = _q;                                                           \
 236   } while (0)
 237 /* Like udiv_qrnnd_preinv, but for for any value D.  DNORM is D shifted left
 238    so that its most significant bit is set.  LGUP is ceil(log2(D)).  */
 239 #define udiv_qrnnd_preinv2gen(q, r, nh, nl, d, di, dnorm, lgup) \
 240   do {                                                                  \
 241     mp_limb n2, n10, n1, nadj, q1;                                      \
 242     mp_limb _xh, _xl;                                                   \
 243     n2 = ((nh) << (BITS_PER_MP_LIMB - (lgup))) + ((nl) >> 1 >> (l - 1));\
 244     n10 = (nl) << (BITS_PER_MP_LIMB - (lgup));                          \
 245     n1 = ((mp_limb_signed) n10 >> (BITS_PER_MP_LIMB - 1));              \
 246     nadj = n10 + (n1 & (dnorm));                                        \
 247     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 248     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 249     q1 = ~(n2 + _xh);                                                   \
 250     umul_ppmm (_xh, _xl, q1, d);                                        \
 251     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 252     _xh -= (d);                                                         \
 253     (r) = _xl + ((d) & _xh);                                            \
 254     (q) = _xh - q1;                                                     \
 255   } while (0)
 256 /* Exactly like udiv_qrnnd_preinv, but branch-free.  It is not clear which
 257    version to use.  */
 258 #define udiv_qrnnd_preinv2norm(q, r, nh, nl, d, di) \
 259   do {                                                                  \
 260     mp_limb n2, n10, n1, nadj, q1;                                      \
 261     mp_limb _xh, _xl;                                                   \
 262     n2 = (nh);                                                          \
 263     n10 = (nl);                                                         \
 264     n1 = ((mp_limb_signed) n10 >> (BITS_PER_MP_LIMB - 1));              \
 265     nadj = n10 + (n1 & (d));                                            \
 266     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 267     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 268     q1 = ~(n2 + _xh);                                                   \
 269     umul_ppmm (_xh, _xl, q1, d);                                        \
 270     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 271     _xh -= (d);                                                         \
 272     (r) = _xl + ((d) & _xh);                                            \
 273     (q) = _xh - q1;                                                     \
 274   } while (0)
 275
 276 #if defined (__GNUC__)
 277 /* Define stuff for longlong.h.  */
 278 typedef unsigned int UQItype    __attribute__ ((mode (QI)));
 279 typedef          int SItype     __attribute__ ((mode (SI)));
 280 typedef unsigned int USItype    __attribute__ ((mode (SI)));
 281 typedef          int DItype     __attribute__ ((mode (DI)));
 282 typedef unsigned int UDItype    __attribute__ ((mode (DI)));
 283 #else
 284 typedef unsigned char UQItype;
 285 typedef          long SItype;
 286 typedef unsigned long USItype;
 287 #endif
 288
 289 typedef mp_limb UWtype;
 290 typedef unsigned int UHWtype;
 291 #define W_TYPE_SIZE BITS_PER_MP_LIMB
 292
 293
 294 #ifndef IEEE_DOUBLE_BIG_ENDIAN
 295 #define IEEE_DOUBLE_BIG_ENDIAN 1
 296 #endif
 297
 298 #if IEEE_DOUBLE_BIG_ENDIAN
 299 union ieee_double_extract
 300 {
 301   struct
 302     {
 303       unsigned long sig:1;
 304       unsigned long exp:11;
 305       unsigned long manh:20;
 306       unsigned long manl:32;
 307     } s;
 308   double d;
 309 };
 310 #else
 311 union ieee_double_extract
 312 {
 313   struct
 314     {
 315       unsigned long manl:32;
 316       unsigned long manh:20;
 317       unsigned long exp:11;
 318       unsigned long sig:1;
 319     } s;
 320   double d;
 321 };
 322 #endif