1 /* 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 #include "softmath_private.h" 13 #include <inttypes.h> 14 15 static const double 16 bp[] = {1.0, 1.5,}, 17 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ 18 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ 19 Zero[] = {0.0, -0.0,}, 20 zero = 0.0, 21 one = 1.0, 22 two = 2.0, 23 two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ 24 huge = 1.0e300, 25 tiny = 1.0e-300, 26 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ 27 L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ 28 L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ 29 L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ 30 L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ 31 L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ 32 L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ 33 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 34 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 35 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 36 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 37 P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ 38 lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 39 lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ 40 lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ 41 ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ 42 cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ 43 cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ 44 cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ 45 ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ 46 ivln2_h= 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ 47 ivln2_l= 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ 48 49 typedef unsigned int u_int32_t; 50 51 typedef union 52 { 53 double value; 54 struct 55 { 56 u_int32_t lsw; 57 u_int32_t msw; 58 } parts; 59 } ieee_double_shape_type; 60 61 /* Get two 32 bit ints from a double. */ 62 63 #define EXTRACT_WORDS(ix0,ix1,d) \ 64 do { \ 65 ieee_double_shape_type ew_u; \ 66 ew_u.value = (d); \ 67 (ix0) = ew_u.parts.msw; \ 68 (ix1) = ew_u.parts.lsw; \ 69 } while (0) 70 71 /* Get the most significant 32 bit int from a double. */ 72 73 #define GET_HIGH_WORD(i,d) \ 74 do { \ 75 ieee_double_shape_type gh_u; \ 76 gh_u.value = (d); \ 77 (i) = gh_u.parts.msw; \ 78 } while (0) 79 80 /* Get the less significant 32 bit int from a double. */ 81 82 #define GET_LOW_WORD(i,d) \ 83 do { \ 84 ieee_double_shape_type gl_u; \ 85 gl_u.value = (d); \ 86 (i) = gl_u.parts.lsw; \ 87 } while (0) 88 89 /* Set a double from two 32 bit ints. */ 90 91 #define INSERT_WORDS(d,ix0,ix1) \ 92 do { \ 93 ieee_double_shape_type iw_u; \ 94 iw_u.parts.msw = (ix0); \ 95 iw_u.parts.lsw = (ix1); \ 96 (d) = iw_u.value; \ 97 } while (0) 98 99 /* Set the more significant 32 bits of a double from an int. */ 100 101 #define SET_HIGH_WORD(d,v) \ 102 do { \ 103 ieee_double_shape_type sh_u; \ 104 sh_u.value = (d); \ 105 sh_u.parts.msw = (v); \ 106 (d) = sh_u.value; \ 107 } while (0) 108 109 /* Set the less significant 32 bits of a double from an int. */ 110 111 #define SET_LOW_WORD(d,v) \ 112 do { \ 113 ieee_double_shape_type sl_u; \ 114 sl_u.value = (d); \ 115 sl_u.parts.lsw = (v); \ 116 (d) = sl_u.value; \ 117 } while (0) 118