overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// PowerDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// long double powl( long double base, long double exponent );
29//
30
31#include "math.h"
32#include "fp_private.h"
33#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
34#include "DD.h"
35
36extern long double rintl(long double);
37
38
39// Floating-point constants
40
41static const Double kInfinityu = {{0x7ff00000,0x0}};
42
43long double _ExpInnerLD(double, double, double, double *, int);
44long double _LogInnerLD(double, double, double, double *, int);
45static long double PowerInnerLD(double, double, double, double);
46
47long double powl( long double base, long double exponent )
48{
49 DblDbl u, v;
50 double basehi, baselo, exphi, explo;
51 int isExpOddInt;
52 uint32_t expExp;
53 double fpenv;
54 long double ldTemp;
55
56 u.f = base;
57 basehi = u.d[0];
58 baselo = u.d[1];
59 v.f = exponent;
60 exphi = v.d[0];
61 explo = v.d[1];
62
63 if (((v.i[0] & 0x000fffffu) | v.i[1] | (v.i[2] & 0x7fffffffu) | v.i[3]) == 0) {
64 expExp = v.i[0] & 0xfff00000u;
65 if (expExp == 0x40000000) return base*base; // if exponent == 2.0
66 else if (exponent == 0.0) return 1.0; // if exponent == 0.0
67 else if (expExp == 0x3ff00000) return base; // if exponent == 1.0
68 else if (expExp == 0xbff00000) return 1.0/base; // if exponent == -1.0
69 }
70
71 FEGETENVD(fpenv);
72 FESETENVD(0.0);
73
74 if ((v.i[0] & 0x7ff00000u) < 0x7ff00000u) { // exponent is finite
75 if ((u.i[0] & 0x7ff00000u) < 0x7ff00000u) { // base is finite
76 if (basehi > 0.0) { // base positive
77 ldTemp = PowerInnerLD(basehi, baselo, exphi, explo);
78 FESETENVD(fpenv);
79 return ldTemp;
80 }
81 if (basehi < 0.0) { // base negative
82 if (rintl(exponent) != exponent) { // exponent is non-integer
83 u.d[0] = NAN;
84 u.d[1] = 0.0;
85 FESETENVD(fpenv);
86 return u.f;
87 }
88 u.f = PowerInnerLD(-basehi, -baselo, exphi, explo);
89 if (rintl(0.5*exponent) != 0.5*exponent) // exponent is odd
90 u.f = -u.f;
91 FESETENVD(fpenv);
92 return u.f;
93 }
94 else { // base is 0.0:
95 isExpOddInt = ((rintl(exponent) == exponent) ?
96 (rintl(0.5*exponent) != 0.5*exponent) : 0);
97 if (exphi > 0.0) {
98 ldTemp = (isExpOddInt) ? base : 0.0;
99 FESETENVD(fpenv);
100 return ldTemp;
101 }
102 else { // exponent < 0.0
103 u.d[0] = (isExpOddInt) ? 1.0/basehi : 1.0/fabs(basehi);
104 u.d[1] = 0.0;
105 FESETENVD(fpenv);
106 return u.f;
107 }
108 }
109 } // Exponent is finite, Base is not.
110 else if (basehi != basehi) // Base is NaN, return NaN
111 return base;
112 else // base is +-Inf
113 {
114 if (basehi > 0.0) // base is +inf
115 {
116 if (exphi > 0.0) // pow(+inf, >0)
117 return (long double) INFINITY;
118 else // pow(+inf, <0)
119 return 0.0L;
120 }
121 else // base is -inf
122 {
123 // If exponent is an integer, and exponent/2 is not, then exponent is an odd integer.
124 isExpOddInt = ((rintl(exponent) == exponent) && (rintl(0.5*exponent) != 0.5*exponent));
125 if (exphi > 0.0) {
126 // previously signs were incorrect in this expression; we should be returning
127 // -Inf when the exponent is an odd integer, +Inf otherwise.
128 ldTemp = (isExpOddInt) ? (long double) -INFINITY : (long double) INFINITY;
129 FESETENVD(fpenv);
130 return ldTemp;
131 }
132 else { // exponent < 0.0
133 // Return -0 if exponent is an odd integer, +0 otherwise.
134 u.d[0] = (isExpOddInt) ? -0.0 : 0.0;
135 u.d[1] = 0.0;
136 FESETENVD(fpenv);
137 return u.f;
138 }
139 }
140 }
141 }
142 else if (exphi != exphi) // Exponent is a NaN
143 {
144 if (base == 1.0L)
145 return base;
146 else
147 return base + exponent;
148 }
149 else // exponent is +-Inf
150 {
151 if (basehi != basehi)
152 return base;
153
154 ldTemp = fabsl(base);
155 if ( (exphi > 0.0 && ldTemp > 1.0L) || (exphi < 0.0 && ldTemp < 1.0L) )
156 return (long double) INFINITY;
157 else if ( (exphi > 0.0 && ldTemp < 1.0L) || (exphi < 0.0 && ldTemp > 1.0L) )
158 return 0.0L;
159 else
160 return 1.0L;
161 }
162}
163
164//***************************************************************************
165//
166// FUNCTION: static long double PowerInnerLD(double alpha, double beta,
167// double beta, double gamma);
168//
169// DESCRIPTION: Computes the base to the exponent power where: alpha/beta
170// are the head/tail of base, and gamma/delta are head/tail
171// of exponent . This routine is called internally by the
172// long double Power function family. It assumes that
173// base is finite and > 0, exponent is finite.
174//
175//***************************************************************************/
176
177static long double PowerInnerLD(double alpha, double beta, double gamma, double delta)
178{
179 DblDbl u;
180 double extra, ahi, amid, amid2, alow, amid3;
181 double amid4, ahi1, amid5, amid6, ahi2, amid7, amid8;
182
183 u.f = _LogInnerLD(alpha, beta, 0.0, &extra, 1);
184 ahi = __FMUL( u.d[0], gamma );
185 amid = __FMSUB( u.d[0], gamma, ahi ); // u.d[0] * gamma - ahi;
186 amid2 = __FMUL( u.d[1], gamma );
187 alow = __FMADD( gamma, extra, __FMSUB( u.d[1], gamma, amid2 ) ); // u.d[1] * gamma - amid2 + gamma * extra;
188 amid3 = __FMUL( u.d[0], delta );
189 alow = __FMADD( delta, u.d[1], __FMSUB( u.d[0], delta, amid3) + alow ); // u.d[0] * delta - amid3 + alow + delta * u.d[1];
190 amid4 = __FADD( amid, amid2 );
191 ahi1 = __FADD( ahi, amid4 );
192 amid5 = ahi - ahi1 + amid4;
193 if (fabs(amid) > fabs(amid2))
194 alow = amid - amid4 + amid2 + alow;
195 else
196 alow = amid2 - amid4 + amid + alow;
197 amid6 = __FADD( amid3, amid5 );
198 ahi2 = __FADD( ahi1, amid6 );
199 if (fabs(amid3) > fabs(amid5))
200 alow = amid3 - amid6 + amid5 + alow;
201 else
202 alow = amid5 - amid6 + amid3 + alow;
203 amid7 = ahi1 - ahi2 + amid6;
204 amid8 = __FADD( amid7, alow );
205 alow = amid7 - amid8 + alow;
206 if (fabs(ahi) == kInfinityu.f) {
207 if (ahi > 0.0) {
208 u.d[0] = kInfinityu.f;
209 u.d[1] = 0.0;
210 }
211 else {
212 u.d[0] = 0.0;
213 u.d[1] = 0.0;
214 }
215 }
216 else {
217 u.f = _ExpInnerLD(ahi2, amid8, alow, &extra, 3);
218 if ((u.i[0] & 0x7ff00000) == 0x7ff00000)
219 u.d[1] = 0.0;
220 }
221 return u.f;
222}
223#endif