newlib/libm/common/s_logb.c

   1 /* 2009 for Newlib:  Sun's s_ilogb.c converted to be s_logb.c.  */
   2 /* @(#)s_ilogb.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13 /*
  14 FUNCTION
  15        <<logb>>, <<logbf>>---get exponent of floating-point number
  16 INDEX
  17         logb
  18 INDEX
  19         logbf
  20
  21 SYNOPSIS
  22         #include <math.h>
  23         double logb(double <[x]>);
  24         float logbf(float <[x]>);
  25
  26 DESCRIPTION
  27 The <<logb>> functions extract the exponent of <[x]>, as a signed integer value
  28 in floating-point format.  If <[x]> is subnormal it is treated as though it were
  29 normalized; thus, for positive finite <[x]>,
  30 @ifnottex
  31 1 <= (<[x]> * FLT_RADIX to the power (-logb(<[x]>))) < FLT_RADIX.
  32 @end ifnottex
  33 @tex
  34 $1 \leq ( x \cdot FLT\_RADIX ^ {-logb(x)} ) < FLT\_RADIX$.
  35 @end tex
  36 A domain error may occur if the argument is zero.
  37 In this floating-point implementation, FLT_RADIX is 2.  Which also means
  38 that for finite <[x]>, <<logb>>(<[x]>) = <<floor>>(<<log2>>(<<fabs>>(<[x]>))).
  39
  40 All nonzero, normal numbers can be described as
  41 @ifnottex
  42 <[m]> * 2**<[p]>, where 1.0 <= <[m]> < 2.0.
  43 @end ifnottex
  44 @tex
  45 $m \cdot 2^p$, where $1.0 \leq m < 2.0$.
  46 @end tex
  47 The <<logb>> functions examine the argument <[x]>, and return <[p]>.
  48 The <<frexp>> functions are similar to the <<logb>> functions, but
  49 returning <[m]> adjusted to the interval [.5, 1) or 0, and <[p]>+1.
  50
  51 RETURNS
  52 @comment Formatting note:  "$@" forces a new line
  53 When <[x]> is:@*
  54 +inf or -inf, +inf is returned;@*
  55 NaN, NaN is returned;@*
  56 0, -inf is returned, and the divide-by-zero exception is raised;@*
  57 otherwise, the <<logb>> functions return the signed exponent of <[x]>.
  58
  59 PORTABILITY
  60 ANSI C, POSIX
  61
  62 SEEALSO
  63 frexp, ilogb
  64 */
  65
  66 /* double logb(double x)
  67  * return the binary exponent of non-zero x
  68  * logb(0) = -inf, raise divide-by-zero floating point exception
  69  * logb(+inf|-inf) = +inf (no signal is raised)
  70  * logb(NaN) = NaN (no signal is raised)
  71  * Per C99 recommendation, a NaN argument is returned unchanged.
  72  */
  73
  74 #include "fdlibm.h"
  75
  76 #ifndef _DOUBLE_IS_32BITS
  77
  78 double
  79 #ifdef __STDC__
  80 logb(double x)
  81 #else
  82 logb(x)
  83 double x;
  84 #endif
  85 {
  86         __int32_t hx,lx,ix;
  87
  88         EXTRACT_WORDS(hx,lx,x);
  89         hx &= 0x7fffffff;               /* high |x| */
  90         if(hx<0x00100000) {             /* 0 or subnormal */
  91             if((hx|lx)==0)  {
  92                 double  xx;
  93                 /* arg==0:  return -inf and raise divide-by-zero exception */
  94                 INSERT_WORDS(xx,hx,lx); /* +0.0 */
  95                 return -1./xx;  /* logb(0) = -inf */
  96                 }
  97             else                        /* subnormal x */
  98                 if(hx==0) {
  99                     for (ix = -1043; lx>0; lx<<=1) ix -=1;
 100                 } else {
 101                     for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
 102                 }
 103             return (double) ix;
 104         }
 105         else if (hx<0x7ff00000) return (hx>>20)-1023;   /* normal # */
 106         else if (hx>0x7ff00000 || lx)  return x;        /* x==NaN */
 107         else  return HUGE_VAL;  /* x==inf (+ or -) */
 108 }
 109
 110 #endif /* _DOUBLE_IS_32BITS */