187d2adbcSAndrew Rist /**************************************************************
287d2adbcSAndrew Rist *
387d2adbcSAndrew Rist * Licensed to the Apache Software Foundation (ASF) under one
487d2adbcSAndrew Rist * or more contributor license agreements. See the NOTICE file
587d2adbcSAndrew Rist * distributed with this work for additional information
687d2adbcSAndrew Rist * regarding copyright ownership. The ASF licenses this file
787d2adbcSAndrew Rist * to you under the Apache License, Version 2.0 (the
887d2adbcSAndrew Rist * "License"); you may not use this file except in compliance
987d2adbcSAndrew Rist * with the License. You may obtain a copy of the License at
1087d2adbcSAndrew Rist *
1187d2adbcSAndrew Rist * http://www.apache.org/licenses/LICENSE-2.0
1287d2adbcSAndrew Rist *
1387d2adbcSAndrew Rist * Unless required by applicable law or agreed to in writing,
1487d2adbcSAndrew Rist * software distributed under the License is distributed on an
1587d2adbcSAndrew Rist * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
1687d2adbcSAndrew Rist * KIND, either express or implied. See the License for the
1787d2adbcSAndrew Rist * specific language governing permissions and limitations
1887d2adbcSAndrew Rist * under the License.
1987d2adbcSAndrew Rist *
2087d2adbcSAndrew Rist *************************************************************/
2187d2adbcSAndrew Rist
2287d2adbcSAndrew Rist
23cdf0e10cSrcweir
24cdf0e10cSrcweir // MARKER(update_precomp.py): autogen include statement, do not remove
25cdf0e10cSrcweir #include "precompiled_sal.hxx"
26cdf0e10cSrcweir
27cdf0e10cSrcweir #include "rtl/math.h"
28cdf0e10cSrcweir
29cdf0e10cSrcweir #include "osl/diagnose.h"
30cdf0e10cSrcweir #include "rtl/alloc.h"
31cdf0e10cSrcweir #include "rtl/math.hxx"
32cdf0e10cSrcweir #include "rtl/strbuf.h"
33cdf0e10cSrcweir #include "rtl/string.h"
34cdf0e10cSrcweir #include "rtl/ustrbuf.h"
35cdf0e10cSrcweir #include "rtl/ustring.h"
36cdf0e10cSrcweir #include "sal/mathconf.h"
37cdf0e10cSrcweir #include "sal/types.h"
38cdf0e10cSrcweir
39cdf0e10cSrcweir #include <algorithm>
40cdf0e10cSrcweir #include <float.h>
41cdf0e10cSrcweir #include <limits.h>
42cdf0e10cSrcweir #include <math.h>
43cdf0e10cSrcweir #include <stdlib.h>
44cdf0e10cSrcweir
45cdf0e10cSrcweir
46cdf0e10cSrcweir static int const n10Count = 16;
47cdf0e10cSrcweir static double const n10s[2][n10Count] = {
48cdf0e10cSrcweir { 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8,
49cdf0e10cSrcweir 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16 },
50cdf0e10cSrcweir { 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8,
51cdf0e10cSrcweir 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16 }
52cdf0e10cSrcweir };
53cdf0e10cSrcweir
54cdf0e10cSrcweir // return pow(10.0,nExp) optimized for exponents in the interval [-16,16]
getN10Exp(int nExp)55cdf0e10cSrcweir static double getN10Exp( int nExp )
56cdf0e10cSrcweir {
57cdf0e10cSrcweir if ( nExp < 0 )
58cdf0e10cSrcweir {
59cdf0e10cSrcweir if ( -nExp <= n10Count )
60cdf0e10cSrcweir return n10s[1][-nExp-1];
61cdf0e10cSrcweir else
62cdf0e10cSrcweir return pow( 10.0, static_cast<double>( nExp ) );
63cdf0e10cSrcweir }
64cdf0e10cSrcweir else if ( nExp > 0 )
65cdf0e10cSrcweir {
66cdf0e10cSrcweir if ( nExp <= n10Count )
67cdf0e10cSrcweir return n10s[0][nExp-1];
68cdf0e10cSrcweir else
69cdf0e10cSrcweir return pow( 10.0, static_cast<double>( nExp ) );
70cdf0e10cSrcweir }
71cdf0e10cSrcweir else // ( nExp == 0 )
72cdf0e10cSrcweir return 1.0;
73cdf0e10cSrcweir }
74cdf0e10cSrcweir
75cdf0e10cSrcweir /** Approximation algorithm for erf for 0 < x < 0.65. */
lcl_Erf0065(double x,double & fVal)76cdf0e10cSrcweir void lcl_Erf0065( double x, double& fVal )
77cdf0e10cSrcweir {
78cdf0e10cSrcweir static const double pn[] = {
79cdf0e10cSrcweir 1.12837916709551256,
80cdf0e10cSrcweir 1.35894887627277916E-1,
81cdf0e10cSrcweir 4.03259488531795274E-2,
82cdf0e10cSrcweir 1.20339380863079457E-3,
83cdf0e10cSrcweir 6.49254556481904354E-5
84cdf0e10cSrcweir };
85cdf0e10cSrcweir static const double qn[] = {
86cdf0e10cSrcweir 1.00000000000000000,
87cdf0e10cSrcweir 4.53767041780002545E-1,
88cdf0e10cSrcweir 8.69936222615385890E-2,
89cdf0e10cSrcweir 8.49717371168693357E-3,
90cdf0e10cSrcweir 3.64915280629351082E-4
91cdf0e10cSrcweir };
92cdf0e10cSrcweir double fPSum = 0.0;
93cdf0e10cSrcweir double fQSum = 0.0;
94cdf0e10cSrcweir double fXPow = 1.0;
95cdf0e10cSrcweir for ( unsigned int i = 0; i <= 4; ++i )
96cdf0e10cSrcweir {
97cdf0e10cSrcweir fPSum += pn[i]*fXPow;
98cdf0e10cSrcweir fQSum += qn[i]*fXPow;
99cdf0e10cSrcweir fXPow *= x*x;
100cdf0e10cSrcweir }
101cdf0e10cSrcweir fVal = x * fPSum / fQSum;
102cdf0e10cSrcweir }
103cdf0e10cSrcweir
104cdf0e10cSrcweir /** Approximation algorithm for erfc for 0.65 < x < 6.0. */
lcl_Erfc0600(double x,double & fVal)105cdf0e10cSrcweir void lcl_Erfc0600( double x, double& fVal )
106cdf0e10cSrcweir {
107cdf0e10cSrcweir double fPSum = 0.0;
108cdf0e10cSrcweir double fQSum = 0.0;
109cdf0e10cSrcweir double fXPow = 1.0;
110cdf0e10cSrcweir const double *pn;
111cdf0e10cSrcweir const double *qn;
112cdf0e10cSrcweir
113cdf0e10cSrcweir if ( x < 2.2 )
114cdf0e10cSrcweir {
115cdf0e10cSrcweir static const double pn22[] = {
116cdf0e10cSrcweir 9.99999992049799098E-1,
117cdf0e10cSrcweir 1.33154163936765307,
118cdf0e10cSrcweir 8.78115804155881782E-1,
119cdf0e10cSrcweir 3.31899559578213215E-1,
120cdf0e10cSrcweir 7.14193832506776067E-2,
121cdf0e10cSrcweir 7.06940843763253131E-3
122cdf0e10cSrcweir };
123cdf0e10cSrcweir static const double qn22[] = {
124cdf0e10cSrcweir 1.00000000000000000,
125cdf0e10cSrcweir 2.45992070144245533,
126cdf0e10cSrcweir 2.65383972869775752,
127cdf0e10cSrcweir 1.61876655543871376,
128cdf0e10cSrcweir 5.94651311286481502E-1,
129cdf0e10cSrcweir 1.26579413030177940E-1,
130cdf0e10cSrcweir 1.25304936549413393E-2
131cdf0e10cSrcweir };
132cdf0e10cSrcweir pn = pn22;
133cdf0e10cSrcweir qn = qn22;
134cdf0e10cSrcweir }
135cdf0e10cSrcweir else /* if ( x < 6.0 ) this is true, but the compiler does not know */
136cdf0e10cSrcweir {
137cdf0e10cSrcweir static const double pn60[] = {
138cdf0e10cSrcweir 9.99921140009714409E-1,
139cdf0e10cSrcweir 1.62356584489366647,
140cdf0e10cSrcweir 1.26739901455873222,
141cdf0e10cSrcweir 5.81528574177741135E-1,
142cdf0e10cSrcweir 1.57289620742838702E-1,
143cdf0e10cSrcweir 2.25716982919217555E-2
144cdf0e10cSrcweir };
145cdf0e10cSrcweir static const double qn60[] = {
146cdf0e10cSrcweir 1.00000000000000000,
147cdf0e10cSrcweir 2.75143870676376208,
148cdf0e10cSrcweir 3.37367334657284535,
149cdf0e10cSrcweir 2.38574194785344389,
150cdf0e10cSrcweir 1.05074004614827206,
151cdf0e10cSrcweir 2.78788439273628983E-1,
152cdf0e10cSrcweir 4.00072964526861362E-2
153cdf0e10cSrcweir };
154cdf0e10cSrcweir pn = pn60;
155cdf0e10cSrcweir qn = qn60;
156cdf0e10cSrcweir }
157cdf0e10cSrcweir
158cdf0e10cSrcweir for ( unsigned int i = 0; i < 6; ++i )
159cdf0e10cSrcweir {
160cdf0e10cSrcweir fPSum += pn[i]*fXPow;
161cdf0e10cSrcweir fQSum += qn[i]*fXPow;
162cdf0e10cSrcweir fXPow *= x;
163cdf0e10cSrcweir }
164cdf0e10cSrcweir fQSum += qn[6]*fXPow;
165cdf0e10cSrcweir fVal = exp( -1.0*x*x )* fPSum / fQSum;
166cdf0e10cSrcweir }
167cdf0e10cSrcweir
168cdf0e10cSrcweir /** Approximation algorithm for erfc for 6.0 < x < 26.54 (but used for all
169cdf0e10cSrcweir x > 6.0). */
lcl_Erfc2654(double x,double & fVal)170cdf0e10cSrcweir void lcl_Erfc2654( double x, double& fVal )
171cdf0e10cSrcweir {
172cdf0e10cSrcweir static const double pn[] = {
173cdf0e10cSrcweir 5.64189583547756078E-1,
174cdf0e10cSrcweir 8.80253746105525775,
175cdf0e10cSrcweir 3.84683103716117320E1,
176cdf0e10cSrcweir 4.77209965874436377E1,
177cdf0e10cSrcweir 8.08040729052301677
178cdf0e10cSrcweir };
179cdf0e10cSrcweir static const double qn[] = {
180cdf0e10cSrcweir 1.00000000000000000,
181cdf0e10cSrcweir 1.61020914205869003E1,
182cdf0e10cSrcweir 7.54843505665954743E1,
183cdf0e10cSrcweir 1.12123870801026015E2,
184cdf0e10cSrcweir 3.73997570145040850E1
185cdf0e10cSrcweir };
186cdf0e10cSrcweir
187cdf0e10cSrcweir double fPSum = 0.0;
188cdf0e10cSrcweir double fQSum = 0.0;
189cdf0e10cSrcweir double fXPow = 1.0;
190cdf0e10cSrcweir
191cdf0e10cSrcweir for ( unsigned int i = 0; i <= 4; ++i )
192cdf0e10cSrcweir {
193cdf0e10cSrcweir fPSum += pn[i]*fXPow;
194cdf0e10cSrcweir fQSum += qn[i]*fXPow;
195cdf0e10cSrcweir fXPow /= x*x;
196cdf0e10cSrcweir }
197cdf0e10cSrcweir fVal = exp(-1.0*x*x)*fPSum / (x*fQSum);
198cdf0e10cSrcweir }
199cdf0e10cSrcweir
200cdf0e10cSrcweir namespace {
201cdf0e10cSrcweir
202cdf0e10cSrcweir double const nKorrVal[] = {
203cdf0e10cSrcweir 0, 9e-1, 9e-2, 9e-3, 9e-4, 9e-5, 9e-6, 9e-7, 9e-8,
204cdf0e10cSrcweir 9e-9, 9e-10, 9e-11, 9e-12, 9e-13, 9e-14, 9e-15
205cdf0e10cSrcweir };
206cdf0e10cSrcweir
207cdf0e10cSrcweir struct StringTraits
208cdf0e10cSrcweir {
209cdf0e10cSrcweir typedef sal_Char Char;
210cdf0e10cSrcweir
211cdf0e10cSrcweir typedef rtl_String String;
212cdf0e10cSrcweir
createString__anon0f9323920111::StringTraits213cdf0e10cSrcweir static inline void createString(rtl_String ** pString,
214cdf0e10cSrcweir sal_Char const * pChars, sal_Int32 nLen)
215cdf0e10cSrcweir {
216cdf0e10cSrcweir rtl_string_newFromStr_WithLength(pString, pChars, nLen);
217cdf0e10cSrcweir }
218cdf0e10cSrcweir
createBuffer__anon0f9323920111::StringTraits219cdf0e10cSrcweir static inline void createBuffer(rtl_String ** pBuffer,
220cdf0e10cSrcweir sal_Int32 * pCapacity)
221cdf0e10cSrcweir {
222cdf0e10cSrcweir rtl_string_new_WithLength(pBuffer, *pCapacity);
223cdf0e10cSrcweir }
224cdf0e10cSrcweir
appendChar__anon0f9323920111::StringTraits225cdf0e10cSrcweir static inline void appendChar(rtl_String ** pBuffer, sal_Int32 * pCapacity,
226cdf0e10cSrcweir sal_Int32 * pOffset, sal_Char cChar)
227cdf0e10cSrcweir {
228cdf0e10cSrcweir rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, &cChar, 1);
229cdf0e10cSrcweir ++*pOffset;
230cdf0e10cSrcweir }
231cdf0e10cSrcweir
appendChars__anon0f9323920111::StringTraits232cdf0e10cSrcweir static inline void appendChars(rtl_String ** pBuffer, sal_Int32 * pCapacity,
233cdf0e10cSrcweir sal_Int32 * pOffset, sal_Char const * pChars,
234cdf0e10cSrcweir sal_Int32 nLen)
235cdf0e10cSrcweir {
236cdf0e10cSrcweir rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, pChars, nLen);
237cdf0e10cSrcweir *pOffset += nLen;
238cdf0e10cSrcweir }
239cdf0e10cSrcweir
appendAscii__anon0f9323920111::StringTraits240cdf0e10cSrcweir static inline void appendAscii(rtl_String ** pBuffer, sal_Int32 * pCapacity,
241cdf0e10cSrcweir sal_Int32 * pOffset, sal_Char const * pStr,
242cdf0e10cSrcweir sal_Int32 nLen)
243cdf0e10cSrcweir {
244cdf0e10cSrcweir rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, pStr, nLen);
245cdf0e10cSrcweir *pOffset += nLen;
246cdf0e10cSrcweir }
247cdf0e10cSrcweir };
248cdf0e10cSrcweir
249cdf0e10cSrcweir struct UStringTraits
250cdf0e10cSrcweir {
251cdf0e10cSrcweir typedef sal_Unicode Char;
252cdf0e10cSrcweir
253cdf0e10cSrcweir typedef rtl_uString String;
254cdf0e10cSrcweir
createString__anon0f9323920111::UStringTraits255cdf0e10cSrcweir static inline void createString(rtl_uString ** pString,
256cdf0e10cSrcweir sal_Unicode const * pChars, sal_Int32 nLen)
257cdf0e10cSrcweir {
258cdf0e10cSrcweir rtl_uString_newFromStr_WithLength(pString, pChars, nLen);
259cdf0e10cSrcweir }
260cdf0e10cSrcweir
createBuffer__anon0f9323920111::UStringTraits261cdf0e10cSrcweir static inline void createBuffer(rtl_uString ** pBuffer,
262cdf0e10cSrcweir sal_Int32 * pCapacity)
263cdf0e10cSrcweir {
264cdf0e10cSrcweir rtl_uString_new_WithLength(pBuffer, *pCapacity);
265cdf0e10cSrcweir }
266cdf0e10cSrcweir
appendChar__anon0f9323920111::UStringTraits267cdf0e10cSrcweir static inline void appendChar(rtl_uString ** pBuffer, sal_Int32 * pCapacity,
268cdf0e10cSrcweir sal_Int32 * pOffset, sal_Unicode cChar)
269cdf0e10cSrcweir {
270cdf0e10cSrcweir rtl_uStringbuffer_insert(pBuffer, pCapacity, *pOffset, &cChar, 1);
271cdf0e10cSrcweir ++*pOffset;
272cdf0e10cSrcweir }
273cdf0e10cSrcweir
appendChars__anon0f9323920111::UStringTraits274cdf0e10cSrcweir static inline void appendChars(rtl_uString ** pBuffer,
275cdf0e10cSrcweir sal_Int32 * pCapacity, sal_Int32 * pOffset,
276cdf0e10cSrcweir sal_Unicode const * pChars, sal_Int32 nLen)
277cdf0e10cSrcweir {
278cdf0e10cSrcweir rtl_uStringbuffer_insert(pBuffer, pCapacity, *pOffset, pChars, nLen);
279cdf0e10cSrcweir *pOffset += nLen;
280cdf0e10cSrcweir }
281cdf0e10cSrcweir
appendAscii__anon0f9323920111::UStringTraits282cdf0e10cSrcweir static inline void appendAscii(rtl_uString ** pBuffer,
283cdf0e10cSrcweir sal_Int32 * pCapacity, sal_Int32 * pOffset,
284cdf0e10cSrcweir sal_Char const * pStr, sal_Int32 nLen)
285cdf0e10cSrcweir {
286cdf0e10cSrcweir rtl_uStringbuffer_insert_ascii(pBuffer, pCapacity, *pOffset, pStr,
287cdf0e10cSrcweir nLen);
288cdf0e10cSrcweir *pOffset += nLen;
289cdf0e10cSrcweir }
290cdf0e10cSrcweir };
291cdf0e10cSrcweir
292cdf0e10cSrcweir
293cdf0e10cSrcweir // Solaris C++ 5.2 compiler has problems when "StringT ** pResult" is
294cdf0e10cSrcweir // "typename T::String ** pResult" instead:
295cdf0e10cSrcweir template< typename T, typename StringT >
doubleToString(StringT ** pResult,sal_Int32 * pResultCapacity,sal_Int32 nResultOffset,double fValue,rtl_math_StringFormat eFormat,sal_Int32 nDecPlaces,typename T::Char cDecSeparator,sal_Int32 const * pGroups,typename T::Char cGroupSeparator,bool bEraseTrailingDecZeros)296cdf0e10cSrcweir inline void doubleToString(StringT ** pResult,
297cdf0e10cSrcweir sal_Int32 * pResultCapacity, sal_Int32 nResultOffset,
298cdf0e10cSrcweir double fValue, rtl_math_StringFormat eFormat,
299cdf0e10cSrcweir sal_Int32 nDecPlaces, typename T::Char cDecSeparator,
300cdf0e10cSrcweir sal_Int32 const * pGroups,
301cdf0e10cSrcweir typename T::Char cGroupSeparator,
302cdf0e10cSrcweir bool bEraseTrailingDecZeros)
303cdf0e10cSrcweir {
304cdf0e10cSrcweir static double const nRoundVal[] = {
305cdf0e10cSrcweir 5.0e+0, 0.5e+0, 0.5e-1, 0.5e-2, 0.5e-3, 0.5e-4, 0.5e-5, 0.5e-6,
306cdf0e10cSrcweir 0.5e-7, 0.5e-8, 0.5e-9, 0.5e-10,0.5e-11,0.5e-12,0.5e-13,0.5e-14
307cdf0e10cSrcweir };
308cdf0e10cSrcweir
309cdf0e10cSrcweir // sign adjustment, instead of testing for fValue<0.0 this will also fetch
310cdf0e10cSrcweir // -0.0
311cdf0e10cSrcweir bool bSign = rtl::math::isSignBitSet( fValue );
312cdf0e10cSrcweir if( bSign )
313cdf0e10cSrcweir fValue = -fValue;
314cdf0e10cSrcweir
315cdf0e10cSrcweir if ( rtl::math::isNan( fValue ) )
316cdf0e10cSrcweir {
317cdf0e10cSrcweir // #i112652# XMLSchema-2
318cdf0e10cSrcweir sal_Int32 nCapacity = RTL_CONSTASCII_LENGTH("NaN");
319cdf0e10cSrcweir if (pResultCapacity == 0)
320cdf0e10cSrcweir {
321cdf0e10cSrcweir pResultCapacity = &nCapacity;
322cdf0e10cSrcweir T::createBuffer(pResult, pResultCapacity);
323cdf0e10cSrcweir nResultOffset = 0;
324cdf0e10cSrcweir }
325cdf0e10cSrcweir T::appendAscii(pResult, pResultCapacity, &nResultOffset,
326cdf0e10cSrcweir RTL_CONSTASCII_STRINGPARAM("NaN"));
327cdf0e10cSrcweir
328cdf0e10cSrcweir return;
329cdf0e10cSrcweir }
330cdf0e10cSrcweir
331cdf0e10cSrcweir bool bHuge = fValue == HUGE_VAL; // g++ 3.0.1 requires it this way...
332cdf0e10cSrcweir if ( bHuge || rtl::math::isInf( fValue ) )
333cdf0e10cSrcweir {
334cdf0e10cSrcweir // #i112652# XMLSchema-2
335cdf0e10cSrcweir sal_Int32 nCapacity = RTL_CONSTASCII_LENGTH("-INF");
336cdf0e10cSrcweir if (pResultCapacity == 0)
337cdf0e10cSrcweir {
338cdf0e10cSrcweir pResultCapacity = &nCapacity;
339cdf0e10cSrcweir T::createBuffer(pResult, pResultCapacity);
340cdf0e10cSrcweir nResultOffset = 0;
341cdf0e10cSrcweir }
342cdf0e10cSrcweir if ( bSign )
343cdf0e10cSrcweir T::appendAscii(pResult, pResultCapacity, &nResultOffset,
344cdf0e10cSrcweir RTL_CONSTASCII_STRINGPARAM("-"));
345cdf0e10cSrcweir T::appendAscii(pResult, pResultCapacity, &nResultOffset,
346cdf0e10cSrcweir RTL_CONSTASCII_STRINGPARAM("INF"));
347cdf0e10cSrcweir
348cdf0e10cSrcweir return;
349cdf0e10cSrcweir }
350cdf0e10cSrcweir
351cdf0e10cSrcweir // find the exponent
352cdf0e10cSrcweir int nExp = 0;
353cdf0e10cSrcweir if ( fValue > 0.0 )
354cdf0e10cSrcweir {
355cdf0e10cSrcweir nExp = static_cast< int >( floor( log10( fValue ) ) );
356cdf0e10cSrcweir fValue /= getN10Exp( nExp );
357cdf0e10cSrcweir }
358cdf0e10cSrcweir
359cdf0e10cSrcweir switch ( eFormat )
360cdf0e10cSrcweir {
361cdf0e10cSrcweir case rtl_math_StringFormat_Automatic :
362cdf0e10cSrcweir { // E or F depending on exponent magnitude
363cdf0e10cSrcweir int nPrec;
364cdf0e10cSrcweir if ( nExp <= -15 || nExp >= 15 ) // #58531# was <-16, >16
365cdf0e10cSrcweir {
366cdf0e10cSrcweir nPrec = 14;
367cdf0e10cSrcweir eFormat = rtl_math_StringFormat_E;
368cdf0e10cSrcweir }
369cdf0e10cSrcweir else
370cdf0e10cSrcweir {
371cdf0e10cSrcweir if ( nExp < 14 )
372cdf0e10cSrcweir {
373cdf0e10cSrcweir nPrec = 15 - nExp - 1;
374cdf0e10cSrcweir eFormat = rtl_math_StringFormat_F;
375cdf0e10cSrcweir }
376cdf0e10cSrcweir else
377cdf0e10cSrcweir {
378cdf0e10cSrcweir nPrec = 15;
379cdf0e10cSrcweir eFormat = rtl_math_StringFormat_F;
380cdf0e10cSrcweir }
381cdf0e10cSrcweir }
382cdf0e10cSrcweir if ( nDecPlaces == rtl_math_DecimalPlaces_Max )
383cdf0e10cSrcweir nDecPlaces = nPrec;
384cdf0e10cSrcweir }
385cdf0e10cSrcweir break;
386cdf0e10cSrcweir case rtl_math_StringFormat_G :
387cdf0e10cSrcweir { // G-Point, similar to sprintf %G
388cdf0e10cSrcweir if ( nDecPlaces == rtl_math_DecimalPlaces_DefaultSignificance )
389cdf0e10cSrcweir nDecPlaces = 6;
390cdf0e10cSrcweir if ( nExp < -4 || nExp >= nDecPlaces )
391cdf0e10cSrcweir {
392cdf0e10cSrcweir nDecPlaces = std::max< sal_Int32 >( 1, nDecPlaces - 1 );
393cdf0e10cSrcweir eFormat = rtl_math_StringFormat_E;
394cdf0e10cSrcweir }
395cdf0e10cSrcweir else
396cdf0e10cSrcweir {
397cdf0e10cSrcweir nDecPlaces = std::max< sal_Int32 >( 0, nDecPlaces - nExp - 1 );
398cdf0e10cSrcweir eFormat = rtl_math_StringFormat_F;
399cdf0e10cSrcweir }
400cdf0e10cSrcweir }
401cdf0e10cSrcweir break;
402cdf0e10cSrcweir default:
403cdf0e10cSrcweir break;
404cdf0e10cSrcweir }
405cdf0e10cSrcweir
406cdf0e10cSrcweir sal_Int32 nDigits = nDecPlaces + 1;
407cdf0e10cSrcweir
408cdf0e10cSrcweir if( eFormat == rtl_math_StringFormat_F )
409cdf0e10cSrcweir nDigits += nExp;
410cdf0e10cSrcweir
411cdf0e10cSrcweir // Round the number
412cdf0e10cSrcweir if( nDigits >= 0 )
413cdf0e10cSrcweir {
414cdf0e10cSrcweir if( ( fValue += nRoundVal[ nDigits > 15 ? 15 : nDigits ] ) >= 10 )
415cdf0e10cSrcweir {
416cdf0e10cSrcweir fValue = 1.0;
417cdf0e10cSrcweir nExp++;
418cdf0e10cSrcweir if( eFormat == rtl_math_StringFormat_F )
419cdf0e10cSrcweir nDigits++;
420cdf0e10cSrcweir }
421cdf0e10cSrcweir }
422cdf0e10cSrcweir
423cdf0e10cSrcweir static sal_Int32 const nBufMax = 256;
424cdf0e10cSrcweir typename T::Char aBuf[nBufMax];
425cdf0e10cSrcweir typename T::Char * pBuf;
426cdf0e10cSrcweir sal_Int32 nBuf = static_cast< sal_Int32 >
427cdf0e10cSrcweir ( nDigits <= 0 ? std::max< sal_Int32 >( nDecPlaces, abs(nExp) )
428cdf0e10cSrcweir : nDigits + nDecPlaces ) + 10 + (pGroups ? abs(nDigits) * 2 : 0);
429cdf0e10cSrcweir if ( nBuf > nBufMax )
430cdf0e10cSrcweir {
431cdf0e10cSrcweir pBuf = reinterpret_cast< typename T::Char * >(
432cdf0e10cSrcweir rtl_allocateMemory(nBuf * sizeof (typename T::Char)));
433cdf0e10cSrcweir OSL_ENSURE(pBuf != 0, "Out of memory");
434cdf0e10cSrcweir }
435cdf0e10cSrcweir else
436cdf0e10cSrcweir pBuf = aBuf;
437cdf0e10cSrcweir typename T::Char * p = pBuf;
438cdf0e10cSrcweir if ( bSign )
439cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('-');
440cdf0e10cSrcweir
441cdf0e10cSrcweir bool bHasDec = false;
442cdf0e10cSrcweir
443cdf0e10cSrcweir int nDecPos;
444cdf0e10cSrcweir // Check for F format and number < 1
445cdf0e10cSrcweir if( eFormat == rtl_math_StringFormat_F )
446cdf0e10cSrcweir {
447cdf0e10cSrcweir if( nExp < 0 )
448cdf0e10cSrcweir {
449cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
450cdf0e10cSrcweir if ( nDecPlaces > 0 )
451cdf0e10cSrcweir {
452cdf0e10cSrcweir *p++ = cDecSeparator;
453cdf0e10cSrcweir bHasDec = true;
454cdf0e10cSrcweir }
455cdf0e10cSrcweir sal_Int32 i = ( nDigits <= 0 ? nDecPlaces : -nExp - 1 );
456cdf0e10cSrcweir while( (i--) > 0 )
457cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
458cdf0e10cSrcweir nDecPos = 0;
459cdf0e10cSrcweir }
460cdf0e10cSrcweir else
461cdf0e10cSrcweir nDecPos = nExp + 1;
462cdf0e10cSrcweir }
463cdf0e10cSrcweir else
464cdf0e10cSrcweir nDecPos = 1;
465cdf0e10cSrcweir
466cdf0e10cSrcweir int nGrouping = 0, nGroupSelector = 0, nGroupExceed = 0;
467cdf0e10cSrcweir if ( nDecPos > 1 && pGroups && pGroups[0] && cGroupSeparator )
468cdf0e10cSrcweir {
469cdf0e10cSrcweir while ( nGrouping + pGroups[nGroupSelector] < nDecPos )
470cdf0e10cSrcweir {
471cdf0e10cSrcweir nGrouping += pGroups[ nGroupSelector ];
472cdf0e10cSrcweir if ( pGroups[nGroupSelector+1] )
473cdf0e10cSrcweir {
474cdf0e10cSrcweir if ( nGrouping + pGroups[nGroupSelector+1] >= nDecPos )
475cdf0e10cSrcweir break; // while
476cdf0e10cSrcweir ++nGroupSelector;
477cdf0e10cSrcweir }
478cdf0e10cSrcweir else if ( !nGroupExceed )
479cdf0e10cSrcweir nGroupExceed = nGrouping;
480cdf0e10cSrcweir }
481cdf0e10cSrcweir }
482cdf0e10cSrcweir
483cdf0e10cSrcweir // print the number
484cdf0e10cSrcweir if( nDigits > 0 )
485cdf0e10cSrcweir {
486cdf0e10cSrcweir for ( int i = 0; ; i++ )
487cdf0e10cSrcweir {
488cdf0e10cSrcweir if( i < 15 )
489cdf0e10cSrcweir {
490cdf0e10cSrcweir int nDigit;
491cdf0e10cSrcweir if (nDigits-1 == 0 && i > 0 && i < 14)
492cdf0e10cSrcweir nDigit = static_cast< int >( floor( fValue
493cdf0e10cSrcweir + nKorrVal[15-i] ) );
494cdf0e10cSrcweir else
495cdf0e10cSrcweir nDigit = static_cast< int >( fValue + 1E-15 );
496cdf0e10cSrcweir if (nDigit >= 10)
497cdf0e10cSrcweir { // after-treatment of up-rounding to the next decade
498cdf0e10cSrcweir sal_Int32 sLen = static_cast< long >(p-pBuf)-1;
499cdf0e10cSrcweir if (sLen == -1)
500cdf0e10cSrcweir {
501cdf0e10cSrcweir p = pBuf;
502cdf0e10cSrcweir if ( eFormat == rtl_math_StringFormat_F )
503cdf0e10cSrcweir {
504cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('1');
505cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
506cdf0e10cSrcweir }
507cdf0e10cSrcweir else
508cdf0e10cSrcweir {
509cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('1');
510cdf0e10cSrcweir *p++ = cDecSeparator;
511cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
512cdf0e10cSrcweir nExp++;
513cdf0e10cSrcweir bHasDec = true;
514cdf0e10cSrcweir }
515cdf0e10cSrcweir }
516cdf0e10cSrcweir else
517cdf0e10cSrcweir {
518cdf0e10cSrcweir for (sal_Int32 j = sLen; j >= 0; j--)
519cdf0e10cSrcweir {
520cdf0e10cSrcweir typename T::Char cS = pBuf[j];
521cdf0e10cSrcweir if (cS != cDecSeparator)
522cdf0e10cSrcweir {
523cdf0e10cSrcweir if ( cS != static_cast< typename T::Char >('9'))
524cdf0e10cSrcweir {
525cdf0e10cSrcweir pBuf[j] = ++cS;
526cdf0e10cSrcweir j = -1; // break loop
527cdf0e10cSrcweir }
528cdf0e10cSrcweir else
529cdf0e10cSrcweir {
530cdf0e10cSrcweir pBuf[j]
531cdf0e10cSrcweir = static_cast< typename T::Char >('0');
532cdf0e10cSrcweir if (j == 0)
533cdf0e10cSrcweir {
534cdf0e10cSrcweir if ( eFormat == rtl_math_StringFormat_F)
535cdf0e10cSrcweir { // insert '1'
536cdf0e10cSrcweir typename T::Char * px = p++;
537cdf0e10cSrcweir while ( pBuf < px )
538cdf0e10cSrcweir {
539cdf0e10cSrcweir *px = *(px-1);
540cdf0e10cSrcweir px--;
541cdf0e10cSrcweir }
542cdf0e10cSrcweir pBuf[0] = static_cast<
543cdf0e10cSrcweir typename T::Char >('1');
544cdf0e10cSrcweir }
545cdf0e10cSrcweir else
546cdf0e10cSrcweir {
547cdf0e10cSrcweir pBuf[j] = static_cast<
548cdf0e10cSrcweir typename T::Char >('1');
549cdf0e10cSrcweir nExp++;
550cdf0e10cSrcweir }
551cdf0e10cSrcweir }
552cdf0e10cSrcweir }
553cdf0e10cSrcweir }
554cdf0e10cSrcweir }
555cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
556cdf0e10cSrcweir }
557cdf0e10cSrcweir fValue = 0.0;
558cdf0e10cSrcweir }
559cdf0e10cSrcweir else
560cdf0e10cSrcweir {
561cdf0e10cSrcweir *p++ = static_cast< typename T::Char >(
562cdf0e10cSrcweir nDigit + static_cast< typename T::Char >('0') );
563cdf0e10cSrcweir fValue = ( fValue - nDigit ) * 10.0;
564cdf0e10cSrcweir }
565cdf0e10cSrcweir }
566cdf0e10cSrcweir else
567cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
568cdf0e10cSrcweir if( !--nDigits )
569cdf0e10cSrcweir break; // for
570cdf0e10cSrcweir if( nDecPos )
571cdf0e10cSrcweir {
572cdf0e10cSrcweir if( !--nDecPos )
573cdf0e10cSrcweir {
574cdf0e10cSrcweir *p++ = cDecSeparator;
575cdf0e10cSrcweir bHasDec = true;
576cdf0e10cSrcweir }
577cdf0e10cSrcweir else if ( nDecPos == nGrouping )
578cdf0e10cSrcweir {
579cdf0e10cSrcweir *p++ = cGroupSeparator;
580cdf0e10cSrcweir nGrouping -= pGroups[ nGroupSelector ];
581cdf0e10cSrcweir if ( nGroupSelector && nGrouping < nGroupExceed )
582cdf0e10cSrcweir --nGroupSelector;
583cdf0e10cSrcweir }
584cdf0e10cSrcweir }
585cdf0e10cSrcweir }
586cdf0e10cSrcweir }
587cdf0e10cSrcweir
588cdf0e10cSrcweir if ( !bHasDec && eFormat == rtl_math_StringFormat_F )
589cdf0e10cSrcweir { // nDecPlaces < 0 did round the value
590cdf0e10cSrcweir while ( --nDecPos > 0 )
591cdf0e10cSrcweir { // fill before decimal point
592cdf0e10cSrcweir if ( nDecPos == nGrouping )
593cdf0e10cSrcweir {
594cdf0e10cSrcweir *p++ = cGroupSeparator;
595cdf0e10cSrcweir nGrouping -= pGroups[ nGroupSelector ];
596cdf0e10cSrcweir if ( nGroupSelector && nGrouping < nGroupExceed )
597cdf0e10cSrcweir --nGroupSelector;
598cdf0e10cSrcweir }
599cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('0');
600cdf0e10cSrcweir }
601cdf0e10cSrcweir }
602cdf0e10cSrcweir
603cdf0e10cSrcweir if ( bEraseTrailingDecZeros && bHasDec && p > pBuf )
604cdf0e10cSrcweir {
605cdf0e10cSrcweir while ( *(p-1) == static_cast< typename T::Char >('0') )
606cdf0e10cSrcweir p--;
607cdf0e10cSrcweir if ( *(p-1) == cDecSeparator )
608cdf0e10cSrcweir p--;
609cdf0e10cSrcweir }
610cdf0e10cSrcweir
611cdf0e10cSrcweir // Print the exponent ('E', followed by '+' or '-', followed by exactly
612cdf0e10cSrcweir // three digits). The code in rtl_[u]str_valueOf{Float|Double} relies on
613cdf0e10cSrcweir // this format.
614cdf0e10cSrcweir if( eFormat == rtl_math_StringFormat_E )
615cdf0e10cSrcweir {
616cdf0e10cSrcweir if ( p == pBuf )
617cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('1');
618cdf0e10cSrcweir // maybe no nDigits if nDecPlaces < 0
619cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('E');
620cdf0e10cSrcweir if( nExp < 0 )
621cdf0e10cSrcweir {
622cdf0e10cSrcweir nExp = -nExp;
623cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('-');
624cdf0e10cSrcweir }
625cdf0e10cSrcweir else
626cdf0e10cSrcweir *p++ = static_cast< typename T::Char >('+');
627cdf0e10cSrcweir // if (nExp >= 100 )
628cdf0e10cSrcweir *p++ = static_cast< typename T::Char >(
629cdf0e10cSrcweir nExp / 100 + static_cast< typename T::Char >('0') );
630cdf0e10cSrcweir nExp %= 100;
631cdf0e10cSrcweir *p++ = static_cast< typename T::Char >(
632cdf0e10cSrcweir nExp / 10 + static_cast< typename T::Char >('0') );
633cdf0e10cSrcweir *p++ = static_cast< typename T::Char >(
634cdf0e10cSrcweir nExp % 10 + static_cast< typename T::Char >('0') );
635cdf0e10cSrcweir }
636cdf0e10cSrcweir
637cdf0e10cSrcweir if (pResultCapacity == 0)
638cdf0e10cSrcweir T::createString(pResult, pBuf, p - pBuf);
639cdf0e10cSrcweir else
640cdf0e10cSrcweir T::appendChars(pResult, pResultCapacity, &nResultOffset, pBuf,
641cdf0e10cSrcweir p - pBuf);
642cdf0e10cSrcweir
643cdf0e10cSrcweir if ( pBuf != &aBuf[0] )
644cdf0e10cSrcweir rtl_freeMemory(pBuf);
645cdf0e10cSrcweir }
646cdf0e10cSrcweir
647cdf0e10cSrcweir }
648cdf0e10cSrcweir
rtl_math_doubleToString(rtl_String ** pResult,sal_Int32 * pResultCapacity,sal_Int32 nResultOffset,double fValue,rtl_math_StringFormat eFormat,sal_Int32 nDecPlaces,sal_Char cDecSeparator,sal_Int32 const * pGroups,sal_Char cGroupSeparator,sal_Bool bEraseTrailingDecZeros)649cdf0e10cSrcweir void SAL_CALL rtl_math_doubleToString(rtl_String ** pResult,
650cdf0e10cSrcweir sal_Int32 * pResultCapacity,
651cdf0e10cSrcweir sal_Int32 nResultOffset, double fValue,
652cdf0e10cSrcweir rtl_math_StringFormat eFormat,
653cdf0e10cSrcweir sal_Int32 nDecPlaces,
654cdf0e10cSrcweir sal_Char cDecSeparator,
655cdf0e10cSrcweir sal_Int32 const * pGroups,
656cdf0e10cSrcweir sal_Char cGroupSeparator,
657cdf0e10cSrcweir sal_Bool bEraseTrailingDecZeros)
658cdf0e10cSrcweir SAL_THROW_EXTERN_C()
659cdf0e10cSrcweir {
660cdf0e10cSrcweir doubleToString< StringTraits, StringTraits::String >(
661cdf0e10cSrcweir pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
662cdf0e10cSrcweir cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
663cdf0e10cSrcweir }
664cdf0e10cSrcweir
rtl_math_doubleToUString(rtl_uString ** pResult,sal_Int32 * pResultCapacity,sal_Int32 nResultOffset,double fValue,rtl_math_StringFormat eFormat,sal_Int32 nDecPlaces,sal_Unicode cDecSeparator,sal_Int32 const * pGroups,sal_Unicode cGroupSeparator,sal_Bool bEraseTrailingDecZeros)665cdf0e10cSrcweir void SAL_CALL rtl_math_doubleToUString(rtl_uString ** pResult,
666cdf0e10cSrcweir sal_Int32 * pResultCapacity,
667cdf0e10cSrcweir sal_Int32 nResultOffset, double fValue,
668cdf0e10cSrcweir rtl_math_StringFormat eFormat,
669cdf0e10cSrcweir sal_Int32 nDecPlaces,
670cdf0e10cSrcweir sal_Unicode cDecSeparator,
671cdf0e10cSrcweir sal_Int32 const * pGroups,
672cdf0e10cSrcweir sal_Unicode cGroupSeparator,
673cdf0e10cSrcweir sal_Bool bEraseTrailingDecZeros)
674cdf0e10cSrcweir SAL_THROW_EXTERN_C()
675cdf0e10cSrcweir {
676cdf0e10cSrcweir doubleToString< UStringTraits, UStringTraits::String >(
677cdf0e10cSrcweir pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
678cdf0e10cSrcweir cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
679cdf0e10cSrcweir }
680cdf0e10cSrcweir
681cdf0e10cSrcweir
682cdf0e10cSrcweir namespace {
683cdf0e10cSrcweir
684cdf0e10cSrcweir // if nExp * 10 + nAdd would result in overflow
long10Overflow(long & nExp,int nAdd)685cdf0e10cSrcweir inline bool long10Overflow( long& nExp, int nAdd )
686cdf0e10cSrcweir {
687cdf0e10cSrcweir if ( nExp > (LONG_MAX/10)
688cdf0e10cSrcweir || (nExp == (LONG_MAX/10) && nAdd > (LONG_MAX%10)) )
689cdf0e10cSrcweir {
690cdf0e10cSrcweir nExp = LONG_MAX;
691cdf0e10cSrcweir return true;
692cdf0e10cSrcweir }
693cdf0e10cSrcweir return false;
694cdf0e10cSrcweir }
695cdf0e10cSrcweir
696cdf0e10cSrcweir // We are only concerned about ASCII arabic numerical digits here
697cdf0e10cSrcweir template< typename CharT >
isDigit(CharT c)698cdf0e10cSrcweir inline bool isDigit( CharT c )
699cdf0e10cSrcweir {
700cdf0e10cSrcweir return 0x30 <= c && c <= 0x39;
701cdf0e10cSrcweir }
702cdf0e10cSrcweir
703cdf0e10cSrcweir template< typename CharT >
stringToDouble(CharT const * pBegin,CharT const * pEnd,CharT cDecSeparator,CharT cGroupSeparator,rtl_math_ConversionStatus * pStatus,CharT const ** pParsedEnd)704cdf0e10cSrcweir inline double stringToDouble(CharT const * pBegin, CharT const * pEnd,
705cdf0e10cSrcweir CharT cDecSeparator, CharT cGroupSeparator,
706cdf0e10cSrcweir rtl_math_ConversionStatus * pStatus,
707cdf0e10cSrcweir CharT const ** pParsedEnd)
708cdf0e10cSrcweir {
709cdf0e10cSrcweir double fVal = 0.0;
710cdf0e10cSrcweir rtl_math_ConversionStatus eStatus = rtl_math_ConversionStatus_Ok;
711cdf0e10cSrcweir
712cdf0e10cSrcweir CharT const * p0 = pBegin;
713cdf0e10cSrcweir while (p0 != pEnd && (*p0 == CharT(' ') || *p0 == CharT('\t')))
714cdf0e10cSrcweir ++p0;
715cdf0e10cSrcweir bool bSign;
716cdf0e10cSrcweir if (p0 != pEnd && *p0 == CharT('-'))
717cdf0e10cSrcweir {
718cdf0e10cSrcweir bSign = true;
719cdf0e10cSrcweir ++p0;
720cdf0e10cSrcweir }
721cdf0e10cSrcweir else
722cdf0e10cSrcweir {
723cdf0e10cSrcweir bSign = false;
724cdf0e10cSrcweir if (p0 != pEnd && *p0 == CharT('+'))
725cdf0e10cSrcweir ++p0;
726cdf0e10cSrcweir }
727cdf0e10cSrcweir CharT const * p = p0;
728cdf0e10cSrcweir bool bDone = false;
729cdf0e10cSrcweir
730cdf0e10cSrcweir // #i112652# XMLSchema-2
731cdf0e10cSrcweir if (3 >= (pEnd - p))
732cdf0e10cSrcweir {
733cdf0e10cSrcweir if ((CharT('N') == p[0]) && (CharT('a') == p[1])
734cdf0e10cSrcweir && (CharT('N') == p[2]))
735cdf0e10cSrcweir {
736cdf0e10cSrcweir p += 3;
737cdf0e10cSrcweir rtl::math::setNan( &fVal );
738cdf0e10cSrcweir bDone = true;
739cdf0e10cSrcweir }
740cdf0e10cSrcweir else if ((CharT('I') == p[0]) && (CharT('N') == p[1])
741cdf0e10cSrcweir && (CharT('F') == p[2]))
742cdf0e10cSrcweir {
743cdf0e10cSrcweir p += 3;
744cdf0e10cSrcweir fVal = HUGE_VAL;
745cdf0e10cSrcweir eStatus = rtl_math_ConversionStatus_OutOfRange;
746cdf0e10cSrcweir bDone = true;
747cdf0e10cSrcweir }
748cdf0e10cSrcweir }
749cdf0e10cSrcweir
750cdf0e10cSrcweir if (!bDone) // do not recognize e.g. NaN1.23
751cdf0e10cSrcweir {
752cdf0e10cSrcweir // leading zeros and group separators may be safely ignored
753cdf0e10cSrcweir while (p != pEnd && (*p == CharT('0') || *p == cGroupSeparator))
754cdf0e10cSrcweir ++p;
755cdf0e10cSrcweir
756cdf0e10cSrcweir long nValExp = 0; // carry along exponent of mantissa
757cdf0e10cSrcweir
758cdf0e10cSrcweir // integer part of mantissa
759cdf0e10cSrcweir for (; p != pEnd; ++p)
760cdf0e10cSrcweir {
761cdf0e10cSrcweir CharT c = *p;
762cdf0e10cSrcweir if (isDigit(c))
763cdf0e10cSrcweir {
764cdf0e10cSrcweir fVal = fVal * 10.0 + static_cast< double >( c - CharT('0') );
765cdf0e10cSrcweir ++nValExp;
766cdf0e10cSrcweir }
767cdf0e10cSrcweir else if (c != cGroupSeparator)
768cdf0e10cSrcweir break;
769cdf0e10cSrcweir }
770cdf0e10cSrcweir
771cdf0e10cSrcweir // fraction part of mantissa
772cdf0e10cSrcweir if (p != pEnd && *p == cDecSeparator)
773cdf0e10cSrcweir {
774cdf0e10cSrcweir ++p;
775cdf0e10cSrcweir double fFrac = 0.0;
776cdf0e10cSrcweir long nFracExp = 0;
777cdf0e10cSrcweir while (p != pEnd && *p == CharT('0'))
778cdf0e10cSrcweir {
779cdf0e10cSrcweir --nFracExp;
780cdf0e10cSrcweir ++p;
781cdf0e10cSrcweir }
782cdf0e10cSrcweir if ( nValExp == 0 )
783cdf0e10cSrcweir nValExp = nFracExp - 1; // no integer part => fraction exponent
784cdf0e10cSrcweir // one decimal digit needs ld(10) ~= 3.32 bits
785cdf0e10cSrcweir static const int nSigs = (DBL_MANT_DIG / 3) + 1;
786cdf0e10cSrcweir int nDigs = 0;
787cdf0e10cSrcweir for (; p != pEnd; ++p)
788cdf0e10cSrcweir {
789cdf0e10cSrcweir CharT c = *p;
790cdf0e10cSrcweir if (!isDigit(c))
791cdf0e10cSrcweir break;
792cdf0e10cSrcweir if ( nDigs < nSigs )
793cdf0e10cSrcweir { // further digits (more than nSigs) don't have any
794cdf0e10cSrcweir // significance
795cdf0e10cSrcweir fFrac = fFrac * 10.0 + static_cast<double>(c - CharT('0'));
796cdf0e10cSrcweir --nFracExp;
797cdf0e10cSrcweir ++nDigs;
798cdf0e10cSrcweir }
799cdf0e10cSrcweir }
800cdf0e10cSrcweir if ( fFrac != 0.0 )
801cdf0e10cSrcweir fVal += rtl::math::pow10Exp( fFrac, nFracExp );
802cdf0e10cSrcweir else if ( nValExp < 0 )
803cdf0e10cSrcweir nValExp = 0; // no digit other than 0 after decimal point
804cdf0e10cSrcweir }
805cdf0e10cSrcweir
806cdf0e10cSrcweir if ( nValExp > 0 )
807cdf0e10cSrcweir --nValExp; // started with offset +1 at the first mantissa digit
808cdf0e10cSrcweir
809cdf0e10cSrcweir // Exponent
810cdf0e10cSrcweir if (p != p0 && p != pEnd && (*p == CharT('E') || *p == CharT('e')))
811cdf0e10cSrcweir {
812cdf0e10cSrcweir ++p;
813cdf0e10cSrcweir bool bExpSign;
814cdf0e10cSrcweir if (p != pEnd && *p == CharT('-'))
815cdf0e10cSrcweir {
816cdf0e10cSrcweir bExpSign = true;
817cdf0e10cSrcweir ++p;
818cdf0e10cSrcweir }
819cdf0e10cSrcweir else
820cdf0e10cSrcweir {
821cdf0e10cSrcweir bExpSign = false;
822cdf0e10cSrcweir if (p != pEnd && *p == CharT('+'))
823cdf0e10cSrcweir ++p;
824cdf0e10cSrcweir }
825cdf0e10cSrcweir if ( fVal == 0.0 )
826cdf0e10cSrcweir { // no matter what follows, zero stays zero, but carry on the
827cdf0e10cSrcweir // offset
828cdf0e10cSrcweir while (p != pEnd && isDigit(*p))
829cdf0e10cSrcweir ++p;
830cdf0e10cSrcweir }
831cdf0e10cSrcweir else
832cdf0e10cSrcweir {
833cdf0e10cSrcweir bool bOverFlow = false;
834cdf0e10cSrcweir long nExp = 0;
835cdf0e10cSrcweir for (; p != pEnd; ++p)
836cdf0e10cSrcweir {
837cdf0e10cSrcweir CharT c = *p;
838cdf0e10cSrcweir if (!isDigit(c))
839cdf0e10cSrcweir break;
840cdf0e10cSrcweir int i = c - CharT('0');
841cdf0e10cSrcweir if ( long10Overflow( nExp, i ) )
842cdf0e10cSrcweir bOverFlow = true;
843cdf0e10cSrcweir else
844cdf0e10cSrcweir nExp = nExp * 10 + i;
845cdf0e10cSrcweir }
846cdf0e10cSrcweir if ( nExp )
847cdf0e10cSrcweir {
848cdf0e10cSrcweir if ( bExpSign )
849cdf0e10cSrcweir nExp = -nExp;
850cdf0e10cSrcweir long nAllExp = ( bOverFlow ? 0 : nExp + nValExp );
851cdf0e10cSrcweir if ( nAllExp > DBL_MAX_10_EXP || (bOverFlow && !bExpSign) )
852cdf0e10cSrcweir { // overflow
853cdf0e10cSrcweir fVal = HUGE_VAL;
854cdf0e10cSrcweir eStatus = rtl_math_ConversionStatus_OutOfRange;
855cdf0e10cSrcweir }
856cdf0e10cSrcweir else if ((nAllExp < DBL_MIN_10_EXP) ||
857cdf0e10cSrcweir (bOverFlow && bExpSign) )
858cdf0e10cSrcweir { // underflow
859cdf0e10cSrcweir fVal = 0.0;
860cdf0e10cSrcweir eStatus = rtl_math_ConversionStatus_OutOfRange;
861cdf0e10cSrcweir }
862cdf0e10cSrcweir else if ( nExp > DBL_MAX_10_EXP || nExp < DBL_MIN_10_EXP )
863cdf0e10cSrcweir { // compensate exponents
864cdf0e10cSrcweir fVal = rtl::math::pow10Exp( fVal, -nValExp );
865cdf0e10cSrcweir fVal = rtl::math::pow10Exp( fVal, nAllExp );
866cdf0e10cSrcweir }
867cdf0e10cSrcweir else
868cdf0e10cSrcweir fVal = rtl::math::pow10Exp( fVal, nExp ); // normal
869cdf0e10cSrcweir }
870cdf0e10cSrcweir }
871cdf0e10cSrcweir }
872cdf0e10cSrcweir else if (p - p0 == 2 && p != pEnd && p[0] == CharT('#')
873cdf0e10cSrcweir && p[-1] == cDecSeparator && p[-2] == CharT('1'))
874cdf0e10cSrcweir {
875cdf0e10cSrcweir if (pEnd - p >= 4 && p[1] == CharT('I') && p[2] == CharT('N')
876cdf0e10cSrcweir && p[3] == CharT('F'))
877cdf0e10cSrcweir {
878cdf0e10cSrcweir // "1.#INF", "+1.#INF", "-1.#INF"
879cdf0e10cSrcweir p += 4;
880cdf0e10cSrcweir fVal = HUGE_VAL;
881cdf0e10cSrcweir eStatus = rtl_math_ConversionStatus_OutOfRange;
882cdf0e10cSrcweir // Eat any further digits:
883cdf0e10cSrcweir while (p != pEnd && isDigit(*p))
884cdf0e10cSrcweir ++p;
885cdf0e10cSrcweir }
886cdf0e10cSrcweir else if (pEnd - p >= 4 && p[1] == CharT('N') && p[2] == CharT('A')
887cdf0e10cSrcweir && p[3] == CharT('N'))
888cdf0e10cSrcweir {
889cdf0e10cSrcweir // "1.#NAN", "+1.#NAN", "-1.#NAN"
890cdf0e10cSrcweir p += 4;
891cdf0e10cSrcweir rtl::math::setNan( &fVal );
892cdf0e10cSrcweir if (bSign)
893cdf0e10cSrcweir {
894cdf0e10cSrcweir union {
895cdf0e10cSrcweir double sd;
896cdf0e10cSrcweir sal_math_Double md;
897cdf0e10cSrcweir } m;
898cdf0e10cSrcweir m.sd = fVal;
899cdf0e10cSrcweir m.md.w32_parts.msw |= 0x80000000; // create negative NaN
900cdf0e10cSrcweir fVal = m.sd;
901cdf0e10cSrcweir bSign = false; // don't negate again
902cdf0e10cSrcweir }
903cdf0e10cSrcweir // Eat any further digits:
904cdf0e10cSrcweir while (p != pEnd && isDigit(*p))
905cdf0e10cSrcweir ++p;
906cdf0e10cSrcweir }
907cdf0e10cSrcweir }
908cdf0e10cSrcweir }
909cdf0e10cSrcweir
910cdf0e10cSrcweir // overflow also if more than DBL_MAX_10_EXP digits without decimal
911cdf0e10cSrcweir // separator, or 0. and more than DBL_MIN_10_EXP digits, ...
912cdf0e10cSrcweir bool bHuge = fVal == HUGE_VAL; // g++ 3.0.1 requires it this way...
913cdf0e10cSrcweir if ( bHuge )
914cdf0e10cSrcweir eStatus = rtl_math_ConversionStatus_OutOfRange;
915cdf0e10cSrcweir
916cdf0e10cSrcweir if ( bSign )
917cdf0e10cSrcweir fVal = -fVal;
918cdf0e10cSrcweir
919cdf0e10cSrcweir if (pStatus != 0)
920cdf0e10cSrcweir *pStatus = eStatus;
921cdf0e10cSrcweir if (pParsedEnd != 0)
922cdf0e10cSrcweir *pParsedEnd = p == p0 ? pBegin : p;
923cdf0e10cSrcweir
924cdf0e10cSrcweir return fVal;
925cdf0e10cSrcweir }
926cdf0e10cSrcweir
927cdf0e10cSrcweir }
928cdf0e10cSrcweir
rtl_math_stringToDouble(sal_Char const * pBegin,sal_Char const * pEnd,sal_Char cDecSeparator,sal_Char cGroupSeparator,rtl_math_ConversionStatus * pStatus,sal_Char const ** pParsedEnd)929cdf0e10cSrcweir double SAL_CALL rtl_math_stringToDouble(sal_Char const * pBegin,
930cdf0e10cSrcweir sal_Char const * pEnd,
931cdf0e10cSrcweir sal_Char cDecSeparator,
932cdf0e10cSrcweir sal_Char cGroupSeparator,
933cdf0e10cSrcweir rtl_math_ConversionStatus * pStatus,
934cdf0e10cSrcweir sal_Char const ** pParsedEnd)
935cdf0e10cSrcweir SAL_THROW_EXTERN_C()
936cdf0e10cSrcweir {
937cdf0e10cSrcweir return stringToDouble(pBegin, pEnd, cDecSeparator, cGroupSeparator, pStatus,
938cdf0e10cSrcweir pParsedEnd);
939cdf0e10cSrcweir }
940cdf0e10cSrcweir
rtl_math_uStringToDouble(sal_Unicode const * pBegin,sal_Unicode const * pEnd,sal_Unicode cDecSeparator,sal_Unicode cGroupSeparator,rtl_math_ConversionStatus * pStatus,sal_Unicode const ** pParsedEnd)941cdf0e10cSrcweir double SAL_CALL rtl_math_uStringToDouble(sal_Unicode const * pBegin,
942cdf0e10cSrcweir sal_Unicode const * pEnd,
943cdf0e10cSrcweir sal_Unicode cDecSeparator,
944cdf0e10cSrcweir sal_Unicode cGroupSeparator,
945cdf0e10cSrcweir rtl_math_ConversionStatus * pStatus,
946cdf0e10cSrcweir sal_Unicode const ** pParsedEnd)
947cdf0e10cSrcweir SAL_THROW_EXTERN_C()
948cdf0e10cSrcweir {
949cdf0e10cSrcweir return stringToDouble(pBegin, pEnd, cDecSeparator, cGroupSeparator, pStatus,
950cdf0e10cSrcweir pParsedEnd);
951cdf0e10cSrcweir }
952cdf0e10cSrcweir
rtl_math_round(double fValue,int nDecPlaces,enum rtl_math_RoundingMode eMode)953cdf0e10cSrcweir double SAL_CALL rtl_math_round(double fValue, int nDecPlaces,
954cdf0e10cSrcweir enum rtl_math_RoundingMode eMode)
955cdf0e10cSrcweir SAL_THROW_EXTERN_C()
956cdf0e10cSrcweir {
957cdf0e10cSrcweir OSL_ASSERT(nDecPlaces >= -20 && nDecPlaces <= 20);
958cdf0e10cSrcweir
959cdf0e10cSrcweir if ( fValue == 0.0 )
960cdf0e10cSrcweir return fValue;
961cdf0e10cSrcweir
962cdf0e10cSrcweir // sign adjustment
963cdf0e10cSrcweir bool bSign = rtl::math::isSignBitSet( fValue );
964cdf0e10cSrcweir if ( bSign )
965cdf0e10cSrcweir fValue = -fValue;
966cdf0e10cSrcweir
967cdf0e10cSrcweir double fFac = 0;
968cdf0e10cSrcweir if ( nDecPlaces != 0 )
969cdf0e10cSrcweir {
970cdf0e10cSrcweir // max 20 decimals, we don't have unlimited precision
971cdf0e10cSrcweir // #38810# and no overflow on fValue*=fFac
972cdf0e10cSrcweir if ( nDecPlaces < -20 || 20 < nDecPlaces || fValue > (DBL_MAX / 1e20) )
973cdf0e10cSrcweir return bSign ? -fValue : fValue;
974cdf0e10cSrcweir
975cdf0e10cSrcweir fFac = getN10Exp( nDecPlaces );
976cdf0e10cSrcweir fValue *= fFac;
977cdf0e10cSrcweir }
978cdf0e10cSrcweir //else //! uninitialized fFac, not needed
979cdf0e10cSrcweir
980cdf0e10cSrcweir switch ( eMode )
981cdf0e10cSrcweir {
982cdf0e10cSrcweir case rtl_math_RoundingMode_Corrected :
983cdf0e10cSrcweir {
984cdf0e10cSrcweir int nExp; // exponent for correction
985cdf0e10cSrcweir if ( fValue > 0.0 )
986cdf0e10cSrcweir nExp = static_cast<int>( floor( log10( fValue ) ) );
987cdf0e10cSrcweir else
988cdf0e10cSrcweir nExp = 0;
989cdf0e10cSrcweir int nIndex = 15 - nExp;
990cdf0e10cSrcweir if ( nIndex > 15 )
991cdf0e10cSrcweir nIndex = 15;
992cdf0e10cSrcweir else if ( nIndex <= 1 )
993cdf0e10cSrcweir nIndex = 0;
994cdf0e10cSrcweir fValue = floor( fValue + 0.5 + nKorrVal[nIndex] );
995cdf0e10cSrcweir }
996cdf0e10cSrcweir break;
997cdf0e10cSrcweir case rtl_math_RoundingMode_Down :
998cdf0e10cSrcweir fValue = rtl::math::approxFloor( fValue );
999cdf0e10cSrcweir break;
1000cdf0e10cSrcweir case rtl_math_RoundingMode_Up :
1001cdf0e10cSrcweir fValue = rtl::math::approxCeil( fValue );
1002cdf0e10cSrcweir break;
1003cdf0e10cSrcweir case rtl_math_RoundingMode_Floor :
1004cdf0e10cSrcweir fValue = bSign ? rtl::math::approxCeil( fValue )
1005cdf0e10cSrcweir : rtl::math::approxFloor( fValue );
1006cdf0e10cSrcweir break;
1007cdf0e10cSrcweir case rtl_math_RoundingMode_Ceiling :
1008cdf0e10cSrcweir fValue = bSign ? rtl::math::approxFloor( fValue )
1009cdf0e10cSrcweir : rtl::math::approxCeil( fValue );
1010cdf0e10cSrcweir break;
1011cdf0e10cSrcweir case rtl_math_RoundingMode_HalfDown :
1012cdf0e10cSrcweir {
1013cdf0e10cSrcweir double f = floor( fValue );
1014cdf0e10cSrcweir fValue = ((fValue - f) <= 0.5) ? f : ceil( fValue );
1015cdf0e10cSrcweir }
1016cdf0e10cSrcweir break;
1017cdf0e10cSrcweir case rtl_math_RoundingMode_HalfUp :
1018cdf0e10cSrcweir {
1019cdf0e10cSrcweir double f = floor( fValue );
1020cdf0e10cSrcweir fValue = ((fValue - f) < 0.5) ? f : ceil( fValue );
1021cdf0e10cSrcweir }
1022cdf0e10cSrcweir break;
1023cdf0e10cSrcweir case rtl_math_RoundingMode_HalfEven :
1024cdf0e10cSrcweir #if defined FLT_ROUNDS
1025cdf0e10cSrcweir /*
1026cdf0e10cSrcweir Use fast version. FLT_ROUNDS may be defined to a function by some compilers!
1027cdf0e10cSrcweir
1028cdf0e10cSrcweir DBL_EPSILON is the smallest fractional number which can be represented,
1029cdf0e10cSrcweir its reciprocal is therefore the smallest number that cannot have a
1030cdf0e10cSrcweir fractional part. Once you add this reciprocal to `x', its fractional part
1031cdf0e10cSrcweir is stripped off. Simply subtracting the reciprocal back out returns `x'
1032cdf0e10cSrcweir without its fractional component.
1033cdf0e10cSrcweir Simple, clever, and elegant - thanks to Ross Cottrell, the original author,
1034cdf0e10cSrcweir who placed it into public domain.
1035cdf0e10cSrcweir
1036cdf0e10cSrcweir volatile: prevent compiler from being too smart
1037cdf0e10cSrcweir */
1038cdf0e10cSrcweir if ( FLT_ROUNDS == 1 )
1039cdf0e10cSrcweir {
1040cdf0e10cSrcweir volatile double x = fValue + 1.0 / DBL_EPSILON;
1041cdf0e10cSrcweir fValue = x - 1.0 / DBL_EPSILON;
1042cdf0e10cSrcweir }
1043cdf0e10cSrcweir else
1044cdf0e10cSrcweir #endif // FLT_ROUNDS
1045cdf0e10cSrcweir {
1046cdf0e10cSrcweir double f = floor( fValue );
1047cdf0e10cSrcweir if ( (fValue - f) != 0.5 )
1048cdf0e10cSrcweir fValue = floor( fValue + 0.5 );
1049cdf0e10cSrcweir else
1050cdf0e10cSrcweir {
1051cdf0e10cSrcweir double g = f / 2.0;
1052cdf0e10cSrcweir fValue = (g == floor( g )) ? f : (f + 1.0);
1053cdf0e10cSrcweir }
1054cdf0e10cSrcweir }
1055cdf0e10cSrcweir break;
1056cdf0e10cSrcweir default:
1057cdf0e10cSrcweir OSL_ASSERT(false);
1058cdf0e10cSrcweir break;
1059cdf0e10cSrcweir }
1060cdf0e10cSrcweir
1061cdf0e10cSrcweir if ( nDecPlaces != 0 )
1062cdf0e10cSrcweir fValue /= fFac;
1063cdf0e10cSrcweir
1064cdf0e10cSrcweir return bSign ? -fValue : fValue;
1065cdf0e10cSrcweir }
1066cdf0e10cSrcweir
1067cdf0e10cSrcweir
rtl_math_pow10Exp(double fValue,int nExp)1068cdf0e10cSrcweir double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THROW_EXTERN_C()
1069cdf0e10cSrcweir {
1070cdf0e10cSrcweir return fValue * getN10Exp( nExp );
1071cdf0e10cSrcweir }
1072cdf0e10cSrcweir
1073cdf0e10cSrcweir
rtl_math_approxValue(double fValue)1074cdf0e10cSrcweir double SAL_CALL rtl_math_approxValue( double fValue ) SAL_THROW_EXTERN_C()
1075cdf0e10cSrcweir {
1076cdf0e10cSrcweir if (fValue == 0.0 || fValue == HUGE_VAL || !::rtl::math::isFinite( fValue))
1077cdf0e10cSrcweir // We don't handle these conditions. Bail out.
1078cdf0e10cSrcweir return fValue;
1079cdf0e10cSrcweir
1080cdf0e10cSrcweir double fOrigValue = fValue;
1081cdf0e10cSrcweir
1082cdf0e10cSrcweir bool bSign = ::rtl::math::isSignBitSet( fValue);
1083cdf0e10cSrcweir if (bSign)
1084cdf0e10cSrcweir fValue = -fValue;
1085cdf0e10cSrcweir
1086cdf0e10cSrcweir int nExp = static_cast<int>( floor( log10( fValue)));
1087cdf0e10cSrcweir nExp = 14 - nExp;
1088cdf0e10cSrcweir double fExpValue = getN10Exp( nExp);
1089cdf0e10cSrcweir
1090cdf0e10cSrcweir fValue *= fExpValue;
1091cdf0e10cSrcweir // If the original value was near DBL_MIN we got an overflow. Restore and
1092cdf0e10cSrcweir // bail out.
1093cdf0e10cSrcweir if (!rtl::math::isFinite( fValue))
1094cdf0e10cSrcweir return fOrigValue;
1095cdf0e10cSrcweir fValue = rtl_math_round( fValue, 0, rtl_math_RoundingMode_Corrected);
1096cdf0e10cSrcweir fValue /= fExpValue;
1097cdf0e10cSrcweir // If the original value was near DBL_MAX we got an overflow. Restore and
1098cdf0e10cSrcweir // bail out.
1099cdf0e10cSrcweir if (!rtl::math::isFinite( fValue))
1100cdf0e10cSrcweir return fOrigValue;
1101cdf0e10cSrcweir
1102cdf0e10cSrcweir return bSign ? -fValue : fValue;
1103cdf0e10cSrcweir }
1104cdf0e10cSrcweir
1105cdf0e10cSrcweir
rtl_math_expm1(double fValue)1106cdf0e10cSrcweir double SAL_CALL rtl_math_expm1( double fValue ) SAL_THROW_EXTERN_C()
1107cdf0e10cSrcweir {
1108cdf0e10cSrcweir double fe = exp( fValue );
1109cdf0e10cSrcweir if (fe == 1.0)
1110cdf0e10cSrcweir return fValue;
1111cdf0e10cSrcweir if (fe-1.0 == -1.0)
1112cdf0e10cSrcweir return -1.0;
1113cdf0e10cSrcweir return (fe-1.0) * fValue / log(fe);
1114cdf0e10cSrcweir }
1115cdf0e10cSrcweir
rtl_math_powr(double fValue,double fExp)1116*cda7f8b3SPedro Giffuni double SAL_CALL rtl_math_powr( double fValue, double fExp ) SAL_THROW_EXTERN_C()
1117*cda7f8b3SPedro Giffuni {
1118*cda7f8b3SPedro Giffuni if ((fValue == 0.0 && fExp == 0.0) ||
1119*cda7f8b3SPedro Giffuni (rtl::math::isInf( fExp ) && !rtl::math::isSignBitSet( fExp )) ||
1120*cda7f8b3SPedro Giffuni (rtl::math::isInf( fValue ) && !rtl::math::isSignBitSet( fValue )))
1121*cda7f8b3SPedro Giffuni {
1122*cda7f8b3SPedro Giffuni double fResult;
1123*cda7f8b3SPedro Giffuni ::rtl::math::setNan( &fResult );
1124*cda7f8b3SPedro Giffuni return fResult;
1125*cda7f8b3SPedro Giffuni }
1126*cda7f8b3SPedro Giffuni return pow(fValue, fExp);
1127*cda7f8b3SPedro Giffuni }
1128*cda7f8b3SPedro Giffuni
1129cdf0e10cSrcweir
rtl_math_log1p(double fValue)1130cdf0e10cSrcweir double SAL_CALL rtl_math_log1p( double fValue ) SAL_THROW_EXTERN_C()
1131cdf0e10cSrcweir {
1132cdf0e10cSrcweir // Use volatile because a compiler may be too smart "optimizing" the
1133cdf0e10cSrcweir // condition such that in certain cases the else path was called even if
1134cdf0e10cSrcweir // (fp==1.0) was true, where the term (fp-1.0) then resulted in 0.0 and
1135cdf0e10cSrcweir // hence the entire expression resulted in NaN.
1136cdf0e10cSrcweir // Happened with g++ 3.4.1 and an input value of 9.87E-18
1137cdf0e10cSrcweir volatile double fp = 1.0 + fValue;
1138cdf0e10cSrcweir if (fp == 1.0)
1139cdf0e10cSrcweir return fValue;
1140cdf0e10cSrcweir else
1141cdf0e10cSrcweir return log(fp) * fValue / (fp-1.0);
1142cdf0e10cSrcweir }
1143cdf0e10cSrcweir
1144cdf0e10cSrcweir
rtl_math_atanh(double fValue)1145cdf0e10cSrcweir double SAL_CALL rtl_math_atanh( double fValue ) SAL_THROW_EXTERN_C()
1146cdf0e10cSrcweir {
1147cdf0e10cSrcweir return 0.5 * rtl_math_log1p( 2.0 * fValue / (1.0-fValue) );
1148cdf0e10cSrcweir }
1149cdf0e10cSrcweir
1150cdf0e10cSrcweir
1151cdf0e10cSrcweir /** Parent error function (erf) that calls different algorithms based on the
1152cdf0e10cSrcweir value of x. It takes care of cases where x is negative as erf is an odd
1153cdf0e10cSrcweir function i.e. erf(-x) = -erf(x).
1154cdf0e10cSrcweir
1155cdf0e10cSrcweir Kramer, W., and Blomquist, F., 2000, Algorithms with Guaranteed Error Bounds
1156cdf0e10cSrcweir for the Error Function and the Complementary Error Function
1157cdf0e10cSrcweir
1158cdf0e10cSrcweir http://www.math.uni-wuppertal.de/wrswt/literatur_en.html
1159cdf0e10cSrcweir
1160cdf0e10cSrcweir @author Kohei Yoshida <kohei@openoffice.org>
1161cdf0e10cSrcweir
1162cdf0e10cSrcweir @see #i55735#
1163cdf0e10cSrcweir */
rtl_math_erf(double x)1164cdf0e10cSrcweir double SAL_CALL rtl_math_erf( double x ) SAL_THROW_EXTERN_C()
1165cdf0e10cSrcweir {
1166cdf0e10cSrcweir if( x == 0.0 )
1167cdf0e10cSrcweir return 0.0;
1168cdf0e10cSrcweir
1169cdf0e10cSrcweir bool bNegative = false;
1170cdf0e10cSrcweir if ( x < 0.0 )
1171cdf0e10cSrcweir {
1172cdf0e10cSrcweir x = fabs( x );
1173cdf0e10cSrcweir bNegative = true;
1174cdf0e10cSrcweir }
1175cdf0e10cSrcweir
1176cdf0e10cSrcweir double fErf = 1.0;
1177cdf0e10cSrcweir if ( x < 1.0e-10 )
1178cdf0e10cSrcweir fErf = (double) (x*1.1283791670955125738961589031215452L);
1179cdf0e10cSrcweir else if ( x < 0.65 )
1180cdf0e10cSrcweir lcl_Erf0065( x, fErf );
1181cdf0e10cSrcweir else
1182cdf0e10cSrcweir fErf = 1.0 - rtl_math_erfc( x );
1183cdf0e10cSrcweir
1184cdf0e10cSrcweir if ( bNegative )
1185cdf0e10cSrcweir fErf *= -1.0;
1186cdf0e10cSrcweir
1187cdf0e10cSrcweir return fErf;
1188cdf0e10cSrcweir }
1189cdf0e10cSrcweir
1190cdf0e10cSrcweir
1191cdf0e10cSrcweir /** Parent complementary error function (erfc) that calls different algorithms
1192cdf0e10cSrcweir based on the value of x. It takes care of cases where x is negative as erfc
1193cdf0e10cSrcweir satisfies relationship erfc(-x) = 2 - erfc(x). See the comment for Erf(x)
1194cdf0e10cSrcweir for the source publication.
1195cdf0e10cSrcweir
1196cdf0e10cSrcweir @author Kohei Yoshida <kohei@openoffice.org>
1197cdf0e10cSrcweir
1198cdf0e10cSrcweir @see #i55735#, moved from module scaddins (#i97091#)
1199cdf0e10cSrcweir
1200cdf0e10cSrcweir */
rtl_math_erfc(double x)1201cdf0e10cSrcweir double SAL_CALL rtl_math_erfc( double x ) SAL_THROW_EXTERN_C()
1202cdf0e10cSrcweir {
1203cdf0e10cSrcweir if ( x == 0.0 )
1204cdf0e10cSrcweir return 1.0;
1205cdf0e10cSrcweir
1206cdf0e10cSrcweir bool bNegative = false;
1207cdf0e10cSrcweir if ( x < 0.0 )
1208cdf0e10cSrcweir {
1209cdf0e10cSrcweir x = fabs( x );
1210cdf0e10cSrcweir bNegative = true;
1211cdf0e10cSrcweir }
1212cdf0e10cSrcweir
1213cdf0e10cSrcweir double fErfc = 0.0;
1214cdf0e10cSrcweir if ( x >= 0.65 )
1215cdf0e10cSrcweir {
1216cdf0e10cSrcweir if ( x < 6.0 )
1217cdf0e10cSrcweir lcl_Erfc0600( x, fErfc );
1218cdf0e10cSrcweir else
1219cdf0e10cSrcweir lcl_Erfc2654( x, fErfc );
1220cdf0e10cSrcweir }
1221cdf0e10cSrcweir else
1222cdf0e10cSrcweir fErfc = 1.0 - rtl_math_erf( x );
1223cdf0e10cSrcweir
1224cdf0e10cSrcweir if ( bNegative )
1225cdf0e10cSrcweir fErfc = 2.0 - fErfc;
1226cdf0e10cSrcweir
1227cdf0e10cSrcweir return fErfc;
1228cdf0e10cSrcweir }
1229cdf0e10cSrcweir
1230cdf0e10cSrcweir /** improved accuracy of asinh for |x| large and for x near zero
1231cdf0e10cSrcweir @see #i97605#
1232cdf0e10cSrcweir */
rtl_math_asinh(double fX)1233cdf0e10cSrcweir double SAL_CALL rtl_math_asinh( double fX ) SAL_THROW_EXTERN_C()
1234cdf0e10cSrcweir {
1235cdf0e10cSrcweir double fSign = 1.0;
1236cdf0e10cSrcweir if ( fX == 0.0 )
1237cdf0e10cSrcweir return 0.0;
1238cdf0e10cSrcweir else
1239cdf0e10cSrcweir {
1240cdf0e10cSrcweir if ( fX < 0.0 )
1241cdf0e10cSrcweir {
1242cdf0e10cSrcweir fX = - fX;
1243cdf0e10cSrcweir fSign = -1.0;
1244cdf0e10cSrcweir }
1245cdf0e10cSrcweir if ( fX < 0.125 )
1246cdf0e10cSrcweir return fSign * rtl_math_log1p( fX + fX*fX / (1.0 + sqrt( 1.0 + fX*fX)));
1247cdf0e10cSrcweir else if ( fX < 1.25e7 )
1248cdf0e10cSrcweir return fSign * log( fX + sqrt( 1.0 + fX*fX));
1249cdf0e10cSrcweir else
1250cdf0e10cSrcweir return fSign * log( 2.0*fX);
1251cdf0e10cSrcweir }
1252cdf0e10cSrcweir }
1253cdf0e10cSrcweir
1254cdf0e10cSrcweir /** improved accuracy of acosh for x large and for x near 1
1255cdf0e10cSrcweir @see #i97605#
1256cdf0e10cSrcweir */
rtl_math_acosh(double fX)1257cdf0e10cSrcweir double SAL_CALL rtl_math_acosh( double fX ) SAL_THROW_EXTERN_C()
1258cdf0e10cSrcweir {
1259cdf0e10cSrcweir volatile double fZ = fX - 1.0;
1260cdf0e10cSrcweir if ( fX < 1.0 )
1261cdf0e10cSrcweir {
1262cdf0e10cSrcweir double fResult;
1263cdf0e10cSrcweir ::rtl::math::setNan( &fResult );
1264cdf0e10cSrcweir return fResult;
1265cdf0e10cSrcweir }
1266cdf0e10cSrcweir else if ( fX == 1.0 )
1267cdf0e10cSrcweir return 0.0;
1268cdf0e10cSrcweir else if ( fX < 1.1 )
1269cdf0e10cSrcweir return rtl_math_log1p( fZ + sqrt( fZ*fZ + 2.0*fZ));
1270cdf0e10cSrcweir else if ( fX < 1.25e7 )
1271cdf0e10cSrcweir return log( fX + sqrt( fX*fX - 1.0));
1272cdf0e10cSrcweir else
1273cdf0e10cSrcweir return log( 2.0*fX);
1274cdf0e10cSrcweir }
1275