stdlib/gmp-impl.h

   1 /* Include file for internal GNU MP types and definitions.
   2
   3 Copyright (C) 1991, 1993, 1994, 1995, 1996 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., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 /* When using gcc, make sure to use its builtin alloca.  */
  23 #if ! defined (alloca) && defined (__GNUC__)
  24 #define alloca __builtin_alloca
  25 #define HAVE_ALLOCA
  26 #endif
  27
  28 /* When using cc, do whatever necessary to allow use of alloca.  For many
  29    machines, this means including alloca.h.  IBM's compilers need a #pragma
  30    in "each module that needs to use alloca".  */
  31 #if ! defined (alloca)
  32 /* We need lots of variants for MIPS, to cover all versions and perversions
  33    of OSes for MIPS.  */
  34 #if defined (__mips) || defined (MIPSEL) || defined (MIPSEB) \
  35  || defined (_MIPSEL) || defined (_MIPSEB) || defined (__sgi) \
  36  || defined (__alpha) || defined (__sparc) || defined (sparc) \
  37  || defined (__ksr__)
  38 #include <alloca.h>
  39 #define HAVE_ALLOCA
  40 #endif
  41 #if defined (_IBMR2)
  42 #pragma alloca
  43 #define HAVE_ALLOCA
  44 #endif
  45 #if defined (__DECC)
  46 #define alloca(x) __ALLOCA(x)
  47 #define HAVE_ALLOCA
  48 #endif
  49 #endif
  50
  51 #if ! defined (HAVE_ALLOCA) || USE_STACK_ALLOC
  52 #include "stack-alloc.h"
  53 #else
  54 #define TMP_DECL(m)
  55 #define TMP_ALLOC(x) alloca(x)
  56 #define TMP_MARK(m)
  57 #define TMP_FREE(m)
  58 #endif
  59
  60 #ifndef NULL
  61 #define NULL ((void *) 0)
  62 #endif
  63
  64 #if ! defined (__GNUC__)
  65 #define inline                  /* Empty */
  66 #endif
  67
  68 #define ABS(x) (x >= 0 ? x : -x)
  69 #define MIN(l,o) ((l) < (o) ? (l) : (o))
  70 #define MAX(h,i) ((h) > (i) ? (h) : (i))
  71
  72 /* Field access macros.  */
  73 #define SIZ(x) ((x)->_mp_size)
  74 #define ABSIZ(x) ABS (SIZ (x))
  75 #define PTR(x) ((x)->_mp_d)
  76 #define EXP(x) ((x)->_mp_exp)
  77 #define PREC(x) ((x)->_mp_prec)
  78 #define ALLOC(x) ((x)->_mp_alloc)
  79
  80 #include "gmp-mparam.h"
  81 /* #include "longlong.h" */
  82
  83 #if defined (__STDC__)  || defined (__cplusplus)
  84 void *malloc (size_t);
  85 void *realloc (void *, size_t);
  86 void free (void *);
  87
  88 extern void *   (*_mp_allocate_func) (size_t);
  89 extern void *   (*_mp_reallocate_func) (void *, size_t, size_t);
  90 extern void     (*_mp_free_func) (void *, size_t);
  91
  92 void *_mp_default_allocate (size_t);
  93 void *_mp_default_reallocate (void *, size_t, size_t);
  94 void _mp_default_free (void *, size_t);
  95
  96 #else
  97
  98 #define const                   /* Empty */
  99 #define signed                  /* Empty */
 100
 101 void *malloc ();
 102 void *realloc ();
 103 void free ();
 104
 105 extern void *   (*_mp_allocate_func) ();
 106 extern void *   (*_mp_reallocate_func) ();
 107 extern void     (*_mp_free_func) ();
 108
 109 void *_mp_default_allocate ();
 110 void *_mp_default_reallocate ();
 111 void _mp_default_free ();
 112 #endif
 113
 114 /* Copy NLIMBS *limbs* from SRC to DST.  */
 115 #define MPN_COPY_INCR(DST, SRC, NLIMBS) \
 116   do {                                                                  \
 117     mp_size_t __i;                                                      \
 118     for (__i = 0; __i < (NLIMBS); __i++)                                \
 119       (DST)[__i] = (SRC)[__i];                                          \
 120   } while (0)
 121 #define MPN_COPY_DECR(DST, SRC, NLIMBS) \
 122   do {                                                                  \
 123     mp_size_t __i;                                                      \
 124     for (__i = (NLIMBS) - 1; __i >= 0; __i--)                           \
 125       (DST)[__i] = (SRC)[__i];                                          \
 126   } while (0)
 127 #define MPN_COPY MPN_COPY_INCR
 128
 129 /* Zero NLIMBS *limbs* AT DST.  */
 130 #define MPN_ZERO(DST, NLIMBS) \
 131   do {                                                                  \
 132     mp_size_t __i;                                                      \
 133     for (__i = 0; __i < (NLIMBS); __i++)                                \
 134       (DST)[__i] = 0;                                                   \
 135   } while (0)
 136
 137 #define MPN_NORMALIZE(DST, NLIMBS) \
 138   do {                                                                  \
 139     while (NLIMBS > 0)                                                  \
 140       {                                                                 \
 141         if ((DST)[(NLIMBS) - 1] != 0)                                   \
 142           break;                                                        \
 143         NLIMBS--;                                                       \
 144       }                                                                 \
 145   } while (0)
 146 #define MPN_NORMALIZE_NOT_ZERO(DST, NLIMBS) \
 147   do {                                                                  \
 148     while (1)                                                           \
 149       {                                                                 \
 150         if ((DST)[(NLIMBS) - 1] != 0)                                   \
 151           break;                                                        \
 152         NLIMBS--;                                                       \
 153       }                                                                 \
 154   } while (0)
 155
 156 /* Initialize the MP_INT X with space for NLIMBS limbs.
 157    X should be a temporary variable, and it will be automatically
 158    cleared out when the running function returns.
 159    We use __x here to make it possible to accept both mpz_ptr and mpz_t
 160    arguments.  */
 161 #define MPZ_TMP_INIT(X, NLIMBS) \
 162   do {                                                                  \
 163     mpz_ptr __x = (X);                                                  \
 164     __x->_mp_alloc = (NLIMBS);                                          \
 165     __x->_mp_d = (mp_ptr) TMP_ALLOC ((NLIMBS) * BYTES_PER_MP_LIMB);     \
 166   } while (0)
 167
 168 #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
 169   do {                                                                  \
 170     if ((size) < KARATSUBA_THRESHOLD)                                   \
 171       impn_mul_n_basecase (prodp, up, vp, size);                        \
 172     else                                                                \
 173       impn_mul_n (prodp, up, vp, size, tspace);                 \
 174   } while (0);
 175 #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
 176   do {                                                                  \
 177     if ((size) < KARATSUBA_THRESHOLD)                                   \
 178       impn_sqr_n_basecase (prodp, up, size);                            \
 179     else                                                                \
 180       impn_sqr_n (prodp, up, size, tspace);                             \
 181   } while (0);
 182
 183 /* Structure for conversion between internal binary format and
 184    strings in base 2..36.  */
 185 struct bases
 186 {
 187   /* Number of digits in the conversion base that always fits in an mp_limb_t.
 188      For example, for base 10 on a machine where a mp_limb_t has 32 bits this
 189      is 9, since 10**9 is the largest number that fits into a mp_limb_t.  */
 190   int chars_per_limb;
 191
 192   /* log(2)/log(conversion_base) */
 193   float chars_per_bit_exactly;
 194
 195   /* base**chars_per_limb, i.e. the biggest number that fits a word, built by
 196      factors of base.  Exception: For 2, 4, 8, etc, big_base is log2(base),
 197      i.e. the number of bits used to represent each digit in the base.  */
 198   mp_limb_t big_base;
 199
 200   /* A BITS_PER_MP_LIMB bit approximation to 1/big_base, represented as a
 201      fixed-point number.  Instead of dividing by big_base an application can
 202      choose to multiply by big_base_inverted.  */
 203   mp_limb_t big_base_inverted;
 204 };
 205
 206 extern const struct bases __mp_bases[];
 207 extern mp_size_t __gmp_default_fp_limb_precision;
 208
 209 /* Divide the two-limb number in (NH,,NL) by D, with DI being the largest
 210    limb not larger than (2**(2*BITS_PER_MP_LIMB))/D - (2**BITS_PER_MP_LIMB).
 211    If this would yield overflow, DI should be the largest possible number
 212    (i.e., only ones).  For correct operation, the most significant bit of D
 213    has to be set.  Put the quotient in Q and the remainder in R.  */
 214 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
 215   do {                                                                  \
 216     mp_limb_t _q, _ql, _r;                                              \
 217     mp_limb_t _xh, _xl;                                                 \
 218     umul_ppmm (_q, _ql, (nh), (di));                                    \
 219     _q += (nh);                 /* DI is 2**BITS_PER_MP_LIMB too small */\
 220     umul_ppmm (_xh, _xl, _q, (d));                                      \
 221     sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl);                         \
 222     if (_xh != 0)                                                       \
 223       {                                                                 \
 224         sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                          \
 225         _q += 1;                                                        \
 226         if (_xh != 0)                                                   \
 227           {                                                             \
 228             sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                      \
 229             _q += 1;                                                    \
 230           }                                                             \
 231       }                                                                 \
 232     if (_r >= (d))                                                      \
 233       {                                                                 \
 234         _r -= (d);                                                      \
 235         _q += 1;                                                        \
 236       }                                                                 \
 237     (r) = _r;                                                           \
 238     (q) = _q;                                                           \
 239   } while (0)
 240 /* Like udiv_qrnnd_preinv, but for for any value D.  DNORM is D shifted left
 241    so that its most significant bit is set.  LGUP is ceil(log2(D)).  */
 242 #define udiv_qrnnd_preinv2gen(q, r, nh, nl, d, di, dnorm, lgup) \
 243   do {                                                                  \
 244     mp_limb_t n2, n10, n1, nadj, q1;                                    \
 245     mp_limb_t _xh, _xl;                                                 \
 246     n2 = ((nh) << (BITS_PER_MP_LIMB - (lgup))) + ((nl) >> 1 >> (l - 1));\
 247     n10 = (nl) << (BITS_PER_MP_LIMB - (lgup));                          \
 248     n1 = ((mp_limb_signed_t) n10 >> (BITS_PER_MP_LIMB - 1));            \
 249     nadj = n10 + (n1 & (dnorm));                                        \
 250     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 251     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 252     q1 = ~(n2 + _xh);                                                   \
 253     umul_ppmm (_xh, _xl, q1, d);                                        \
 254     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 255     _xh -= (d);                                                         \
 256     (r) = _xl + ((d) & _xh);                                            \
 257     (q) = _xh - q1;                                                     \
 258   } while (0)
 259 /* Exactly like udiv_qrnnd_preinv, but branch-free.  It is not clear which
 260    version to use.  */
 261 #define udiv_qrnnd_preinv2norm(q, r, nh, nl, d, di) \
 262   do {                                                                  \
 263     mp_limb_t n2, n10, n1, nadj, q1;                                    \
 264     mp_limb_t _xh, _xl;                                                 \
 265     n2 = (nh);                                                          \
 266     n10 = (nl);                                                         \
 267     n1 = ((mp_limb_signed_t) n10 >> (BITS_PER_MP_LIMB - 1));            \
 268     nadj = n10 + (n1 & (d));                                            \
 269     umul_ppmm (_xh, _xl, di, n2 - n1);                                  \
 270     add_ssaaaa (_xh, _xl, _xh, _xl, 0, nadj);                           \
 271     q1 = ~(n2 + _xh);                                                   \
 272     umul_ppmm (_xh, _xl, q1, d);                                        \
 273     add_ssaaaa (_xh, _xl, _xh, _xl, nh, nl);                            \
 274     _xh -= (d);                                                         \
 275     (r) = _xl + ((d) & _xh);                                            \
 276     (q) = _xh - q1;                                                     \
 277   } while (0)
 278
 279 #if defined (__GNUC__)
 280 /* Define stuff for longlong.h.  */
 281 typedef unsigned int UQItype    __attribute__ ((mode (QI)));
 282 typedef          int SItype     __attribute__ ((mode (SI)));
 283 typedef unsigned int USItype    __attribute__ ((mode (SI)));
 284 typedef          int DItype     __attribute__ ((mode (DI)));
 285 typedef unsigned int UDItype    __attribute__ ((mode (DI)));
 286 #else
 287 typedef unsigned char UQItype;
 288 typedef          long SItype;
 289 typedef unsigned long USItype;
 290 #endif
 291
 292 typedef mp_limb_t UWtype;
 293 typedef unsigned int UHWtype;
 294 #define W_TYPE_SIZE BITS_PER_MP_LIMB
 295
 296 /* Internal mpn calls */
 297 #define impn_mul_n_basecase     __MPN(impn_mul_n_basecase)
 298 #define impn_mul_n              __MPN(impn_mul_n)
 299 #define impn_sqr_n_basecase     __MPN(impn_sqr_n_basecase)
 300 #define impn_sqr_n              __MPN(impn_sqr_n)
 301
 302 #ifndef IEEE_DOUBLE_BIG_ENDIAN
 303 #define IEEE_DOUBLE_BIG_ENDIAN 1
 304 #endif
 305
 306 #if IEEE_DOUBLE_BIG_ENDIAN
 307 union ieee_double_extract
 308 {
 309   struct
 310     {
 311       unsigned int sig:1;
 312       unsigned int exp:11;
 313       unsigned int manh:20;
 314       unsigned int manl:32;
 315     } s;
 316   double d;
 317 };
 318 #else
 319 union ieee_double_extract
 320 {
 321   struct
 322     {
 323       unsigned int manl:32;
 324       unsigned int manh:20;
 325       unsigned int exp:11;
 326       unsigned int sig:1;
 327     } s;
 328   double d;
 329 };
 330 #endif