14887Schin /*********************************************************************** 24887Schin * * 34887Schin * This software is part of the ast package * 4*12068SRoger.Faulkner@Oracle.COM * Copyright (c) 1985-2010 AT&T Intellectual Property * 54887Schin * and is licensed under the * 64887Schin * Common Public License, Version 1.0 * 78462SApril.Chin@Sun.COM * by AT&T Intellectual Property * 84887Schin * * 94887Schin * A copy of the License is available at * 104887Schin * http://www.opensource.org/licenses/cpl1.0.txt * 114887Schin * (with md5 checksum 059e8cd6165cb4c31e351f2b69388fd9) * 124887Schin * * 134887Schin * Information and Software Systems Research * 144887Schin * AT&T Research * 154887Schin * Florham Park NJ * 164887Schin * * 174887Schin * Glenn Fowler <gsf@research.att.com> * 184887Schin * David Korn <dgk@research.att.com> * 194887Schin * Phong Vo <kpv@research.att.com> * 204887Schin * * 214887Schin ***********************************************************************/ 224887Schin #pragma prototyped 234887Schin 244887Schin /* 254887Schin * frexpl/ldexpl implementation 264887Schin */ 274887Schin 284887Schin #include <ast.h> 294887Schin 304887Schin #include "FEATURE/float" 314887Schin 324887Schin #if _lib_frexpl && _lib_ldexpl 334887Schin 344887Schin NoN(frexpl) 354887Schin 364887Schin #else 374887Schin 384887Schin #ifndef LDBL_MAX_EXP 394887Schin #define LDBL_MAX_EXP DBL_MAX_EXP 404887Schin #endif 414887Schin 424887Schin #if defined(_ast_fltmax_exp_index) && defined(_ast_fltmax_exp_shift) 434887Schin 444887Schin #define INIT() _ast_fltmax_exp_t _pow_ 454887Schin #define pow2(i) (_pow_.f=1,_pow_.e[_ast_fltmax_exp_index]+=((i)<<_ast_fltmax_exp_shift),_pow_.f) 464887Schin 474887Schin #else 484887Schin 494887Schin static _ast_fltmax_t pow2tab[LDBL_MAX_EXP + 1]; 504887Schin 514887Schin static int 524887Schin init(void) 534887Schin { 544887Schin register int x; 554887Schin _ast_fltmax_t g; 564887Schin 574887Schin g = 1; 584887Schin for (x = 0; x < elementsof(pow2tab); x++) 594887Schin { 604887Schin pow2tab[x] = g; 614887Schin g *= 2; 624887Schin } 634887Schin return 0; 644887Schin } 654887Schin 664887Schin #define INIT() (pow2tab[0]?0:init()) 674887Schin 684887Schin #define pow2(i) (pow2tab[i]) 694887Schin 704887Schin #endif 714887Schin 724887Schin #if !_lib_frexpl 734887Schin 744887Schin #undef frexpl 754887Schin 764887Schin extern _ast_fltmax_t 774887Schin frexpl(_ast_fltmax_t f, int* p) 784887Schin { 794887Schin register int k; 804887Schin register int x; 814887Schin _ast_fltmax_t g; 824887Schin 834887Schin INIT(); 844887Schin 854887Schin /* 864887Schin * normalize 874887Schin */ 884887Schin 894887Schin x = k = LDBL_MAX_EXP / 2; 904887Schin if (f < 1) 914887Schin { 924887Schin g = 1.0L / f; 934887Schin for (;;) 944887Schin { 954887Schin k = (k + 1) / 2; 964887Schin if (g < pow2(x)) 974887Schin x -= k; 984887Schin else if (k == 1 && g < pow2(x+1)) 994887Schin break; 1004887Schin else 1014887Schin x += k; 1024887Schin } 1034887Schin if (g == pow2(x)) 1044887Schin x--; 1054887Schin x = -x; 1064887Schin } 1074887Schin else if (f > 1) 1084887Schin { 1094887Schin for (;;) 1104887Schin { 1114887Schin k = (k + 1) / 2; 1124887Schin if (f > pow2(x)) 1134887Schin x += k; 1144887Schin else if (k == 1 && f > pow2(x-1)) 1154887Schin break; 1164887Schin else 1174887Schin x -= k; 1184887Schin } 1194887Schin if (f == pow2(x)) 1204887Schin x++; 1214887Schin } 1224887Schin else 1234887Schin x = 1; 1244887Schin *p = x; 1254887Schin 1264887Schin /* 1274887Schin * shift 1284887Schin */ 1294887Schin 1304887Schin x = -x; 1314887Schin if (x < 0) 1324887Schin f /= pow2(-x); 1334887Schin else if (x < LDBL_MAX_EXP) 1344887Schin f *= pow2(x); 1354887Schin else 1364887Schin f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1)); 1374887Schin return f; 1384887Schin } 1394887Schin 1404887Schin #endif 1414887Schin 1424887Schin #if !_lib_ldexpl 1434887Schin 1444887Schin #undef ldexpl 1454887Schin 1464887Schin extern _ast_fltmax_t 1474887Schin ldexpl(_ast_fltmax_t f, register int x) 1484887Schin { 1494887Schin INIT(); 1504887Schin if (x < 0) 1514887Schin f /= pow2(-x); 1524887Schin else if (x < LDBL_MAX_EXP) 1534887Schin f *= pow2(x); 1544887Schin else 1554887Schin f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1)); 1564887Schin return f; 1574887Schin } 1584887Schin 1594887Schin #endif 1604887Schin 1614887Schin #endif 162