diff -r c55016431358 -r 0a7b44b10206 symport/e32/euser/maths/um_pow.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/symport/e32/euser/maths/um_pow.cpp Thu Jun 25 15:59:54 2009 +0100 @@ -0,0 +1,374 @@ +// Copyright (c) 1995-2009 Nokia Corporation and/or its subsidiary(-ies). +// All rights reserved. +// This component and the accompanying materials are made available +// under the terms of the License "Symbian Foundation License v1.0" +// which accompanies this distribution, and is available +// at the URL "http://www.symbianfoundation.org/legal/sfl-v10.html". +// +// Initial Contributors: +// Nokia Corporation - initial contribution. +// +// Contributors: +// +// Description: +// e32\euser\maths\um_pow.cpp +// Raise to the power. +// +// + +#include "um_std.h" + +#if defined(__USE_VFP_MATH) && !defined(__CPU_HAS_VFP) +#error __USE_VFP_MATH was defined but not __CPU_HAS_VFP - impossible combination, check variant.mmh +#endif + + +#ifndef __USE_VFP_MATH + +LOCAL_D const TUint32 ArtanhCoeffs[] = + { + 0x5C17F0BC,0xB8AA3B29,0x80010000, // polynomial approximation to (4/ln2)artanh(x) + 0xD01FDDD8,0xF6384EE1,0x7FFF0000, // for |x| <= (sqr2-1)/(sqr2+1) + 0x7D0DDC69,0x93BB6287,0x7FFF0000, + 0x6564D4F5,0xD30BB153,0x7FFE0000, + 0x1546C858,0xA4258A33,0x7FFE0000, + 0xCCE50DA9,0x864D28DF,0x7FFE0000, + 0x8E1A5DBB,0xE35271A0,0x7FFD0000, + 0xF5A67D92,0xC3A36B08,0x7FFD0000, + 0x62D53E02,0xC4A1FFAC,0x7FFD0000 + }; + +LOCAL_D const TUint32 TwoToxCoeffs[] = + { + 0x00000000,0x80000000,0x7FFF0000, // polynomial approximation to 2^(x/8) for + 0xD1CF79AC,0xB17217F7,0x7FFB0000, // 0<=x<=1 + 0x162CF72B,0xF5FDEFFC,0x7FF60000, + 0x23EC0D04,0xE35846B8,0x7FF10000, + 0xBDB408D7,0x9D955B7E,0x7FEC0000, + 0xFDD8A678,0xAEC3FE73,0x7FE60000, + 0xBD6E3950,0xA184E90A,0x7FE00000, + 0xC1054DA3,0xFFB259D8,0x7FD90000, + 0x70893DE4,0xB8BEDE2F,0x7FD30000 + }; + +LOCAL_D const TUint32 TwoToNover8[] = + { + 0xEA8BD6E7,0x8B95C1E3,0x7FFF0000, // 2^0.125 + 0x8DB8A96F,0x9837F051,0x7FFF0000, // 2^0.250 + 0xB15138EA,0xA5FED6A9,0x7FFF0000, // 2^0.375 + 0xF9DE6484,0xB504F333,0x7FFF0000, // 2^0.500 + 0x5506DADD,0xC5672A11,0x7FFF0000, // 2^0.625 + 0xD69D6AF4,0xD744FCCA,0x7FFF0000, // 2^0.750 + 0xDD24392F,0xEAC0C6E7,0x7FFF0000 // 2^0.875 + }; + +LOCAL_D const TUint32 Sqr2data[] = {0xF9DE6484,0xB504F333,0x7FFF0000}; // sqr2 +LOCAL_D const TUint32 Sqr2Invdata[] = {0xF9DE6484,0xB504F333,0x7FFE0000}; // 1/sqr2 +LOCAL_D const TUint32 Onedata[] = {0x00000000,0x80000000,0x7FFF0000}; // 1.0 + +LOCAL_C void Log2(TRealX& y, TRealX& x) + { + // Calculate log2(x) and write to y + // Result to 64-bit precision to allow accurate powers + // Algorithm: + // log2(aSrc)=log2(2^e.m) e=exponent of aSrc, m=mantissa 1<=m<2 + // log2(aSrc)=e+log2(m) + // If e=-1 (0.5<=aSrc<1), let x=aSrc else let x=mantissa(aSrc) + // If x>Sqr2, replace x with x/Sqr2 + // If xSqr2) + { + x*=Sqr2Inv; + n++; + } + } + else + { + n=0; + x.iExp=0x7FFE; + if (x1) + y.iExp--; + else + y.iExp=0; + } + +LOCAL_C TInt TwoTox(TRealX& y, TRealX& x) + { + // Calculate 2^x and write result to y. Result to 64 bit precision. + // Algorithm: + // 2^x = 2^int(x).2^frac(x) + // 2^int(x) just adds int(x) to the final result exponent + // Reduce frac(x) to the range [0,0.125] (modulo 0.125) + // Use polynomial to calculate 2^x for 0<=x<=0.125 + // Multiply by 2^(n/8) for n=0,1,2,3,4,5,6,7 to give 2^frac(x) + + if (x.iExp) + x.iExp+=3; + TInt n=(TInt)x; + if (n<16384 && n>-16384) + { + if (x.iSign&1) + n--; + x-=TRealX(n); + Math::PolyX(y,x,8,(const TRealX*)TwoToxCoeffs); + y.iExp=TUint16(TInt(y.iExp)+(n>>3)); + n&=7; + if (n) + y*= (*(const TRealX*)(TwoToNover8+3*n-3)); + return KErrNone; + } + else + { + if (n<0) + { + y.SetZero(); + return KErrUnderflow; + } + else + { + y.SetInfinite(0); + return KErrOverflow; + } + } + } + + + + +EXPORT_C TInt Math::Pow(TReal &aTrg,const TReal &aSrc,const TReal &aPower) +/** +Calculates the value of x raised to the power of y. + +The behaviour conforms to that specified for pow() in the +ISO C Standard ISO/IEC 9899 (Annex F), although floating-point exceptions +are not supported. + +@param aTrg A reference containing the result. +@param aSrc The x argument of the function. +@param aPower The y argument of the function. + +@return KErrNone if successful; + KErrOverflow if the result is +/- infinity; + KErrUnderflow if the result is too small to be represented; + KErrArgument if the result is not a number (NaN). +*/ +// +// Evaluates aSrc raised to the power aPower and places the result in aTrg. +// For non-special values algorithm is aTrg=2^(aPower*log2(aSrc)) +// + { + TRealX x,p; + + TInt ret2=p.Set(aPower); + // pow(x, +/-0) -> 1 for any x, even a NaN + if (p.IsZero()) + { + aTrg=1.0; + return KErrNone; + } + + TInt ret1=x.Set(aSrc); + if (ret1==KErrArgument || ret2==KErrArgument) + { + // pow(+1, y) -> 1 for any y, even a NaN + // XXX First test should not be necessary, but on WINS + // aSrc == 1.0 is true when aSrc is NaN. + if (ret1 != KErrArgument && aSrc == 1.0) + { + aTrg=aSrc; + return KErrNone; + } + SetNaN(aTrg); + return KErrArgument; + } + + // Infinite power + if (ret2==KErrOverflow) + { + // figure out which of these cases we have: + // + // pow(x, -INF) -> +INF for |x| < 1 } flag = 0 + // pow(x, +INF) -> +INF for |x| > 1 } + // pow(x, -INF) -> +0 for |x| > 1 } flag = 1 + // pow(x, +INF) -> +0 for |x| < 1 } + // + // flag = 2 => |x| == 1.0 + // + TInt flag=2; + if (Abs(aSrc)>1.0) + flag=p.iSign&1; + if (Abs(aSrc)<1.0) + flag=1-(p.iSign&1); + if (flag==0) + { + SetInfinite(aTrg,0); + return KErrOverflow; + } + if (flag==1) + { + SetZero(aTrg,0); + return KErrNone; + } + if (Abs(aSrc)==1.0) + { + // pow(-1, +/-INF) -> 1 + aTrg=1.0; + return KErrNone; + } + // This should never happen (i.e. aSrc is NaN, which + // should be taken care of above) + SetNaN(aTrg); + return KErrArgument; + } + + // Negative Base raised to a power + TInt odd=1; + if (x.iSign & 1) + { + TReal pint; + Math::Int(pint,aPower); + if (aPower-pint) // Checks that if aSrc is less than zero, then aPower is integral + { + // pow(-INF, y) -> +0 for y < 0 and not an odd integer + // pow(-INF, y) -> +INF for y > 0 and not an odd integer + // Since we're here, aPower is not integral, so can't be odd, either + if (ret1 == KErrOverflow) + { + if (aPower < 0) + { + SetZero(aTrg); + return KErrNone; + } + else + { + SetInfinite(aTrg,0); + return KErrOverflow; + } + } + SetNaN(aTrg); + return KErrArgument; + } + TReal powerby2=aPower*0.5; + Math::Int(pint,powerby2); + if (powerby2-pint) + odd=(-1); + x.iSign=0; + } + + // Zero or infinity raised to a power + if (x.IsZero() || ret1==KErrOverflow) + { + if (x.IsZero() && p.IsZero()) + { + aTrg=1.0; + return KErrNone; + } + TInt sign=(odd==-1 ? 1 : 0); + if ((x.IsZero() && (p.iSign&1)==0) || (ret1==KErrOverflow && (p.iSign&1))) + { + SetZero(aTrg,sign); + return KErrNone; + } + else + { + SetInfinite(aTrg,sign); + return KErrOverflow; + } + } + + TRealX y; + Log2(y,x); + x=y*p; // this cannot overflow or underflow + TInt r=TwoTox(y,x); + if (odd<0) + y.iSign=1; + TInt r2=y.GetTReal(aTrg); + return (r==KErrNone)?r2:r; + } + +#else // __USE_VFP_MATH + +// definitions come from RVCT math library +extern "C" TReal pow(TReal,TReal); + +EXPORT_C TInt Math::Pow(TReal &aTrg,const TReal &aSrc,const TReal &aPower) + { + aTrg = pow(aSrc,aPower); + if (Math::IsZero(aTrg) && !Math::IsZero(aSrc) && !Math::IsInfinite(aSrc) && !Math::IsInfinite(aPower)) + return KErrUnderflow; + if (Math::IsFinite(aTrg)) + return KErrNone; + if (Math::IsZero(aPower)) // pow(x, +/-0) -> 1 for any x, even a NaN + { + aTrg = 1.0; + return KErrNone; + } + if (Math::IsInfinite(aTrg)) + return KErrOverflow; + if (aSrc==1.0) // pow(+1, y) -> 1 for any y, even a NaN + { + aTrg=aSrc; + return KErrNone; + } + if (Math::IsInfinite(aPower)) + { + if (aSrc == -1.0) // pow(-1, +/-INF) -> 1 + { + aTrg = 1.0; + return KErrNone; + } + if (((Abs(aSrc) < 1) && (aPower < 0)) || // pow(x, -INF) -> +INF for |x| < 1 + ((Abs(aSrc) > 1) && (aPower > 0))) // pow(x, +INF) -> +INF for |x| > 1 + { + SetInfinite(aTrg,0); + return KErrOverflow; + } + } + // pow(-INF, y) -> +INF for y > 0 and not an odd integer + if (Math::IsInfinite(aSrc) && (aSrc < 0) && (aPower > 0)) + { + TBool odd = EFalse; + TReal pint; + Math::Int(pint, aPower); + if (aPower == pint) + { + TReal halfPower = aPower * 0.5; + Math::Int(pint, halfPower); + if (halfPower != pint) + odd = ETrue; + } + if (odd == EFalse) + { + SetInfinite(aTrg,0); + return KErrOverflow; + } + } + + // Otherwise... + SetNaN(aTrg); + return KErrArgument; + } + +#endif