lib/frexpl.c

   1 /* Emulation for frexpl.
   2    Contributed by Paolo Bonzini
   3
   4    Copyright 2002, 2003 Free Software Foundation, Inc.
   5
   6    This file is part of gnulib.
   7
   8    This program is free software; you can redistribute it and/or modify
   9    it under the terms of the GNU General Public License as published by
  10    the Free Software Foundation; either version 2, or (at your option)
  11    any later version.
  12
  13    This program is distributed in the hope that it will be useful,
  14    but WITHOUT ANY WARRANTY; without even the implied warranty of
  15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16    GNU General Public License for more details.
  17
  18    You should have received a copy of the GNU General Public License along
  19    with this program; if not, write to the Free Software Foundation,
  20    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
  21
  22 #include <float.h>
  23 #include <math.h>
  24
  25 #include "mathl.h"
  26
  27 /* Binary search.  Quite inefficient but portable. */
  28 long double
  29 frexpl(long double x, int *exp)
  30 {
  31   long double exponents[20], *next;
  32   int exponent, bit;
  33
  34   /* Check for zero, nan and infinity. */
  35   if (x != x || x + x == x )
  36     {
  37       *exp = 0;
  38       return x;
  39     }
  40
  41   if (x < 0)
  42     return -frexpl(-x, exp);
  43
  44   exponent = 0;
  45   if (x > 1.0)
  46     {
  47       for (next = exponents, exponents[0] = 2.0L, bit = 1;
  48            *next <= x + x;
  49            bit <<= 1, next[1] = next[0] * next[0], next++);
  50
  51       for (; next >= exponents; bit >>= 1, next--)
  52         if (x + x >= *next)
  53           {
  54             x /= *next;
  55             exponent |= bit;
  56           }
  57
  58     }
  59
  60   else if (x < 0.5)
  61     {
  62       for (next = exponents, exponents[0] = 0.5L, bit = 1;
  63            *next > x;
  64            bit <<= 1, next[1] = next[0] * next[0], next++);
  65
  66       for (; next >= exponents; bit >>= 1, next--)
  67         if (x < *next)
  68           {
  69             x /= *next;
  70             exponent |= bit;
  71           }
  72
  73       exponent = -exponent;
  74     }
  75
  76   *exp = exponent;
  77   return x;
  78 }
  79
  80 #if 0
  81 int
  82 main()
  83 {
  84   long double x;
  85   int y;
  86   x = frexpl(0.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
  87   x = frexpl(1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
  88   x = frexpl(-1.0L / 0.0L, &y); printf ("%.6Lg %d\n", x, y);
  89   x = frexpl(0.5L, &y); printf ("%.6Lg %d\n", x, y);
  90   x = frexpl(0.75L, &y); printf ("%.6Lg %d\n", x, y);
  91   x = frexpl(3.6L, &y); printf ("%.6Lg %d\n", x, y);
  92   x = frexpl(17.8L, &y); printf ("%.6Lg %d\n", x, y);
  93   x = frexpl(8.0L, &y); printf ("%.6Lg %d\n", x, y);
  94   x = frexpl(0.3L, &y); printf ("%.6Lg %d\n", x, y);
  95   x = frexpl(0.49L, &y); printf ("%.6Lg %d\n", x, y);
  96   x = frexpl(0.049L, &y); printf ("%.6Lg %d\n", x, y);
  97   x = frexpl(0.0245L, &y); printf ("%.6Lg %d\n", x, y);
  98   x = frexpl(0.0625L, &y); printf ("%.6Lg %d\n", x, y);
  99 }
 100 #endif
 101