[444] | 1 | |
---|
| 2 | /* @(#)z_pow.c 1.0 98/08/13 */ |
---|
| 3 | |
---|
| 4 | /* |
---|
| 5 | FUNCTION |
---|
| 6 | <<pow>>, <<powf>>---x to the power y |
---|
| 7 | INDEX |
---|
| 8 | pow |
---|
| 9 | INDEX |
---|
| 10 | powf |
---|
| 11 | |
---|
| 12 | |
---|
| 13 | SYNOPSIS |
---|
| 14 | #include <math.h> |
---|
| 15 | double pow(double <[x]>, double <[y]>); |
---|
| 16 | float pow(float <[x]>, float <[y]>); |
---|
| 17 | |
---|
| 18 | DESCRIPTION |
---|
| 19 | <<pow>> and <<powf>> calculate <[x]> raised to the exponent <[y]>. |
---|
| 20 | @tex |
---|
| 21 | (That is, $x^y$.) |
---|
| 22 | @end tex |
---|
| 23 | |
---|
| 24 | RETURNS |
---|
| 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 | |
---|
| 36 | PORTABILITY |
---|
| 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 | |
---|
| 45 | double 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 |
---|