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 * frexpl/ldexpl implementation 26 */ 27 28#include <ast.h> 29 30#include "FEATURE/float" 31 32#if _lib_frexpl && _lib_ldexpl 33 34NoN(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 49static _ast_fltmax_t pow2tab[LDBL_MAX_EXP + 1]; 50 51static int 52init(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 76extern _ast_fltmax_t 77frexpl(_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 146extern _ast_fltmax_t 147ldexpl(_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