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