sysdeps/m68k/m680x0/fpu/e_pow.c

   1 /* Copyright (C) 1997-2014 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 Lesser General Public
   6    License as published by the Free Software Foundation; either
   7    version 2.1 of the 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    Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public
  15    License along with the GNU C Library.  If not, see
  16    <http://www.gnu.org/licenses/>.  */
  17
  18 #include <math.h>
  19 #include <math_private.h>
  20 #include "mathimpl.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   if (y_cond & __M81_COND_NAN)
  44     return x == 1.0 ? x : x + y;
  45
  46   x_cond = __m81_test (x);
  47   if (x_cond & __M81_COND_NAN)
  48     return x + y;
  49
  50   if (y_cond & __M81_COND_INF)
  51     {
  52       ax = s(fabs) (x);
  53       if (ax == 1.0)
  54         return ax;
  55       if (ax > 1.0)
  56         return y_cond & __M81_COND_NEG ? 0 : y;
  57       else
  58         return y_cond & __M81_COND_NEG ? -y : 0;
  59     }
  60
  61   if (s(fabs) (y) == 1.0)
  62     return y_cond & __M81_COND_NEG ? 1 / x : x;
  63
  64   if (y == 2)
  65     return x * x;
  66   if (y == 0.5 && !(x_cond & __M81_COND_NEG))
  67     return m81(__ieee754_sqrt) (x);
  68
  69   if (x == 10.0)
  70     {
  71       __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y));
  72       return z;
  73     }
  74   if (x == 2.0)
  75     {
  76       __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y));
  77       return z;
  78     }
  79
  80   ax = s(fabs) (x);
  81   if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1.0)
  82     {
  83       z = ax;
  84       if (y_cond & __M81_COND_NEG)
  85         z = 1 / z;
  86       if (x_cond & __M81_COND_NEG)
  87         {
  88           if (y != m81(__rint) (y))
  89             {
  90               if (x == -1)
  91                 z = (z - z) / (z - z);
  92             }
  93           else
  94             goto maybe_negate;
  95         }
  96       return z;
  97     }
  98
  99   if (x_cond & __M81_COND_NEG)
 100     {
 101       if (y == m81(__rint) (y))
 102         {
 103           z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x));
 104         maybe_negate:
 105           /* We always use the long double format, since y is already in
 106              this format and rounding won't change the result.  */
 107           {
 108             int32_t exponent;
 109             u_int32_t i0, i1;
 110             GET_LDOUBLE_WORDS (exponent, i0, i1, y);
 111             exponent = (exponent & 0x7fff) - 0x3fff;
 112             if (exponent <= 31
 113                 ? i0 & (1 << (31 - exponent))
 114                 : (exponent <= 63
 115                    && i1 & (1 << (63 - exponent))))
 116               z = -z;
 117           }
 118         }
 119       else
 120         z = (y - y) / (y - y);
 121     }
 122   else
 123     z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x));
 124   return z;
 125 }
 126 strong_alias (s(__ieee754_pow), CONCATX (s(__pow), _finite))