1 /*********************************************************************** 2 * * 3 * This software is part of the ast package * 4 * Copyright (c) 1985-2009 AT&T Intellectual Property * 5 * and is licensed under the * 6 * Common Public License, Version 1.0 * 7 * by AT&T Intellectual Property * 8 * * 9 * A copy of the License is available at * 10 * http://www.opensource.org/licenses/cpl1.0.txt * 11 * (with md5 checksum 059e8cd6165cb4c31e351f2b69388fd9) * 12 * * 13 * Information and Software Systems Research * 14 * AT&T Research * 15 * Florham Park NJ * 16 * * 17 * Glenn Fowler <gsf@research.att.com> * 18 * David Korn <dgk@research.att.com> * 19 * Phong Vo <kpv@research.att.com> * 20 * * 21 ***********************************************************************/ 22 #pragma prototyped 23 24 /* 25 * frexp/ldexp implementation 26 */ 27 28 #include <ast.h> 29 30 #include "FEATURE/float" 31 32 #if _lib_frexp && _lib_ldexp 33 34 NoN(frexp) 35 36 #else 37 38 #if defined(_ast_dbl_exp_index) && defined(_ast_dbl_exp_shift) 39 40 #define INIT() _ast_dbl_exp_t _pow_ 41 #define pow2(i) (_pow_.f=1,_pow_.e[_ast_dbl_exp_index]+=((i)<<_ast_dbl_exp_shift),_pow_.f) 42 43 #else 44 45 static double pow2tab[DBL_MAX_EXP + 1]; 46 47 static int 48 init(void) 49 { 50 register int x; 51 double g; 52 53 g = 1; 54 for (x = 0; x < elementsof(pow2tab); x++) 55 { 56 pow2tab[x] = g; 57 g *= 2; 58 } 59 return 0; 60 } 61 62 #define INIT() (pow2tab[0]?0:init()) 63 64 #define pow2(i) pow2tab[i] 65 66 #endif 67 68 #if !_lib_frexp 69 70 extern double 71 frexp(double f, int* p) 72 { 73 register int k; 74 register int x; 75 double g; 76 77 INIT(); 78 79 /* 80 * normalize 81 */ 82 83 x = k = DBL_MAX_EXP / 2; 84 if (f < 1) 85 { 86 g = 1.0L / f; 87 for (;;) 88 { 89 k = (k + 1) / 2; 90 if (g < pow2(x)) 91 x -= k; 92 else if (k == 1 && g < pow2(x+1)) 93 break; 94 else 95 x += k; 96 } 97 if (g == pow2(x)) 98 x--; 99 x = -x; 100 } 101 else if (f > 1) 102 { 103 for (;;) 104 { 105 k = (k + 1) / 2; 106 if (f > pow2(x)) 107 x += k; 108 else if (k == 1 && f > pow2(x-1)) 109 break; 110 else 111 x -= k; 112 } 113 if (f == pow2(x)) 114 x++; 115 } 116 else 117 x = 1; 118 *p = x; 119 120 /* 121 * shift 122 */ 123 124 x = -x; 125 if (x < 0) 126 f /= pow2(-x); 127 else if (x < DBL_MAX_EXP) 128 f *= pow2(x); 129 else 130 f = (f * pow2(DBL_MAX_EXP - 1)) * pow2(x - (DBL_MAX_EXP - 1)); 131 return f; 132 } 133 134 #endif 135 136 #if !_lib_ldexp 137 138 extern double 139 ldexp(double f, register int x) 140 { 141 INIT(); 142 if (x < 0) 143 f /= pow2(-x); 144 else if (x < DBL_MAX_EXP) 145 f *= pow2(x); 146 else 147 f = (f * pow2(DBL_MAX_EXP - 1)) * pow2(x - (DBL_MAX_EXP - 1)); 148 return f; 149 } 150 151 #endif 152 153 #endif 154