1 /**************************************************************** 2 * 3 * The author of this software is David M. Gay. 4 * 5 * Copyright (c) 1991, 2000 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 <config.h> 30 #include "ieeefp.h" 31 32 // #include <math.h> 33 // #include <float.h> 34 // #include <errno.h> 35 36 #if defined HAVE_STDINT_H 37 #include <stdint.h> 38 #elif defined HAVE_INTTYPES_H 39 #include <inttypes.h> 40 #endif 41 42 #if defined HAVE_SYS_TYPES_H 43 #include <sys/types.h> 44 #endif 45 46 #if defined HAVE_SYS_CONFIG_H 47 #include <sys/config.h> 48 #endif 49 50 #ifdef __cplusplus 51 extern "C" { 52 #endif 53 54 /* ISO C99 int type declarations */ 55 56 #if !defined HAVE_INT32_DEFINED && defined HAVE_BSD_INT32_DEFINED 57 typedef u_int32_t uint32_t; 58 #endif 59 60 #if !defined HAVE_BSD_INT32_DEFINED && !defined HAVE_INT32_DEFINED 61 // FIXME -- this could have problems with systems that don't define SI to be 4 62 typedef int int32_t __attribute__((mode(SI))); 63 64 /* This is a blatant hack: on Solaris 2.5, pthread.h defines uint32_t 65 in pthread.h, which we sometimes include. We protect our 66 definition the same way Solaris 2.5 does, to avoid redefining it. */ 67 # ifndef _UINT32_T 68 typedef unsigned int uint32_t __attribute__((mode(SI))); 69 # endif 70 #endif 71 72 /* These typedefs are true for the targets running Java. */ 73 74 #ifdef __IEEE_LITTLE_ENDIAN 75 #define IEEE_8087 76 #endif 77 78 #ifdef __IEEE_BIG_ENDIAN 79 #define IEEE_MC68k 80 #endif 81 82 #ifdef __Z8000__ 83 #define Just_16 84 #endif 85 86 #ifdef DEBUG 87 #include "stdio.h" 88 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 89 #endif 90 91 92 #ifdef Unsigned_Shifts 93 #define Sign_Extend(a,b) if (b < 0) a |= (uint32_t)0xffff0000; 94 #else 95 #define Sign_Extend(a,b) /*no-op*/ 96 #endif 97 98 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 99 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. 100 #endif 101 102 /* If we are going to examine or modify specific bits in a double using 103 the word0 and/or word1 macros, then we must wrap the double inside 104 a union. This is necessary to avoid undefined behavior according to 105 the ANSI C spec. */ 106 union double_union 107 { 108 double d; 109 uint32_t i[2]; 110 }; 111 112 #ifdef IEEE_8087 113 #define word0(x) (x.i[1]) 114 #define word1(x) (x.i[0]) 115 #else 116 #define word0(x) (x.i[0]) 117 #define word1(x) (x.i[1]) 118 #endif 119 120 /* The following definition of Storeinc is appropriate for MIPS processors. 121 * An alternative that might be better on some machines is 122 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 123 */ 124 #if defined(IEEE_8087) + defined(VAX) 125 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 126 ((unsigned short *)a)[0] = (unsigned short)c, a++) 127 #else 128 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 129 ((unsigned short *)a)[1] = (unsigned short)c, a++) 130 #endif 131 132 /* #define P DBL_MANT_DIG */ 133 /* Ten_pmax = floor(P*log(2)/log(5)) */ 134 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 135 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 136 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 137 138 #if defined(IEEE_8087) + defined(IEEE_MC68k) 139 #if defined (_DOUBLE_IS_32BITS) 140 #define Exp_shift 23 141 #define Exp_shift1 23 142 #define Exp_msk1 ((uint32_t)0x00800000L) 143 #define Exp_msk11 ((uint32_t)0x00800000L) 144 #define Exp_mask ((uint32_t)0x7f800000L) 145 #define P 24 146 #define Bias 127 147 #if 0 148 #define IEEE_Arith /* it is, but the code doesn't handle IEEE singles yet */ 149 #endif 150 #define Emin (-126) 151 #define Exp_1 ((uint32_t)0x3f800000L) 152 #define Exp_11 ((uint32_t)0x3f800000L) 153 #define Ebits 8 154 #define Frac_mask ((uint32_t)0x007fffffL) 155 #define Frac_mask1 ((uint32_t)0x007fffffL) 156 #define Ten_pmax 10 157 #define Sign_bit ((uint32_t)0x80000000L) 158 #define Ten_pmax 10 159 #define Bletch 2 160 #define Bndry_mask ((uint32_t)0x007fffffL) 161 #define Bndry_mask1 ((uint32_t)0x007fffffL) 162 #define LSB 1 163 #define Sign_bit ((uint32_t)0x80000000L) 164 #define Log2P 1 165 #define Tiny0 0 166 #define Tiny1 1 167 #define Quick_max 5 168 #define Int_max 6 169 #define Infinite(x) (word0(x) == ((uint32_t)0x7f800000L)) 170 #undef word0 171 #undef word1 172 173 #define word0(x) (x.i[0]) 174 #define word1(x) 0 175 #else 176 177 #define Exp_shift 20 178 #define Exp_shift1 20 179 #define Exp_msk1 ((uint32_t)0x100000L) 180 #define Exp_msk11 ((uint32_t)0x100000L) 181 #define Exp_mask ((uint32_t)0x7ff00000L) 182 #define P 53 183 #define Bias 1023 184 #define IEEE_Arith 185 #define Emin (-1022) 186 #define Exp_1 ((uint32_t)0x3ff00000L) 187 #define Exp_11 ((uint32_t)0x3ff00000L) 188 #define Ebits 11 189 #define Frac_mask ((uint32_t)0xfffffL) 190 #define Frac_mask1 ((uint32_t)0xfffffL) 191 #define Ten_pmax 22 192 #define Bletch 0x10 193 #define Bndry_mask ((uint32_t)0xfffffL) 194 #define Bndry_mask1 ((uint32_t)0xfffffL) 195 #define LSB 1 196 #define Sign_bit ((uint32_t)0x80000000L) 197 #define Log2P 1 198 #define Tiny0 0 199 #define Tiny1 1 200 #define Quick_max 14 201 #define Int_max 14 202 #define Infinite(x) (word0(x) == ((uint32_t)0x7ff00000L)) /* sufficient test for here */ 203 #endif 204 205 #else 206 #undef Sudden_Underflow 207 #define Sudden_Underflow 208 #ifdef IBM 209 #define Exp_shift 24 210 #define Exp_shift1 24 211 #define Exp_msk1 ((uint32_t)0x1000000L) 212 #define Exp_msk11 ((uint32_t)0x1000000L) 213 #define Exp_mask ((uint32_t)0x7f000000L) 214 #define P 14 215 #define Bias 65 216 #define Exp_1 ((uint32_t)0x41000000L) 217 #define Exp_11 ((uint32_t)0x41000000L) 218 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ 219 #define Frac_mask ((uint32_t)0xffffffL) 220 #define Frac_mask1 ((uint32_t)0xffffffL) 221 #define Bletch 4 222 #define Ten_pmax 22 223 #define Bndry_mask ((uint32_t)0xefffffL) 224 #define Bndry_mask1 ((uint32_t)0xffffffL) 225 #define LSB 1 226 #define Sign_bit ((uint32_t)0x80000000L) 227 #define Log2P 4 228 #define Tiny0 ((uint32_t)0x100000L) 229 #define Tiny1 0 230 #define Quick_max 14 231 #define Int_max 15 232 #else /* VAX */ 233 #define Exp_shift 23 234 #define Exp_shift1 7 235 #define Exp_msk1 0x80 236 #define Exp_msk11 ((uint32_t)0x800000L) 237 #define Exp_mask ((uint32_t)0x7f80L) 238 #define P 56 239 #define Bias 129 240 #define Exp_1 ((uint32_t)0x40800000L) 241 #define Exp_11 ((uint32_t)0x4080L) 242 #define Ebits 8 243 #define Frac_mask ((uint32_t)0x7fffffL) 244 #define Frac_mask1 ((uint32_t)0xffff007fL) 245 #define Ten_pmax 24 246 #define Bletch 2 247 #define Bndry_mask ((uint32_t)0xffff007fL) 248 #define Bndry_mask1 ((uint32_t)0xffff007fL) 249 #define LSB ((uint32_t)0x10000L) 250 #define Sign_bit ((uint32_t)0x8000L) 251 #define Log2P 1 252 #define Tiny0 0x80 253 #define Tiny1 0 254 #define Quick_max 15 255 #define Int_max 15 256 #endif 257 #endif 258 259 #ifndef IEEE_Arith 260 #define ROUND_BIASED 261 #endif 262 263 #ifdef RND_PRODQUOT 264 #define rounded_product(a,b) a = rnd_prod(a, b) 265 #define rounded_quotient(a,b) a = rnd_quot(a, b) 266 #ifdef KR_headers 267 extern double rnd_prod(), rnd_quot(); 268 #else 269 extern double rnd_prod(double, double), rnd_quot(double, double); 270 #endif 271 #else 272 #define rounded_product(a,b) a *= b 273 #define rounded_quotient(a,b) a /= b 274 #endif 275 276 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 277 #define Big1 ((uint32_t)0xffffffffL) 278 279 #ifndef Just_16 280 /* When Pack_32 is not defined, we store 16 bits per 32-bit long. 281 * This makes some inner loops simpler and sometimes saves work 282 * during multiplications, but it often seems to make things slightly 283 * slower. Hence the default is now to store 32 bits per long. 284 */ 285 286 #ifndef Pack_32 287 #if SIZEOF_VOID_P != 8 288 #define Pack_32 289 #endif 290 #endif 291 #endif 292 293 294 #define MAX_BIGNUMS 16 295 #define MAX_BIGNUM_WDS 32 296 297 struct _Jv_Bigint 298 { 299 struct _Jv_Bigint *_next; 300 int _k, _maxwds, _sign, _wds; 301 unsigned long _x[MAX_BIGNUM_WDS]; 302 }; 303 304 305 #define _PTR void * 306 #define _AND , 307 #define _NOARGS void 308 #define _CONST const 309 #define _VOLATILE volatile 310 #define _SIGNED signed 311 #define _DOTS , ... 312 #define _VOID void 313 #define _EXFUN(name, proto) name proto 314 #define _DEFUN(name, arglist, args) name(args) 315 #define _DEFUN_VOID(name) name(_NOARGS) 316 #define _CAST_VOID (void) 317 318 319 struct _Jv_reent 320 { 321 /* local copy of errno */ 322 int _errno; 323 324 /* used by mprec routines */ 325 struct _Jv_Bigint *_result; 326 int _result_k; 327 struct _Jv_Bigint *_p5s; 328 329 struct _Jv_Bigint _freelist[MAX_BIGNUMS]; 330 int _allocation_map; 331 332 int num; 333 }; 334 335 336 typedef struct _Jv_Bigint _Jv_Bigint; 337 338 #define Balloc _Jv_Balloc 339 #define Bfree _Jv_Bfree 340 #define multadd _Jv_multadd 341 #define s2b _Jv_s2b 342 #define lo0bits _Jv_lo0bits 343 #define hi0bits _Jv_hi0bits 344 #define i2b _Jv_i2b 345 #define mult _Jv_mult 346 #define pow5mult _Jv_pow5mult 347 #define lshift _Jv_lshift 348 #define cmp _Jv__mcmp 349 #define diff _Jv__mdiff 350 #define ulp _Jv_ulp 351 #define b2d _Jv_b2d 352 #define d2b _Jv_d2b 353 #define ratio _Jv_ratio 354 355 #define tens _Jv__mprec_tens 356 #define bigtens _Jv__mprec_bigtens 357 #define tinytens _Jv__mprec_tinytens 358 359 #define _dtoa _Jv_dtoa 360 #define _dtoa_r _Jv_dtoa_r 361 #define _strtod_r _Jv_strtod_r 362 363 extern double _EXFUN(_strtod_r, (struct _Jv_reent *ptr, const char *s00, char **se)); 364 extern char* _EXFUN(_dtoa_r, (struct _Jv_reent *ptr, double d, 365 int mode, int ndigits, int *decpt, int *sign, 366 char **rve, int float_type)); 367 void _EXFUN(_dtoa, (double d, int mode, int ndigits, int *decpt, int *sign, 368 char **rve, char *buf, int float_type)); 369 370 double _EXFUN(ulp,(double x)); 371 double _EXFUN(b2d,(_Jv_Bigint *a , int *e)); 372 _Jv_Bigint * _EXFUN(Balloc,(struct _Jv_reent *p, int k)); 373 void _EXFUN(Bfree,(struct _Jv_reent *p, _Jv_Bigint *v)); 374 _Jv_Bigint * _EXFUN(multadd,(struct _Jv_reent *p, _Jv_Bigint *, int, int)); 375 _Jv_Bigint * _EXFUN(s2b,(struct _Jv_reent *, const char*, int, int, unsigned long)); 376 _Jv_Bigint * _EXFUN(i2b,(struct _Jv_reent *,int)); 377 _Jv_Bigint * _EXFUN(mult, (struct _Jv_reent *, _Jv_Bigint *, _Jv_Bigint *)); 378 _Jv_Bigint * _EXFUN(pow5mult, (struct _Jv_reent *, _Jv_Bigint *, int k)); 379 int _EXFUN(hi0bits,(unsigned long)); 380 int _EXFUN(lo0bits,(unsigned long *)); 381 _Jv_Bigint * _EXFUN(d2b,(struct _Jv_reent *p, double d, int *e, int *bits)); 382 _Jv_Bigint * _EXFUN(lshift,(struct _Jv_reent *p, _Jv_Bigint *b, int k)); 383 _Jv_Bigint * _EXFUN(diff,(struct _Jv_reent *p, _Jv_Bigint *a, _Jv_Bigint *b)); 384 int _EXFUN(cmp,(_Jv_Bigint *a, _Jv_Bigint *b)); 385 386 double _EXFUN(ratio,(_Jv_Bigint *a, _Jv_Bigint *b)); 387 #define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(long) + 2*sizeof(int)) 388 389 #if defined(_DOUBLE_IS_32BITS) && defined(__v800) 390 #define n_bigtens 2 391 #else 392 #define n_bigtens 5 393 #endif 394 395 extern _CONST double tinytens[]; 396 extern _CONST double bigtens[]; 397 extern _CONST double tens[]; 398 399 #ifdef __cplusplus 400 } 401 #endif 402