sysdeps/m68k/fpu/e_pow.c

   1 /* Copyright (C) 1997 Free Software Foundation, Inc.
   2    This file is part of the GNU C Library.
   3
   4    The GNU C Library is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU Library General Public License as
   6    published by the Free Software Foundation; either version 2 of the
   7    License, or (at your option) any later version.
   8
   9    The GNU C Library is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    Library General Public License for more details.
  13
  14    You should have received a copy of the GNU Library General Public
  15    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  16    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17    Boston, MA 02111-1307, USA.  */
  18
  19 #include <math.h>
  20 #include "math_private.h"
  21
  22 #ifndef SUFF
  23 #define SUFF
  24 #endif
  25 #ifndef float_type
  26 #define float_type double
  27 #endif
  28
  29 #define CONCATX(a,b) __CONCAT(a,b)
  30 #define s(name) CONCATX(name,SUFF)
  31 #define m81(func) __m81_u(s(func))
  32
  33 float_type
  34 s(__ieee754_pow) (float_type x, float_type y)
  35 {
  36   float_type z;
  37   float_type ax;
  38   unsigned long x_cond, y_cond;
  39
  40   y_cond = __m81_test (y);
  41   if (y_cond & __M81_COND_ZERO)
  42     return 1.0;
  43
  44   x_cond = __m81_test (x);
  45   if ((x_cond | y_cond) & __M81_COND_NAN)
  46     return x + y;
  47
  48   if (y_cond & __M81_COND_INF)
  49     {
  50       ax = s(fabs) (x);
  51       if (ax == 1)
  52         return y - y;
  53       if (ax > 1)
  54         return y_cond & __M81_COND_NEG ? 0 : y;
  55       else
  56         return y_cond & __M81_COND_NEG ? -y : 0;
  57     }
  58
  59   if (s(fabs) (y) == 1)
  60     return y_cond & __M81_COND_NEG ? 1 / x : x;
  61
  62   if (y == 2)
  63     return x * x;
  64   if (y == 0.5 && !(x_cond & __M81_COND_NEG))
  65     return m81(__ieee754_sqrt) (x);
  66
  67   if (x == 10.0)
  68     {
  69       __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y));
  70       return z;
  71     }
  72   if (x == 2.0)
  73     {
  74       __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y));
  75       return z;
  76     }
  77
  78   ax = s(fabs) (x);
  79   if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1)
  80     {
  81       z = ax;
  82       if (y_cond & __M81_COND_NEG)
  83         z = 1 / z;
  84       if (x_cond & __M81_COND_NEG)
  85         {
  86           if (y != m81(__rint) (y))
  87             {
  88               if (x == -1)
  89                 z = (z - z) / (z - z);
  90             }
  91           else
  92             goto maybe_negate;
  93         }
  94       return z;
  95     }
  96
  97   if (x_cond & __M81_COND_NEG)
  98     {
  99       if (y == m81(__rint) (y))
 100         {
 101           z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x));
 102         maybe_negate:
 103           /* We always use the long double format, since y is already in
 104              this format and rounding won't change the result.  */
 105           {
 106             int32_t exponent;
 107             u_int32_t i0, i1;
 108             GET_LDOUBLE_WORDS (exponent, i0, i1, y);
 109             exponent = (exponent & 0x7fff) - 0x3fff;
 110             if (exponent <= 31
 111                 ? i0 & (1 << (31 - exponent))
 112                 : (exponent <= 63
 113                    && i1 & (1 << (63 - exponent))))
 114               z = -z;
 115           }
 116         }
 117       else
 118         z = (y - y) / (y - y);
 119     }
 120   else
 121     z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x));
 122   return z;
 123 }