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 an mp_limb.
 183      For example, for base 10 on a machine where a mp_limb has 32 bits this
 184      is 9, since 10**9 is the largest number that fits into a mp_limb.  */
 185   int chars_per_limb;
 186
 187   /* log(2)/log(conversion_base) */
 188   float chars_per_bit_exactly;
 189
 190   /* base**chars_per_limb, i.e. the biggest number that fits a word, built by
 191      factors of base.  Exception: For 2, 4, 8, etc, big_base is log2(base),
 192      i.e. the number of bits used to represent each digit in the base.  */
 193   mp_limb big_base;
 194
 195   /* A BITS_PER_MP_LIMB bit approximation to 1/big_base, represented as a
 196      fixed-point number.  Instead of dividing by big_base an application can
 197      choose to multiply by big_base_inverted.  */
 198   mp_limb big_base_inverted;
 199 };
 200
 201 extern const struct bases __mp_bases[];
 202 extern mp_size_t __gmp_default_fp_limb_precision;
 203
 204 /* Divide the two-limb number in (NH,,NL) by D, with DI being the largest
 205    limb not larger than (2**(2*BITS_PER_MP_LIMB))/D - (2**BITS_PER_MP_LIMB).
 206    If this would yield overflow, DI should be the largest possible number
 207    (i.e., only ones).  For correct operation, the most significant bit of D
 208    has to be set.  Put the quotient in Q and the remainder in R.  */
 209 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
 210   do {                                                                  \
 211     mp_limb _q, _ql, _r;                                                \
 212     mp_limb _xh, _xl;                                                   \
 213     umul_ppmm (_q, _ql, (nh), (di));                                    \
 214     _q += (nh);                 /* DI is 2**BITS_PER_MP_LIMB too small */\
 215     umul_ppmm (_xh, _xl, _q, (d));                                      \
 216     sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl);                         \
 217     if (_xh != 0)                                                       \
 218       {                                                                 \
 219         sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                          \
 220         _q += 1;                                                        \
 221         if (_xh != 0)                                                   \
 222           {                                                             \
 223             sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                      \
 224             _q += 1;                                                    \
 225           }                                                             \
 226       }                                                                 \
 227     if (_r >= (d))                                                      \
 228       {                                                                 \
 229         _r -= (d);                                                      \
 230         _q += 1;                                                        \
 231       }                                                                 \
 232     (r) = _r;                                                           \
 233     (q) = _q;                                                           \
 234   } while (0)
 235 /* Like udiv_qrnnd_preinv, but for for any value D.  DNORM is D shifted left
 236    so that its most significant bit is set.  LGUP is ceil(log2(D)).  */
 237 #define udiv_qrnnd_preinv2gen(q, r, nh, nl, d, di, dnorm, lgup) \
 238   do {                                                                  \
 239     mp_limb n2, n10, n1, nadj, q1;                                      \
 240     mp_limb _xh, _xl;                                                   \
 241     n2 = ((nh) << (BITS_PER_MP_LIMB - (lgup))) + ((nl) >> 1 >> (l - 1));\
 242     n10 = (nl) << (BITS_PER_MP_LIMB - (lgup));                          \
 243     n1 = ((mp_limb_signed) n10 >> (BITS_PER_MP_LIMB - 1));              \
 244     nadj = n10 + (n1 & (dnorm));                                        \
 245     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 246     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 247     q1 = ~(n2 + _xh);                                                   \
 248     umul_ppmm (_xh, _xl, q1, d);                                        \
 249     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 250     _xh -= (d);                                                         \
 251     (r) = _xl + ((d) & _xh);                                            \
 252     (q) = _xh - q1;                                                     \
 253   } while (0)
 254 /* Exactly like udiv_qrnnd_preinv, but branch-free.  It is not clear which
 255    version to use.  */
 256 #define udiv_qrnnd_preinv2norm(q, r, nh, nl, d, di) \
 257   do {                                                                  \
 258     mp_limb n2, n10, n1, nadj, q1;                                      \
 259     mp_limb _xh, _xl;                                                   \
 260     n2 = (nh);                                                          \
 261     n10 = (nl);                                                         \
 262     n1 = ((mp_limb_signed) n10 >> (BITS_PER_MP_LIMB - 1));              \
 263     nadj = n10 + (n1 & (d));                                            \
 264     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 265     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 266     q1 = ~(n2 + _xh);                                                   \
 267     umul_ppmm (_xh, _xl, q1, d);                                        \
 268     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 269     _xh -= (d);                                                         \
 270     (r) = _xl + ((d) & _xh);                                            \
 271     (q) = _xh - q1;                                                     \
 272   } while (0)
 273
 274 #if defined (__GNUC__)
 275 /* Define stuff for longlong.h.  */
 276 typedef unsigned int UQItype    __attribute__ ((mode (QI)));
 277 typedef          int SItype     __attribute__ ((mode (SI)));
 278 typedef unsigned int USItype    __attribute__ ((mode (SI)));
 279 typedef          int DItype     __attribute__ ((mode (DI)));
 280 typedef unsigned int UDItype    __attribute__ ((mode (DI)));
 281 #else
 282 typedef unsigned char UQItype;
 283 typedef          long SItype;
 284 typedef unsigned long USItype;
 285 #endif
 286
 287 typedef mp_limb UWtype;
 288 typedef unsigned int UHWtype;
 289 #define W_TYPE_SIZE BITS_PER_MP_LIMB
 290
 291
 292 #ifndef IEEE_DOUBLE_BIG_ENDIAN
 293 #define IEEE_DOUBLE_BIG_ENDIAN 1
 294 #endif
 295
 296 #if IEEE_DOUBLE_BIG_ENDIAN
 297 union ieee_double_extract
 298 {
 299   struct
 300     {
 301       unsigned long sig:1;
 302       unsigned long exp:11;
 303       unsigned long manh:20;
 304       unsigned long manl:32;
 305     } s;
 306   double d;
 307 };
 308 #else
 309 union ieee_double_extract
 310 {
 311   struct
 312     {
 313       unsigned long manl:32;
 314       unsigned long manh:20;
 315       unsigned long exp:11;
 316       unsigned long sig:1;
 317     } s;
 318   double d;
 319 };
 320 #endif