newlib/libc/stdlib/mprec.h

   1 /****************************************************************
   2  *
   3  * The author of this software is David M. Gay.
   4  *
   5  * Copyright (c) 1991 by AT&T.
   6  *
   7  * Permission to use, copy, modify, and distribute this software for any
   8  * purpose without fee is hereby granted, provided that this entire notice
   9  * is included in all copies of any software which is or includes a copy
  10  * or modification of this software and in all copies of the supporting
  11  * documentation for such software.
  12  *
  13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17  *
  18  ***************************************************************/
  19
  20 /* Please send bug reports to
  21         David M. Gay
  22         AT&T Bell Laboratories, Room 2C-463
  23         600 Mountain Avenue
  24         Murray Hill, NJ 07974-2070
  25         U.S.A.
  26         dmg@research.att.com or research!dmg
  27  */
  28
  29 #include <machine/ieeefp.h>
  30 #include <math.h>
  31 #include <float.h>
  32 #include <errno.h>
  33 #include <assert.h>
  34 #include <sys/config.h>
  35 #include <sys/types.h>
  36 #include "../locale/setlocale.h"
  37
  38 #ifdef __IEEE_LITTLE_ENDIAN
  39 #define IEEE_8087
  40 #endif
  41
  42 #ifdef __IEEE_BIG_ENDIAN
  43 #define IEEE_MC68k
  44 #endif
  45
  46 #ifdef __Z8000__
  47 #define Just_16
  48 #endif
  49
  50 #ifdef DEBUG
  51 #include "stdio.h"
  52 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  53 #endif
  54
  55 #ifdef Unsigned_Shifts
  56 #define Sign_Extend(a,b) if (b < 0) a |= (__uint32_t)0xffff0000;
  57 #else
  58 #define Sign_Extend(a,b) /*no-op*/
  59 #endif
  60
  61 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  62 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  63 #endif
  64
  65 /* If we are going to examine or modify specific bits in a double using
  66    the word0 and/or word1 macros, then we must wrap the double inside
  67    a union.  This is necessary to avoid undefined behavior according to
  68    the ANSI C spec.  */
  69 union double_union
  70 {
  71   double d;
  72   __uint32_t i[2];
  73 };
  74
  75 #ifdef IEEE_8087
  76 #define word0(x) (x.i[1])
  77 #define word1(x) (x.i[0])
  78 #else
  79 #define word0(x) (x.i[0])
  80 #define word1(x) (x.i[1])
  81 #endif
  82
  83 /* The following is taken from gdtoaimp.h for use with new strtod, but
  84    adjusted to avoid invalid type-punning.  */
  85 typedef __int32_t Long;
  86
  87 /* Unfortunately, because __ULong might be a different type than
  88    __uint32_t, we can't re-use union double_union as-is without
  89    further edits in strtod.c.  */
  90 typedef union { double d; __ULong i[2]; } U;
  91
  92 #define dword0(x) word0(x)
  93 #define dword1(x) word1(x)
  94 #define dval(x) (x.d)
  95
  96 #undef SI
  97 #ifdef Sudden_Underflow
  98 #define SI 1
  99 #else
 100 #define SI 0
 101 #endif
 102
 103 #define Storeinc(a,b,c) (*(a)++ = ((b) << 16) | ((c) & 0xffff))
 104
 105 /* #define P DBL_MANT_DIG */
 106 /* Ten_pmax = floor(P*log(2)/log(5)) */
 107 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 108 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 109 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
 110
 111 #if defined(IEEE_8087) + defined(IEEE_MC68k)
 112 #if defined (_DOUBLE_IS_32BITS)
 113 #define Exp_shift   23
 114 #define Exp_shift1  23
 115 #define Exp_msk1    ((__uint32_t)0x00800000L)
 116 #define Exp_msk11   ((__uint32_t)0x00800000L)
 117 #define Exp_mask    ((__uint32_t)0x7f800000L)
 118 #define P           24
 119 #define Bias        127
 120 #define NO_HEX_FP   /* not supported in this case */
 121 #define IEEE_Arith
 122 #define Emin        (-126)
 123 #define Exp_1       ((__uint32_t)0x3f800000L)
 124 #define Exp_11      ((__uint32_t)0x3f800000L)
 125 #define Ebits       8
 126 #define Frac_mask   ((__uint32_t)0x007fffffL)
 127 #define Frac_mask1  ((__uint32_t)0x007fffffL)
 128 #define Ten_pmax    10
 129 #define Sign_bit    ((__uint32_t)0x80000000L)
 130 #define Ten_pmax    10
 131 #define Bletch      2
 132 #define Bndry_mask  ((__uint32_t)0x007fffffL)
 133 #define Bndry_mask1 ((__uint32_t)0x007fffffL)
 134 #define LSB 1
 135 #define Sign_bit    ((__uint32_t)0x80000000L)
 136 #define Log2P       1
 137 #define Tiny0       0
 138 #define Tiny1       1
 139 #define Quick_max   5
 140 #define Int_max     6
 141 #define Infinite(x) (word0(x) == ((__uint32_t)0x7f800000L))
 142 #undef word0
 143 #undef word1
 144 #undef dword0
 145 #undef dword1
 146
 147 #define word0(x) (x.i[0])
 148 #define word1(x) 0
 149 #define dword0(x) word0(x)
 150 #define dword1(x) 0
 151 #else
 152
 153 #define Exp_shift  20
 154 #define Exp_shift1 20
 155 #define Exp_msk1    ((__uint32_t)0x100000L)
 156 #define Exp_msk11   ((__uint32_t)0x100000L)
 157 #define Exp_mask  ((__uint32_t)0x7ff00000L)
 158 #define P 53
 159 #define Bias 1023
 160 #define IEEE_Arith
 161 #define Emin (-1022)
 162 #define Exp_1  ((__uint32_t)0x3ff00000L)
 163 #define Exp_11 ((__uint32_t)0x3ff00000L)
 164 #define Ebits 11
 165 #define Frac_mask  ((__uint32_t)0xfffffL)
 166 #define Frac_mask1 ((__uint32_t)0xfffffL)
 167 #define Ten_pmax 22
 168 #define Bletch 0x10
 169 #define Bndry_mask  ((__uint32_t)0xfffffL)
 170 #define Bndry_mask1 ((__uint32_t)0xfffffL)
 171 #define LSB 1
 172 #define Sign_bit ((__uint32_t)0x80000000L)
 173 #define Log2P 1
 174 #define Tiny0 0
 175 #define Tiny1 1
 176 #define Quick_max 14
 177 #define Int_max 14
 178 #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */
 179
 180 #endif /* !_DOUBLE_IS_32BITS */
 181
 182 #ifndef Flt_Rounds
 183 #ifdef FLT_ROUNDS
 184 #define Flt_Rounds FLT_ROUNDS
 185 #else
 186 #define Flt_Rounds 1
 187 #endif
 188 #endif /*Flt_Rounds*/
 189
 190 #else /* !IEEE_8087 && !IEEE_MC68k */
 191 #undef  Sudden_Underflow
 192 #define Sudden_Underflow
 193 #ifdef IBM
 194 #define Flt_Rounds 0
 195 #define Exp_shift  24
 196 #define Exp_shift1 24
 197 #define Exp_msk1   ((__uint32_t)0x1000000L)
 198 #define Exp_msk11  ((__uint32_t)0x1000000L)
 199 #define Exp_mask  ((__uint32_t)0x7f000000L)
 200 #define P 14
 201 #define Bias 65
 202 #define Exp_1  ((__uint32_t)0x41000000L)
 203 #define Exp_11 ((__uint32_t)0x41000000L)
 204 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
 205 #define Frac_mask  ((__uint32_t)0xffffffL)
 206 #define Frac_mask1 ((__uint32_t)0xffffffL)
 207 #define Bletch 4
 208 #define Ten_pmax 22
 209 #define Bndry_mask  ((__uint32_t)0xefffffL)
 210 #define Bndry_mask1 ((__uint32_t)0xffffffL)
 211 #define LSB 1
 212 #define Sign_bit ((__uint32_t)0x80000000L)
 213 #define Log2P 4
 214 #define Tiny0 ((__uint32_t)0x100000L)
 215 #define Tiny1 0
 216 #define Quick_max 14
 217 #define Int_max 15
 218 #else /* VAX */
 219 #define Flt_Rounds 1
 220 #define Exp_shift  23
 221 #define Exp_shift1 7
 222 #define Exp_msk1    0x80
 223 #define Exp_msk11   ((__uint32_t)0x800000L)
 224 #define Exp_mask  ((__uint32_t)0x7f80L)
 225 #define P 56
 226 #define Bias 129
 227 #define Exp_1  ((__uint32_t)0x40800000L)
 228 #define Exp_11 ((__uint32_t)0x4080L)
 229 #define Ebits 8
 230 #define Frac_mask  ((__uint32_t)0x7fffffL)
 231 #define Frac_mask1 ((__uint32_t)0xffff007fL)
 232 #define Ten_pmax 24
 233 #define Bletch 2
 234 #define Bndry_mask  ((__uint32_t)0xffff007fL)
 235 #define Bndry_mask1 ((__uint32_t)0xffff007fL)
 236 #define LSB ((__uint32_t)0x10000L)
 237 #define Sign_bit ((__uint32_t)0x8000L)
 238 #define Log2P 1
 239 #define Tiny0 0x80
 240 #define Tiny1 0
 241 #define Quick_max 15
 242 #define Int_max 15
 243 #endif
 244 #endif
 245
 246 #ifndef IEEE_Arith
 247 #define ROUND_BIASED
 248 #else
 249 #define Scale_Bit 0x10
 250 #if defined(_DOUBLE_IS_32BITS) && defined(__v800)
 251 #define n_bigtens 2
 252 #else
 253 #define n_bigtens 5
 254 #endif
 255 #endif
 256
 257 #ifdef IBM
 258 #define n_bigtens 3
 259 #endif
 260
 261 #ifdef VAX
 262 #define n_bigtens 2
 263 #endif
 264
 265 #ifndef __NO_INFNAN_CHECK
 266 #define INFNAN_CHECK
 267 #endif
 268
 269 #ifdef RND_PRODQUOT
 270 #define rounded_product(a,b) a = rnd_prod(a, b)
 271 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 272 #ifdef KR_headers
 273 extern double rnd_prod(), rnd_quot();
 274 #else
 275 extern double rnd_prod(double, double), rnd_quot(double, double);
 276 #endif
 277 #else
 278 #define rounded_product(a,b) a *= b
 279 #define rounded_quotient(a,b) a /= b
 280 #endif
 281
 282 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 283 #define Big1 ((__uint32_t)0xffffffffL)
 284
 285 #ifndef Just_16
 286 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
 287  * This makes some inner loops simpler and sometimes saves work
 288  * during multiplications, but it often seems to make things slightly
 289  * slower.  Hence the default is now to store 32 bits per long.
 290  */
 291
 292  #ifndef Pack_32
 293   #define Pack_32
 294  #endif
 295 #else  /* Just_16 */
 296  #ifndef Pack_16
 297   #define Pack_16
 298  #endif
 299 #endif /* Just_16 */
 300
 301 #ifdef Pack_32
 302 #define ULbits 32
 303 #define kshift 5
 304 #define kmask 31
 305 #define ALL_ON 0xffffffff
 306 #else
 307 #define ULbits 16
 308 #define kshift 4
 309 #define kmask 15
 310 #define ALL_ON 0xffff
 311 #endif
 312
 313 #ifdef __cplusplus
 314 extern "C" double strtod(const char *s00, char **se);
 315 extern "C" char *dtoa(double d, int mode, int ndigits,
 316                         int *decpt, int *sign, char **rve);
 317 #endif
 318
 319
 320 typedef struct _Bigint _Bigint;
 321
 322 #define Balloc  _Balloc
 323 #define Bfree   _Bfree
 324 #define multadd __multadd
 325 #define s2b     __s2b
 326 #define lo0bits __lo0bits
 327 #define hi0bits __hi0bits
 328 #define i2b     __i2b
 329 #define mult    __multiply
 330 #define pow5mult        __pow5mult
 331 #define lshift  __lshift
 332 #define match   __match
 333 #define cmp     __mcmp
 334 #define diff    __mdiff
 335 #define ulp     __ulp
 336 #define b2d     __b2d
 337 #define d2b     __d2b
 338 #define ratio   __ratio
 339 #define any_on  __any_on
 340 #define gethex  __gethex
 341 #define copybits        __copybits
 342 #define hexnan  __hexnan
 343
 344 #define eBalloc(__reent_ptr, __len) ({ \
 345    void *__ptr = Balloc(__reent_ptr, __len); \
 346    if (__ptr == NULL) \
 347      __assert_func(__FILE__, __LINE__, (char *)0, "Balloc succeeded"); \
 348    __ptr; \
 349    })
 350
 351 #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG)
 352 #define __get_hexdig(x) __hexdig[x] /* NOTE: must evaluate arg only once */
 353 #else /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */
 354 #define __get_hexdig(x) __hexdig_fun(x)
 355 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */
 356
 357 #define tens __mprec_tens
 358 #define bigtens __mprec_bigtens
 359 #define tinytens __mprec_tinytens
 360
 361 struct _reent ;
 362 struct FPI;
 363 double          ulp (double x);
 364 double          b2d (_Bigint *a , int *e);
 365 _Bigint *       Balloc (struct _reent *p, int k);
 366 void            Bfree (struct _reent *p, _Bigint *v);
 367 _Bigint *       multadd (struct _reent *p, _Bigint *, int, int);
 368 _Bigint *       s2b (struct _reent *, const char*, int, int, __ULong);
 369 _Bigint *       i2b (struct _reent *,int);
 370 _Bigint *       mult (struct _reent *, _Bigint *, _Bigint *);
 371 _Bigint *       pow5mult (struct _reent *, _Bigint *, int k);
 372 int             hi0bits (__ULong);
 373 int             lo0bits (__ULong *);
 374 _Bigint *       d2b (struct _reent *p, double d, int *e, int *bits);
 375 _Bigint *       lshift (struct _reent *p, _Bigint *b, int k);
 376 int             match (const char**, char*);
 377 _Bigint *       diff (struct _reent *p, _Bigint *a, _Bigint *b);
 378 int             cmp (_Bigint *a, _Bigint *b);
 379 int             gethex (struct _reent *p, const char **sp, const struct FPI *fpi, Long *exp, _Bigint **bp, int sign, locale_t loc);
 380 double          ratio (_Bigint *a, _Bigint *b);
 381 __ULong         any_on (_Bigint *b, int k);
 382 void            copybits (__ULong *c, int n, _Bigint *b);
 383 double          _strtod_l (struct _reent *ptr, const char *__restrict s00,
 384                            char **__restrict se, locale_t loc);
 385 #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL)
 386 int             _strtorx_l (struct _reent *, const char *, char **, int,
 387                             void *, locale_t);
 388 int             _strtodg_l (struct _reent *p, const char *s00, char **se,
 389                             struct FPI *fpi, Long *exp, __ULong *bits,
 390                             locale_t);
 391 #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */
 392
 393 #if defined(PREFER_SIZE_OVER_SPEED) || defined(__OPTIMIZE_SIZE__) || defined(_SMALL_HEXDIG)
 394 unsigned char __hexdig_fun (unsigned char);
 395 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */
 396 #ifdef INFNAN_CHECK
 397 int             hexnan (const char **sp, const struct FPI *fpi, __ULong *x0);
 398 #endif
 399
 400 #define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int))
 401
 402 extern const double tinytens[];
 403 extern const double bigtens[];
 404 extern const double tens[];
 405 #if !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG)
 406 extern const unsigned char __hexdig[];
 407 #endif /* !defined(PREFER_SIZE_OVER_SPEED) && !defined(__OPTIMIZE_SIZE__) && !defined(_SMALL_HEXDIG) */
 408
 409
 410 double _mprec_log10 (int);