usr/src/lib/libast/common/comp/frexp.c

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