sysdeps/m68k/fpu/e_pow.c

   1 /* Copyright (C) 1997, 1999 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, write to the Free
  16    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  17    02111-1307 USA.  */
  18
  19 #include <math.h>
  20 #include "math_private.h"
  21 #include "mathimpl.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   if (y_cond & __M81_COND_NAN)
  45     return x == 1.0 ? x : x + y;
  46
  47   x_cond = __m81_test (x);
  48   if (x_cond & __M81_COND_NAN)
  49     return x + y;
  50
  51   if (y_cond & __M81_COND_INF)
  52     {
  53       ax = s(fabs) (x);
  54       if (ax == 1.0)
  55         return ax;
  56       if (ax > 1.0)
  57         return y_cond & __M81_COND_NEG ? 0 : y;
  58       else
  59         return y_cond & __M81_COND_NEG ? -y : 0;
  60     }
  61
  62   if (s(fabs) (y) == 1.0)
  63     return y_cond & __M81_COND_NEG ? 1 / x : x;
  64
  65   if (y == 2)
  66     return x * x;
  67   if (y == 0.5 && !(x_cond & __M81_COND_NEG))
  68     return m81(__ieee754_sqrt) (x);
  69
  70   if (x == 10.0)
  71     {
  72       __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y));
  73       return z;
  74     }
  75   if (x == 2.0)
  76     {
  77       __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y));
  78       return z;
  79     }
  80
  81   ax = s(fabs) (x);
  82   if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1.0)
  83     {
  84       z = ax;
  85       if (y_cond & __M81_COND_NEG)
  86         z = 1 / z;
  87       if (x_cond & __M81_COND_NEG)
  88         {
  89           if (y != m81(__rint) (y))
  90             {
  91               if (x == -1)
  92                 z = (z - z) / (z - z);
  93             }
  94           else
  95             goto maybe_negate;
  96         }
  97       return z;
  98     }
  99
 100   if (x_cond & __M81_COND_NEG)
 101     {
 102       if (y == m81(__rint) (y))
 103         {
 104           z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x));
 105         maybe_negate:
 106           /* We always use the long double format, since y is already in
 107              this format and rounding won't change the result.  */
 108           {
 109             int32_t exponent;
 110             u_int32_t i0, i1;
 111             GET_LDOUBLE_WORDS (exponent, i0, i1, y);
 112             exponent = (exponent & 0x7fff) - 0x3fff;
 113             if (exponent <= 31
 114                 ? i0 & (1 << (31 - exponent))
 115                 : (exponent <= 63
 116                    && i1 & (1 << (63 - exponent))))
 117               z = -z;
 118           }
 119         }
 120       else
 121         z = (y - y) / (y - y);
 122     }
 123   else
 124     z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x));
 125   return z;
 126 }