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 * frexpl/ldexpl implementation 26 */ 27 28 #include <ast.h> 29 30 #include "FEATURE/float" 31 32 #if _lib_frexpl && _lib_ldexpl 33 34 NoN(frexpl) 35 36 #else 37 38 #ifndef LDBL_MAX_EXP 39 #define LDBL_MAX_EXP DBL_MAX_EXP 40 #endif 41 42 #if defined(_ast_fltmax_exp_index) && defined(_ast_fltmax_exp_shift) 43 44 #define INIT() _ast_fltmax_exp_t _pow_ 45 #define pow2(i) (_pow_.f=1,_pow_.e[_ast_fltmax_exp_index]+=((i)<<_ast_fltmax_exp_shift),_pow_.f) 46 47 #else 48 49 static _ast_fltmax_t pow2tab[LDBL_MAX_EXP + 1]; 50 51 static int 52 init(void) 53 { 54 register int x; 55 _ast_fltmax_t g; 56 57 g = 1; 58 for (x = 0; x < elementsof(pow2tab); x++) 59 { 60 pow2tab[x] = g; 61 g *= 2; 62 } 63 return 0; 64 } 65 66 #define INIT() (pow2tab[0]?0:init()) 67 68 #define pow2(i) (pow2tab[i]) 69 70 #endif 71 72 #if !_lib_frexpl 73 74 #undef frexpl 75 76 extern _ast_fltmax_t 77 frexpl(_ast_fltmax_t f, int* p) 78 { 79 register int k; 80 register int x; 81 _ast_fltmax_t g; 82 83 INIT(); 84 85 /* 86 * normalize 87 */ 88 89 x = k = LDBL_MAX_EXP / 2; 90 if (f < 1) 91 { 92 g = 1.0L / f; 93 for (;;) 94 { 95 k = (k + 1) / 2; 96 if (g < pow2(x)) 97 x -= k; 98 else if (k == 1 && g < pow2(x+1)) 99 break; 100 else 101 x += k; 102 } 103 if (g == pow2(x)) 104 x--; 105 x = -x; 106 } 107 else if (f > 1) 108 { 109 for (;;) 110 { 111 k = (k + 1) / 2; 112 if (f > pow2(x)) 113 x += k; 114 else if (k == 1 && f > pow2(x-1)) 115 break; 116 else 117 x -= k; 118 } 119 if (f == pow2(x)) 120 x++; 121 } 122 else 123 x = 1; 124 *p = x; 125 126 /* 127 * shift 128 */ 129 130 x = -x; 131 if (x < 0) 132 f /= pow2(-x); 133 else if (x < LDBL_MAX_EXP) 134 f *= pow2(x); 135 else 136 f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1)); 137 return f; 138 } 139 140 #endif 141 142 #if !_lib_ldexpl 143 144 #undef ldexpl 145 146 extern _ast_fltmax_t 147 ldexpl(_ast_fltmax_t f, register int x) 148 { 149 INIT(); 150 if (x < 0) 151 f /= pow2(-x); 152 else if (x < LDBL_MAX_EXP) 153 f *= pow2(x); 154 else 155 f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1)); 156 return f; 157 } 158 159 #endif 160 161 #endif 162