source: trunk/libs/newlib/src/newlib/libm/mathfp/s_pow.c @ 471

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

add newlib,libalmos-mkh, restructure shared_syscalls.h and mini-libc
File size: 3.2 KB
Line 
1
2/* @(#)z_pow.c 1.0 98/08/13 */
3
4/*
5FUNCTION
6        <<pow>>, <<powf>>---x to the power y
7INDEX
8        pow
9INDEX
10        powf
11
12
13SYNOPSIS
14        #include <math.h>
15        double pow(double <[x]>, double <[y]>);
16        float pow(float <[x]>, float <[y]>);
17
18DESCRIPTION
19        <<pow>> and <<powf>> calculate <[x]> raised to the exponent <[y]>.
20        @tex
21        (That is, $x^y$.)
22        @end tex
23
24RETURNS
25        On success, <<pow>> and <<powf>> return the value calculated.
26
27        When the argument values would produce overflow, <<pow>>
28        returns <<HUGE_VAL>> and set <<errno>> to <<ERANGE>>.  If the
29        argument <[x]> passed to <<pow>> or <<powf>> is a negative
30        noninteger, and <[y]> is also not an integer, then <<errno>>
31        is set to <<EDOM>>.  If <[x]> and <[y]> are both 0, then
32        <<pow>> and <<powf>> return <<1>>.
33
34        You can modify error handling for these functions using <<matherr>>.
35
36PORTABILITY
37        <<pow>> is ANSI C. <<powf>> is an extension.  */
38
39#include <float.h>
40#include "fdlibm.h"
41#include "zmath.h"
42
43#ifndef _DOUBLE_IS_32BITS
44
45double pow (double x, double y)
46{
47  double d, k, t, r = 1.0;
48  int n, sign, exponent_is_even_int = 0;
49  __uint32_t px;
50
51  GET_HIGH_WORD (px, x);
52
53  k = modf (y, &d);
54
55  if (k == 0.0)
56    {
57      /* Exponent y is an integer. */
58      if (modf (ldexp (y, -1), &t))
59        {
60          /* y is odd. */
61          exponent_is_even_int = 0;
62        }
63      else
64        {
65          /* y is even. */
66          exponent_is_even_int = 1;
67        }
68    }
69
70  if (x == 0.0)
71    {
72      if (y <= 0.0)
73        errno = EDOM;
74    }
75  else if ((t = y * log (fabs (x))) >= BIGX) 
76    {
77      errno = ERANGE;
78      if (px & 0x80000000) 
79        {
80          /* x is negative. */
81          if (k) 
82            {
83              /* y is not an integer. */
84              errno = EDOM;
85              x = 0.0;
86            }
87          else if (exponent_is_even_int)
88            x = z_infinity.d;
89          else
90            x = -z_infinity.d;
91        }
92      else
93        {
94          x = z_infinity.d;
95        }
96    }
97  else if (t < SMALLX)
98    {
99      errno = ERANGE;
100      x = 0.0;
101    }
102  else 
103    {
104      if ( !k && fabs(d) <= 32767 ) 
105        {
106          n = (int) d;
107
108          if ((sign = (n < 0)))
109            n = -n;
110
111          while ( n > 0 ) 
112            {
113              if ((unsigned int) n % 2) 
114                r *= x;
115              x *= x;
116              n = (unsigned int) n / 2;
117            }
118
119          if (sign)
120            r = 1.0 / r;
121
122          return r;
123        }
124      else 
125        {
126          if ( px & 0x80000000 ) 
127            {
128              /* x is negative. */
129              if ( k ) 
130                {
131                  /* y is not an integer. */
132                  errno = EDOM;
133                  return 0.0;
134                }
135            }
136
137          x = exp (t);
138
139          if (!exponent_is_even_int)
140            {
141              if (px & 0x80000000)
142                {
143                  /* y is an odd integer, and x is negative,
144                     so the result is negative. */
145                  GET_HIGH_WORD (px, x);
146                  px |= 0x80000000;
147                  SET_HIGH_WORD (x, px);
148                }
149            }
150        }
151    }
152
153  return x;
154}
155
156#endif _DOUBLE_IS_32BITS
Note: See TracBrowser for help on using the repository browser.