|
1 // Copyright (c) 1995-2009 Nokia Corporation and/or its subsidiary(-ies). |
|
2 // All rights reserved. |
|
3 // This component and the accompanying materials are made available |
|
4 // under the terms of the License "Symbian Foundation License v1.0" |
|
5 // which accompanies this distribution, and is available |
|
6 // at the URL "http://www.symbianfoundation.org/legal/sfl-v10.html". |
|
7 // |
|
8 // Initial Contributors: |
|
9 // Nokia Corporation - initial contribution. |
|
10 // |
|
11 // Contributors: |
|
12 // |
|
13 // Description: |
|
14 // e32\euser\maths\um_pow.cpp |
|
15 // Raise to the power. |
|
16 // |
|
17 // |
|
18 |
|
19 #include "um_std.h" |
|
20 |
|
21 #if defined(__USE_VFP_MATH) && !defined(__CPU_HAS_VFP) |
|
22 #error __USE_VFP_MATH was defined but not __CPU_HAS_VFP - impossible combination, check variant.mmh |
|
23 #endif |
|
24 |
|
25 |
|
26 #ifndef __USE_VFP_MATH |
|
27 |
|
28 LOCAL_D const TUint32 ArtanhCoeffs[] = |
|
29 { |
|
30 0x5C17F0BC,0xB8AA3B29,0x80010000, // polynomial approximation to (4/ln2)artanh(x) |
|
31 0xD01FDDD8,0xF6384EE1,0x7FFF0000, // for |x| <= (sqr2-1)/(sqr2+1) |
|
32 0x7D0DDC69,0x93BB6287,0x7FFF0000, |
|
33 0x6564D4F5,0xD30BB153,0x7FFE0000, |
|
34 0x1546C858,0xA4258A33,0x7FFE0000, |
|
35 0xCCE50DA9,0x864D28DF,0x7FFE0000, |
|
36 0x8E1A5DBB,0xE35271A0,0x7FFD0000, |
|
37 0xF5A67D92,0xC3A36B08,0x7FFD0000, |
|
38 0x62D53E02,0xC4A1FFAC,0x7FFD0000 |
|
39 }; |
|
40 |
|
41 LOCAL_D const TUint32 TwoToxCoeffs[] = |
|
42 { |
|
43 0x00000000,0x80000000,0x7FFF0000, // polynomial approximation to 2^(x/8) for |
|
44 0xD1CF79AC,0xB17217F7,0x7FFB0000, // 0<=x<=1 |
|
45 0x162CF72B,0xF5FDEFFC,0x7FF60000, |
|
46 0x23EC0D04,0xE35846B8,0x7FF10000, |
|
47 0xBDB408D7,0x9D955B7E,0x7FEC0000, |
|
48 0xFDD8A678,0xAEC3FE73,0x7FE60000, |
|
49 0xBD6E3950,0xA184E90A,0x7FE00000, |
|
50 0xC1054DA3,0xFFB259D8,0x7FD90000, |
|
51 0x70893DE4,0xB8BEDE2F,0x7FD30000 |
|
52 }; |
|
53 |
|
54 LOCAL_D const TUint32 TwoToNover8[] = |
|
55 { |
|
56 0xEA8BD6E7,0x8B95C1E3,0x7FFF0000, // 2^0.125 |
|
57 0x8DB8A96F,0x9837F051,0x7FFF0000, // 2^0.250 |
|
58 0xB15138EA,0xA5FED6A9,0x7FFF0000, // 2^0.375 |
|
59 0xF9DE6484,0xB504F333,0x7FFF0000, // 2^0.500 |
|
60 0x5506DADD,0xC5672A11,0x7FFF0000, // 2^0.625 |
|
61 0xD69D6AF4,0xD744FCCA,0x7FFF0000, // 2^0.750 |
|
62 0xDD24392F,0xEAC0C6E7,0x7FFF0000 // 2^0.875 |
|
63 }; |
|
64 |
|
65 LOCAL_D const TUint32 Sqr2data[] = {0xF9DE6484,0xB504F333,0x7FFF0000}; // sqr2 |
|
66 LOCAL_D const TUint32 Sqr2Invdata[] = {0xF9DE6484,0xB504F333,0x7FFE0000}; // 1/sqr2 |
|
67 LOCAL_D const TUint32 Onedata[] = {0x00000000,0x80000000,0x7FFF0000}; // 1.0 |
|
68 |
|
69 LOCAL_C void Log2(TRealX& y, TRealX& x) |
|
70 { |
|
71 // Calculate log2(x) and write to y |
|
72 // Result to 64-bit precision to allow accurate powers |
|
73 // Algorithm: |
|
74 // log2(aSrc)=log2(2^e.m) e=exponent of aSrc, m=mantissa 1<=m<2 |
|
75 // log2(aSrc)=e+log2(m) |
|
76 // If e=-1 (0.5<=aSrc<1), let x=aSrc else let x=mantissa(aSrc) |
|
77 // If x>Sqr2, replace x with x/Sqr2 |
|
78 // If x<Sqr2/2, replace x with x*Sqr2 |
|
79 // Replace x with (x-1)/(x+1) |
|
80 // Use polynomial to calculate artanh(x) for |x| <= (sqr2-1)/(sqr2+1) |
|
81 // ( use identity ln(x) = 2artanh((x-1)/(x+1)) ) |
|
82 |
|
83 const TRealX& Sqr2=*(const TRealX*)Sqr2data; |
|
84 const TRealX& Sqr2Inv=*(const TRealX*)Sqr2Invdata; |
|
85 const TRealX& One=*(const TRealX*)Onedata; |
|
86 |
|
87 TInt n=(x.iExp-0x7FFF)<<1; |
|
88 x.iExp=0x7FFF; |
|
89 if (n!=-2) |
|
90 { |
|
91 if (x>Sqr2) |
|
92 { |
|
93 x*=Sqr2Inv; |
|
94 n++; |
|
95 } |
|
96 } |
|
97 else |
|
98 { |
|
99 n=0; |
|
100 x.iExp=0x7FFE; |
|
101 if (x<Sqr2Inv) |
|
102 { |
|
103 x*=Sqr2; |
|
104 n--; |
|
105 } |
|
106 } |
|
107 x=(x-One)/(x+One); // ln(x)=2artanh((x-1)/(x+1)) |
|
108 Math::PolyX(y,x*x,8,(const TRealX*)ArtanhCoeffs); |
|
109 y*=x; |
|
110 y+=TRealX(n); |
|
111 if (y.iExp>1) |
|
112 y.iExp--; |
|
113 else |
|
114 y.iExp=0; |
|
115 } |
|
116 |
|
117 LOCAL_C TInt TwoTox(TRealX& y, TRealX& x) |
|
118 { |
|
119 // Calculate 2^x and write result to y. Result to 64 bit precision. |
|
120 // Algorithm: |
|
121 // 2^x = 2^int(x).2^frac(x) |
|
122 // 2^int(x) just adds int(x) to the final result exponent |
|
123 // Reduce frac(x) to the range [0,0.125] (modulo 0.125) |
|
124 // Use polynomial to calculate 2^x for 0<=x<=0.125 |
|
125 // Multiply by 2^(n/8) for n=0,1,2,3,4,5,6,7 to give 2^frac(x) |
|
126 |
|
127 if (x.iExp) |
|
128 x.iExp+=3; |
|
129 TInt n=(TInt)x; |
|
130 if (n<16384 && n>-16384) |
|
131 { |
|
132 if (x.iSign&1) |
|
133 n--; |
|
134 x-=TRealX(n); |
|
135 Math::PolyX(y,x,8,(const TRealX*)TwoToxCoeffs); |
|
136 y.iExp=TUint16(TInt(y.iExp)+(n>>3)); |
|
137 n&=7; |
|
138 if (n) |
|
139 y*= (*(const TRealX*)(TwoToNover8+3*n-3)); |
|
140 return KErrNone; |
|
141 } |
|
142 else |
|
143 { |
|
144 if (n<0) |
|
145 { |
|
146 y.SetZero(); |
|
147 return KErrUnderflow; |
|
148 } |
|
149 else |
|
150 { |
|
151 y.SetInfinite(0); |
|
152 return KErrOverflow; |
|
153 } |
|
154 } |
|
155 } |
|
156 |
|
157 |
|
158 |
|
159 |
|
160 EXPORT_C TInt Math::Pow(TReal &aTrg,const TReal &aSrc,const TReal &aPower) |
|
161 /** |
|
162 Calculates the value of x raised to the power of y. |
|
163 |
|
164 The behaviour conforms to that specified for pow() in the |
|
165 ISO C Standard ISO/IEC 9899 (Annex F), although floating-point exceptions |
|
166 are not supported. |
|
167 |
|
168 @param aTrg A reference containing the result. |
|
169 @param aSrc The x argument of the function. |
|
170 @param aPower The y argument of the function. |
|
171 |
|
172 @return KErrNone if successful; |
|
173 KErrOverflow if the result is +/- infinity; |
|
174 KErrUnderflow if the result is too small to be represented; |
|
175 KErrArgument if the result is not a number (NaN). |
|
176 */ |
|
177 // |
|
178 // Evaluates aSrc raised to the power aPower and places the result in aTrg. |
|
179 // For non-special values algorithm is aTrg=2^(aPower*log2(aSrc)) |
|
180 // |
|
181 { |
|
182 TRealX x,p; |
|
183 |
|
184 TInt ret2=p.Set(aPower); |
|
185 // pow(x, +/-0) -> 1 for any x, even a NaN |
|
186 if (p.IsZero()) |
|
187 { |
|
188 aTrg=1.0; |
|
189 return KErrNone; |
|
190 } |
|
191 |
|
192 TInt ret1=x.Set(aSrc); |
|
193 if (ret1==KErrArgument || ret2==KErrArgument) |
|
194 { |
|
195 // pow(+1, y) -> 1 for any y, even a NaN |
|
196 // XXX First test should not be necessary, but on WINS |
|
197 // aSrc == 1.0 is true when aSrc is NaN. |
|
198 if (ret1 != KErrArgument && aSrc == 1.0) |
|
199 { |
|
200 aTrg=aSrc; |
|
201 return KErrNone; |
|
202 } |
|
203 SetNaN(aTrg); |
|
204 return KErrArgument; |
|
205 } |
|
206 |
|
207 // Infinite power |
|
208 if (ret2==KErrOverflow) |
|
209 { |
|
210 // figure out which of these cases we have: |
|
211 // |
|
212 // pow(x, -INF) -> +INF for |x| < 1 } flag = 0 |
|
213 // pow(x, +INF) -> +INF for |x| > 1 } |
|
214 // pow(x, -INF) -> +0 for |x| > 1 } flag = 1 |
|
215 // pow(x, +INF) -> +0 for |x| < 1 } |
|
216 // |
|
217 // flag = 2 => |x| == 1.0 |
|
218 // |
|
219 TInt flag=2; |
|
220 if (Abs(aSrc)>1.0) |
|
221 flag=p.iSign&1; |
|
222 if (Abs(aSrc)<1.0) |
|
223 flag=1-(p.iSign&1); |
|
224 if (flag==0) |
|
225 { |
|
226 SetInfinite(aTrg,0); |
|
227 return KErrOverflow; |
|
228 } |
|
229 if (flag==1) |
|
230 { |
|
231 SetZero(aTrg,0); |
|
232 return KErrNone; |
|
233 } |
|
234 if (Abs(aSrc)==1.0) |
|
235 { |
|
236 // pow(-1, +/-INF) -> 1 |
|
237 aTrg=1.0; |
|
238 return KErrNone; |
|
239 } |
|
240 // This should never happen (i.e. aSrc is NaN, which |
|
241 // should be taken care of above) |
|
242 SetNaN(aTrg); |
|
243 return KErrArgument; |
|
244 } |
|
245 |
|
246 // Negative Base raised to a power |
|
247 TInt odd=1; |
|
248 if (x.iSign & 1) |
|
249 { |
|
250 TReal pint; |
|
251 Math::Int(pint,aPower); |
|
252 if (aPower-pint) // Checks that if aSrc is less than zero, then aPower is integral |
|
253 { |
|
254 // pow(-INF, y) -> +0 for y < 0 and not an odd integer |
|
255 // pow(-INF, y) -> +INF for y > 0 and not an odd integer |
|
256 // Since we're here, aPower is not integral, so can't be odd, either |
|
257 if (ret1 == KErrOverflow) |
|
258 { |
|
259 if (aPower < 0) |
|
260 { |
|
261 SetZero(aTrg); |
|
262 return KErrNone; |
|
263 } |
|
264 else |
|
265 { |
|
266 SetInfinite(aTrg,0); |
|
267 return KErrOverflow; |
|
268 } |
|
269 } |
|
270 SetNaN(aTrg); |
|
271 return KErrArgument; |
|
272 } |
|
273 TReal powerby2=aPower*0.5; |
|
274 Math::Int(pint,powerby2); |
|
275 if (powerby2-pint) |
|
276 odd=(-1); |
|
277 x.iSign=0; |
|
278 } |
|
279 |
|
280 // Zero or infinity raised to a power |
|
281 if (x.IsZero() || ret1==KErrOverflow) |
|
282 { |
|
283 if (x.IsZero() && p.IsZero()) |
|
284 { |
|
285 aTrg=1.0; |
|
286 return KErrNone; |
|
287 } |
|
288 TInt sign=(odd==-1 ? 1 : 0); |
|
289 if ((x.IsZero() && (p.iSign&1)==0) || (ret1==KErrOverflow && (p.iSign&1))) |
|
290 { |
|
291 SetZero(aTrg,sign); |
|
292 return KErrNone; |
|
293 } |
|
294 else |
|
295 { |
|
296 SetInfinite(aTrg,sign); |
|
297 return KErrOverflow; |
|
298 } |
|
299 } |
|
300 |
|
301 TRealX y; |
|
302 Log2(y,x); |
|
303 x=y*p; // this cannot overflow or underflow |
|
304 TInt r=TwoTox(y,x); |
|
305 if (odd<0) |
|
306 y.iSign=1; |
|
307 TInt r2=y.GetTReal(aTrg); |
|
308 return (r==KErrNone)?r2:r; |
|
309 } |
|
310 |
|
311 #else // __USE_VFP_MATH |
|
312 |
|
313 // definitions come from RVCT math library |
|
314 extern "C" TReal pow(TReal,TReal); |
|
315 |
|
316 EXPORT_C TInt Math::Pow(TReal &aTrg,const TReal &aSrc,const TReal &aPower) |
|
317 { |
|
318 aTrg = pow(aSrc,aPower); |
|
319 if (Math::IsZero(aTrg) && !Math::IsZero(aSrc) && !Math::IsInfinite(aSrc) && !Math::IsInfinite(aPower)) |
|
320 return KErrUnderflow; |
|
321 if (Math::IsFinite(aTrg)) |
|
322 return KErrNone; |
|
323 if (Math::IsZero(aPower)) // pow(x, +/-0) -> 1 for any x, even a NaN |
|
324 { |
|
325 aTrg = 1.0; |
|
326 return KErrNone; |
|
327 } |
|
328 if (Math::IsInfinite(aTrg)) |
|
329 return KErrOverflow; |
|
330 if (aSrc==1.0) // pow(+1, y) -> 1 for any y, even a NaN |
|
331 { |
|
332 aTrg=aSrc; |
|
333 return KErrNone; |
|
334 } |
|
335 if (Math::IsInfinite(aPower)) |
|
336 { |
|
337 if (aSrc == -1.0) // pow(-1, +/-INF) -> 1 |
|
338 { |
|
339 aTrg = 1.0; |
|
340 return KErrNone; |
|
341 } |
|
342 if (((Abs(aSrc) < 1) && (aPower < 0)) || // pow(x, -INF) -> +INF for |x| < 1 |
|
343 ((Abs(aSrc) > 1) && (aPower > 0))) // pow(x, +INF) -> +INF for |x| > 1 |
|
344 { |
|
345 SetInfinite(aTrg,0); |
|
346 return KErrOverflow; |
|
347 } |
|
348 } |
|
349 // pow(-INF, y) -> +INF for y > 0 and not an odd integer |
|
350 if (Math::IsInfinite(aSrc) && (aSrc < 0) && (aPower > 0)) |
|
351 { |
|
352 TBool odd = EFalse; |
|
353 TReal pint; |
|
354 Math::Int(pint, aPower); |
|
355 if (aPower == pint) |
|
356 { |
|
357 TReal halfPower = aPower * 0.5; |
|
358 Math::Int(pint, halfPower); |
|
359 if (halfPower != pint) |
|
360 odd = ETrue; |
|
361 } |
|
362 if (odd == EFalse) |
|
363 { |
|
364 SetInfinite(aTrg,0); |
|
365 return KErrOverflow; |
|
366 } |
|
367 } |
|
368 |
|
369 // Otherwise... |
|
370 SetNaN(aTrg); |
|
371 return KErrArgument; |
|
372 } |
|
373 |
|
374 #endif |