stdlib/gmp-impl.h

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