source: trunk/libs/newlib/src/newlib/libm/mathfp/sf_pow.c @ 503

Last change on this file since 503 was 444, checked in by satin@…, 6 years ago

add newlib,libalmos-mkh, restructure shared_syscalls.h and mini-libc
File size: 2.4 KB
Line 
1
2/* @(#)z_powf.c 1.0 98/08/13 */
3#include <float.h>
4#include "fdlibm.h"
5#include "zmath.h"
6
7float powf (float x, float y)
8{
9  float d, k, t, r = 1.0;
10  int n, sign, exponent_is_even_int = 0;
11  __int32_t px;
12
13  GET_FLOAT_WORD (px, x);
14
15  k = modff (y, &d);
16
17  if (k == 0.0) 
18    {
19      /* Exponent y is an integer. */
20      if (modff (ldexpf (y, -1), &t))
21        {
22          /* y is odd. */
23          exponent_is_even_int = 0;
24        }
25      else
26        {
27          /* y is even. */
28          exponent_is_even_int = 1;
29        }
30    }
31
32  if (x == 0.0)
33    {
34      if (y <= 0.0)
35        errno = EDOM;
36    }
37  else if ((t = y * log (fabsf (x))) >= BIGX) 
38    {
39      errno = ERANGE;
40      if (px & 0x80000000) 
41        {
42          /* x is negative. */
43          if (k) 
44            {
45              /* y is not an integer. */
46              errno = EDOM;
47              x = 0.0;
48            }
49          else if (exponent_is_even_int)
50            x = z_infinity_f.f;
51          else
52            x = -z_infinity_f.f;
53        }
54    else
55      {
56        x = z_infinity_f.f;
57      }
58    }
59  else if (t < SMALLX)
60    {
61      errno = ERANGE;
62      x = 0.0;
63    }
64  else 
65    {
66      if ( !k && fabsf (d) <= 32767 ) 
67        {
68          n = (int) d;
69
70          if ((sign = (n < 0)))
71            n = -n;
72
73          while ( n > 0 ) 
74            {
75              if ((unsigned int) n % 2) 
76                r *= x;
77              x *= x;
78              n = (unsigned int) n / 2;
79            }
80
81          if (sign)
82            r = 1.0 / r;
83
84          return r;
85        }
86      else 
87        {
88          if ( px & 0x80000000 ) 
89            {
90              /* x is negative. */
91              if (k) 
92                {
93                  /* y is not an integer. */
94                  errno = EDOM;
95                  return 0.0;
96                }
97            }
98
99          x = exp (t);
100
101          if (!exponent_is_even_int) 
102            { 
103              if (px & 0x80000000)
104                {
105                  /* y is an odd integer, and x is negative,
106                     so the result is negative. */
107                  GET_FLOAT_WORD (px, x);
108                  px |= 0x80000000;
109                  SET_FLOAT_WORD (x, px);
110                }
111            }
112        }
113    }
114
115  return x;
116}
117
118#ifdef _DOUBLE_IS_32BITS
119
120double pow (double x, double y)
121{
122  return (double) powf ((float) x, (float) y);
123}
124
125#endif /* defined(_DOUBLE_IS_32BITS) */ 
Note: See TracBrowser for help on using the repository browser.