contrib/openbsd_libm/src/s_logbl.c

   1 /*      $OpenBSD: s_logbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $    */
   2 /*
   3  * From: @(#)s_ilogb.c 5.1 93/09/24
   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 #include <sys/types.h>
  15 #include <machine/ieee.h>
  16 #include <float.h>
  17 #include <limits.h>
  18 #include <math.h>
  19
  20 long double
  21 logbl(long double x)
  22 {
  23         union {
  24                 long double e;
  25                 struct ieee_ext bits;
  26         } u;
  27         unsigned long m;
  28         int b;
  29
  30         u.e = x;
  31         if (u.bits.ext_exp == 0) {
  32                 if ((u.bits.ext_fracl
  33 #ifdef EXT_FRACLMBITS
  34                         | u.bits.ext_fraclm
  35 #endif /* EXT_FRACLMBITS */
  36 #ifdef EXT_FRACHMBITS
  37                         | u.bits.ext_frachm
  38 #endif /* EXT_FRACHMBITS */
  39                         | u.bits.ext_frach) == 0) {     /* x == 0 */
  40                         u.bits.ext_sign = 1;
  41                         return (1.0L / u.e);
  42                 }
  43                 /* denormalized */
  44                 if (u.bits.ext_frach == 0
  45 #ifdef EXT_FRACHMBITS
  46                         && u.bits.ext_frachm == 0
  47 #endif
  48                         ) {
  49                         m = 1lu << (EXT_FRACLBITS - 1);
  50                         for (b = EXT_FRACHBITS; !(u.bits.ext_fracl & m); m >>= 1)
  51                                 b++;
  52 #if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS)
  53                         m = 1lu << (EXT_FRACLMBITS - 1);
  54                         for (b += EXT_FRACHMBITS; !(u.bits.ext_fraclm & m);
  55                                 m >>= 1)
  56                                 b++;
  57 #endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */
  58                 } else {
  59                         m = 1lu << (EXT_FRACHBITS - 1);
  60                         for (b = 0; !(u.bits.ext_frach & m); m >>= 1)
  61                                 b++;
  62 #ifdef EXT_FRACHMBITS
  63                         m = 1lu << (EXT_FRACHMBITS - 1);
  64                         for (; !(u.bits.ext_frachm & m); m >>= 1)
  65                                 b++;
  66 #endif /* EXT_FRACHMBITS */
  67                 }
  68 #ifdef EXT_IMPLICIT_NBIT
  69                 b++;
  70 #endif
  71                 return ((long double)(LDBL_MIN_EXP - b - 1));
  72         }
  73         if (u.bits.ext_exp < (LDBL_MAX_EXP << 1) - 1)   /* normal */
  74                 return ((long double)(u.bits.ext_exp - LDBL_MAX_EXP + 1));
  75         else                                            /* +/- inf or nan */
  76                 return (x * x);
  77 }