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 #define __LIBC_INTERNAL_MATH_INLINES
  20 #include <math.h>
  21 #include "math_private.h"
  22
  23 #ifndef SUFF
  24 #define SUFF
  25 #endif
  26 #ifndef float_type
  27 #define float_type double
  28 #endif
  29
  30 #define CONCATX(a,b) __CONCAT(a,b)
  31 #define s(name) CONCATX(name,SUFF)
  32 #define m81(func) __m81_u(s(func))
  33
  34 float_type
  35 s(__ieee754_pow) (float_type x, float_type y)
  36 {
  37   float_type z;
  38   float_type ax;
  39   unsigned long x_cond, y_cond;
  40
  41   y_cond = __m81_test (y);
  42   if (y_cond & __M81_COND_ZERO)
  43     return 1.0;
  44
  45   x_cond = __m81_test (x);
  46   if ((x_cond | y_cond) & __M81_COND_NAN)
  47     return x + y;
  48
  49   if (y_cond & __M81_COND_INF)
  50     {
  51       ax = s(fabs) (x);
  52       if (ax == 1)
  53         return y - y;
  54       if (ax > 1)
  55         return y_cond & __M81_COND_NEG ? 0 : y;
  56       else
  57         return y_cond & __M81_COND_NEG ? -y : 0;
  58     }
  59
  60   if (s(fabs) (y) == 1)
  61     return y_cond & __M81_COND_NEG ? 1 / x : x;
  62
  63   if (y == 2)
  64     return x * x;
  65   if (y == 0.5 && !(x_cond & __M81_COND_NEG))
  66     return m81(__ieee754_sqrt) (x);
  67
  68   if (x == 10.0)
  69     {
  70       __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y));
  71       return z;
  72     }
  73   if (x == 2.0)
  74     {
  75       __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y));
  76       return z;
  77     }
  78
  79   ax = s(fabs) (x);
  80   if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1)
  81     {
  82       z = ax;
  83       if (y_cond & __M81_COND_NEG)
  84         z = 1 / z;
  85       if (x_cond & __M81_COND_NEG)
  86         {
  87           if (y != m81(__rint) (y))
  88             {
  89               if (x == -1)
  90                 z = (z - z) / (z - z);
  91             }
  92           else
  93             goto maybe_negate;
  94         }
  95       return z;
  96     }
  97
  98   if (x_cond & __M81_COND_NEG)
  99     {
 100       if (y == m81(__rint) (y))
 101         {
 102           z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x));
 103         maybe_negate:
 104           /* We always use the long double format, since y is already in
 105              this format and rounding won't change the result.  */
 106           {
 107             int32_t exponent;
 108             u_int32_t i0, i1;
 109             GET_LDOUBLE_WORDS (exponent, i0, i1, y);
 110             exponent = (exponent & 0x7fff) - 0x3fff;
 111             if (exponent <= 31
 112                 ? i0 & (1 << (31 - exponent))
 113                 : (exponent <= 63
 114                    && i1 & (1 << (63 - exponent))))
 115               z = -z;
 116           }
 117         }
 118       else
 119         z = (y - y) / (y - y);
 120     }
 121   else
 122     z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x));
 123   return z;
 124 }