1/*********************************************************************** 2* * 3* This software is part of the ast package * 4* Copyright (c) 1985-2011 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 34NoN(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 45static double pow2tab[DBL_MAX_EXP + 1]; 46 47static int 48init(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 70extern double 71frexp(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 138extern double 139ldexp(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