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