1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
3 * This file is part of the LibreOffice project.
5 * This Source Code Form is subject to the terms of the Mozilla Public
6 * License, v. 2.0. If a copy of the MPL was not distributed with this
7 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 * This file incorporates work covered by the following license notice:
11 * Licensed to the Apache Software Foundation (ASF) under one or more
12 * contributor license agreements. See the NOTICE file distributed
13 * with this work for additional information regarding copyright
14 * ownership. The ASF licenses this file to you under the Apache
15 * License, Version 2.0 (the "License"); you may not use this file
16 * except in compliance with the License. You may obtain a copy of
17 * the License at http://www.apache.org/licenses/LICENSE-2.0 .
22 #include "osl/diagnose.h"
23 #include "rtl/alloc.h"
24 #include "rtl/character.hxx"
25 #include "rtl/math.hxx"
26 #include "rtl/strbuf.h"
27 #include "rtl/string.h"
28 #include "rtl/ustrbuf.h"
29 #include "rtl/ustring.h"
30 #include "sal/mathconf.h"
31 #include "sal/types.h"
40 static int const n10Count
= 16;
41 static double const n10s
[2][n10Count
] = {
42 { 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
,
43 1e9
, 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
},
44 { 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8,
45 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16 }
48 // return pow(10.0,nExp) optimized for exponents in the interval [-16,16]
49 static double getN10Exp( int nExp
)
53 // && -nExp > 0 necessary for std::numeric_limits<int>::min()
54 // because -nExp = nExp
55 if ( -nExp
<= n10Count
&& -nExp
> 0 )
56 return n10s
[1][-nExp
-1];
58 return pow( 10.0, static_cast<double>( nExp
) );
62 if ( nExp
<= n10Count
)
63 return n10s
[0][nExp
-1];
65 return pow( 10.0, static_cast<double>( nExp
) );
71 /** Approximation algorithm for erf for 0 < x < 0.65. */
72 static void lcl_Erf0065( double x
, double& fVal
)
74 static const double pn
[] = {
76 1.35894887627277916E-1,
77 4.03259488531795274E-2,
78 1.20339380863079457E-3,
79 6.49254556481904354E-5
81 static const double qn
[] = {
83 4.53767041780002545E-1,
84 8.69936222615385890E-2,
85 8.49717371168693357E-3,
86 3.64915280629351082E-4
91 for ( unsigned int i
= 0; i
<= 4; ++i
)
97 fVal
= x
* fPSum
/ fQSum
;
100 /** Approximation algorithm for erfc for 0.65 < x < 6.0. */
101 static void lcl_Erfc0600( double x
, double& fVal
)
111 static const double pn22
[] = {
112 9.99999992049799098E-1,
114 8.78115804155881782E-1,
115 3.31899559578213215E-1,
116 7.14193832506776067E-2,
117 7.06940843763253131E-3
119 static const double qn22
[] = {
124 5.94651311286481502E-1,
125 1.26579413030177940E-1,
126 1.25304936549413393E-2
131 else /* if ( x < 6.0 ) this is true, but the compiler does not know */
133 static const double pn60
[] = {
134 9.99921140009714409E-1,
137 5.81528574177741135E-1,
138 1.57289620742838702E-1,
139 2.25716982919217555E-2
141 static const double qn60
[] = {
147 2.78788439273628983E-1,
148 4.00072964526861362E-2
154 for ( unsigned int i
= 0; i
< 6; ++i
)
156 fPSum
+= pn
[i
]*fXPow
;
157 fQSum
+= qn
[i
]*fXPow
;
160 fQSum
+= qn
[6]*fXPow
;
161 fVal
= exp( -1.0*x
*x
)* fPSum
/ fQSum
;
164 /** Approximation algorithm for erfc for 6.0 < x < 26.54 (but used for all
166 static void lcl_Erfc2654( double x
, double& fVal
)
168 static const double pn
[] = {
169 5.64189583547756078E-1,
171 3.84683103716117320E1
,
172 4.77209965874436377E1
,
175 static const double qn
[] = {
177 1.61020914205869003E1
,
178 7.54843505665954743E1
,
179 1.12123870801026015E2
,
180 3.73997570145040850E1
187 for ( unsigned int i
= 0; i
<= 4; ++i
)
189 fPSum
+= pn
[i
]*fXPow
;
190 fQSum
+= qn
[i
]*fXPow
;
193 fVal
= exp(-1.0*x
*x
)*fPSum
/ (x
*fQSum
);
198 double const nKorrVal
[] = {
199 0, 9e-1, 9e-2, 9e-3, 9e-4, 9e-5, 9e-6, 9e-7, 9e-8,
200 9e-9, 9e-10, 9e-11, 9e-12, 9e-13, 9e-14, 9e-15
205 typedef sal_Char Char
;
207 typedef rtl_String String
;
209 static inline void createString(rtl_String
** pString
,
210 sal_Char
const * pChars
, sal_Int32 nLen
)
212 rtl_string_newFromStr_WithLength(pString
, pChars
, nLen
);
215 static inline void createBuffer(rtl_String
** pBuffer
,
216 sal_Int32
* pCapacity
)
218 rtl_string_new_WithLength(pBuffer
, *pCapacity
);
221 static inline void appendChars(rtl_String
** pBuffer
, sal_Int32
* pCapacity
,
222 sal_Int32
* pOffset
, sal_Char
const * pChars
,
225 assert(pChars
!= nullptr);
226 rtl_stringbuffer_insert(pBuffer
, pCapacity
, *pOffset
, pChars
, nLen
);
230 static inline void appendAscii(rtl_String
** pBuffer
, sal_Int32
* pCapacity
,
231 sal_Int32
* pOffset
, sal_Char
const * pStr
,
234 assert(pStr
!= nullptr);
235 rtl_stringbuffer_insert(pBuffer
, pCapacity
, *pOffset
, pStr
, nLen
);
242 typedef sal_Unicode Char
;
244 typedef rtl_uString String
;
246 static inline void createString(rtl_uString
** pString
,
247 sal_Unicode
const * pChars
, sal_Int32 nLen
)
249 rtl_uString_newFromStr_WithLength(pString
, pChars
, nLen
);
252 static inline void createBuffer(rtl_uString
** pBuffer
,
253 sal_Int32
* pCapacity
)
255 rtl_uString_new_WithLength(pBuffer
, *pCapacity
);
258 static inline void appendChars(rtl_uString
** pBuffer
,
259 sal_Int32
* pCapacity
, sal_Int32
* pOffset
,
260 sal_Unicode
const * pChars
, sal_Int32 nLen
)
262 assert(pChars
!= nullptr);
263 rtl_uStringbuffer_insert(pBuffer
, pCapacity
, *pOffset
, pChars
, nLen
);
267 static inline void appendAscii(rtl_uString
** pBuffer
,
268 sal_Int32
* pCapacity
, sal_Int32
* pOffset
,
269 sal_Char
const * pStr
, sal_Int32 nLen
)
271 rtl_uStringbuffer_insert_ascii(pBuffer
, pCapacity
, *pOffset
, pStr
,
277 // Solaris C++ 5.2 compiler has problems when "StringT ** pResult" is
278 // "typename T::String ** pResult" instead:
279 template< typename T
, typename StringT
>
280 inline void doubleToString(StringT
** pResult
,
281 sal_Int32
* pResultCapacity
, sal_Int32 nResultOffset
,
282 double fValue
, rtl_math_StringFormat eFormat
,
283 sal_Int32 nDecPlaces
, typename
T::Char cDecSeparator
,
284 sal_Int32
const * pGroups
,
285 typename
T::Char cGroupSeparator
,
286 bool bEraseTrailingDecZeros
)
288 static double const nRoundVal
[] = {
289 5.0e+0, 0.5e+0, 0.5e-1, 0.5e-2, 0.5e-3, 0.5e-4, 0.5e-5, 0.5e-6,
290 0.5e-7, 0.5e-8, 0.5e-9, 0.5e-10,0.5e-11,0.5e-12,0.5e-13,0.5e-14
293 // sign adjustment, instead of testing for fValue<0.0 this will also fetch
295 bool bSign
= rtl::math::isSignBitSet( fValue
);
299 if ( rtl::math::isNan( fValue
) )
301 // #i112652# XMLSchema-2
302 sal_Int32 nCapacity
= RTL_CONSTASCII_LENGTH("NaN");
303 if (pResultCapacity
== 0)
305 pResultCapacity
= &nCapacity
;
306 T::createBuffer(pResult
, pResultCapacity
);
309 T::appendAscii(pResult
, pResultCapacity
, &nResultOffset
,
310 RTL_CONSTASCII_STRINGPARAM("NaN"));
315 bool bHuge
= fValue
== HUGE_VAL
; // g++ 3.0.1 requires it this way...
316 if ( bHuge
|| rtl::math::isInf( fValue
) )
318 // #i112652# XMLSchema-2
319 sal_Int32 nCapacity
= RTL_CONSTASCII_LENGTH("-INF");
320 if (pResultCapacity
== 0)
322 pResultCapacity
= &nCapacity
;
323 T::createBuffer(pResult
, pResultCapacity
);
327 T::appendAscii(pResult
, pResultCapacity
, &nResultOffset
,
328 RTL_CONSTASCII_STRINGPARAM("-"));
329 T::appendAscii(pResult
, pResultCapacity
, &nResultOffset
,
330 RTL_CONSTASCII_STRINGPARAM("INF"));
339 nExp
= static_cast< int >( floor( log10( fValue
) ) );
340 fValue
/= getN10Exp( nExp
);
345 case rtl_math_StringFormat_Automatic
:
346 { // E or F depending on exponent magnitude
348 if ( nExp
<= -15 || nExp
>= 15 ) // #58531# was <-16, >16
351 eFormat
= rtl_math_StringFormat_E
;
357 nPrec
= 15 - nExp
- 1;
358 eFormat
= rtl_math_StringFormat_F
;
363 eFormat
= rtl_math_StringFormat_F
;
366 if ( nDecPlaces
== rtl_math_DecimalPlaces_Max
)
370 case rtl_math_StringFormat_G
:
371 { // G-Point, similar to sprintf %G
372 if ( nDecPlaces
== rtl_math_DecimalPlaces_DefaultSignificance
)
374 if ( nExp
< -4 || nExp
>= nDecPlaces
)
376 nDecPlaces
= std::max
< sal_Int32
>( 1, nDecPlaces
- 1 );
377 eFormat
= rtl_math_StringFormat_E
;
381 nDecPlaces
= std::max
< sal_Int32
>( 0, nDecPlaces
- nExp
- 1 );
382 eFormat
= rtl_math_StringFormat_F
;
390 sal_Int32 nDigits
= nDecPlaces
+ 1;
392 if( eFormat
== rtl_math_StringFormat_F
)
398 if( ( fValue
+= nRoundVal
[ nDigits
> 15 ? 15 : nDigits
] ) >= 10 )
402 if( eFormat
== rtl_math_StringFormat_F
)
407 static sal_Int32
const nBufMax
= 256;
408 typename
T::Char aBuf
[nBufMax
];
409 typename
T::Char
* pBuf
;
410 sal_Int32 nBuf
= static_cast< sal_Int32
>
411 ( nDigits
<= 0 ? std::max
< sal_Int32
>( nDecPlaces
, abs(nExp
) )
412 : nDigits
+ nDecPlaces
) + 10 + (pGroups
? abs(nDigits
) * 2 : 0);
413 if ( nBuf
> nBufMax
)
415 pBuf
= reinterpret_cast< typename
T::Char
* >(
416 rtl_allocateMemory(nBuf
* sizeof (typename
T::Char
)));
417 OSL_ENSURE(pBuf
!= 0, "Out of memory");
421 typename
T::Char
* p
= pBuf
;
423 *p
++ = static_cast< typename
T::Char
>('-');
425 bool bHasDec
= false;
428 // Check for F format and number < 1
429 if( eFormat
== rtl_math_StringFormat_F
)
433 *p
++ = static_cast< typename
T::Char
>('0');
434 if ( nDecPlaces
> 0 )
436 *p
++ = cDecSeparator
;
439 sal_Int32 i
= ( nDigits
<= 0 ? nDecPlaces
: -nExp
- 1 );
441 *p
++ = static_cast< typename
T::Char
>('0');
450 int nGrouping
= 0, nGroupSelector
= 0, nGroupExceed
= 0;
451 if ( nDecPos
> 1 && pGroups
&& pGroups
[0] && cGroupSeparator
)
453 while ( nGrouping
+ pGroups
[nGroupSelector
] < nDecPos
)
455 nGrouping
+= pGroups
[ nGroupSelector
];
456 if ( pGroups
[nGroupSelector
+1] )
458 if ( nGrouping
+ pGroups
[nGroupSelector
+1] >= nDecPos
)
462 else if ( !nGroupExceed
)
463 nGroupExceed
= nGrouping
;
470 for ( int i
= 0; ; i
++ )
475 if (nDigits
-1 == 0 && i
> 0 && i
< 14)
476 nDigit
= static_cast< int >( floor( fValue
477 + nKorrVal
[15-i
] ) );
479 nDigit
= static_cast< int >( fValue
+ 1E-15 );
481 { // after-treatment of up-rounding to the next decade
482 sal_Int32 sLen
= static_cast< long >(p
-pBuf
)-1;
486 if ( eFormat
== rtl_math_StringFormat_F
)
488 *p
++ = static_cast< typename
T::Char
>('1');
489 *p
++ = static_cast< typename
T::Char
>('0');
493 *p
++ = static_cast< typename
T::Char
>('1');
494 *p
++ = cDecSeparator
;
495 *p
++ = static_cast< typename
T::Char
>('0');
502 for (sal_Int32 j
= sLen
; j
>= 0; j
--)
504 typename
T::Char cS
= pBuf
[j
];
505 if (cS
!= cDecSeparator
)
507 if ( cS
!= static_cast< typename
T::Char
>('9'))
510 j
= -1; // break loop
515 = static_cast< typename
T::Char
>('0');
518 if ( eFormat
== rtl_math_StringFormat_F
)
520 typename
T::Char
* px
= p
++;
526 pBuf
[0] = static_cast<
527 typename
T::Char
>('1');
531 pBuf
[j
] = static_cast<
532 typename
T::Char
>('1');
539 *p
++ = static_cast< typename
T::Char
>('0');
545 *p
++ = static_cast< typename
T::Char
>(
546 nDigit
+ static_cast< typename
T::Char
>('0') );
547 fValue
= ( fValue
- nDigit
) * 10.0;
551 *p
++ = static_cast< typename
T::Char
>('0');
558 *p
++ = cDecSeparator
;
561 else if ( nDecPos
== nGrouping
)
563 *p
++ = cGroupSeparator
;
564 nGrouping
-= pGroups
[ nGroupSelector
];
565 if ( nGroupSelector
&& nGrouping
< nGroupExceed
)
572 if ( !bHasDec
&& eFormat
== rtl_math_StringFormat_F
)
573 { // nDecPlaces < 0 did round the value
574 while ( --nDecPos
> 0 )
575 { // fill before decimal point
576 if ( nDecPos
== nGrouping
)
578 *p
++ = cGroupSeparator
;
579 nGrouping
-= pGroups
[ nGroupSelector
];
580 if ( nGroupSelector
&& nGrouping
< nGroupExceed
)
583 *p
++ = static_cast< typename
T::Char
>('0');
587 if ( bEraseTrailingDecZeros
&& bHasDec
&& p
> pBuf
)
589 while ( *(p
-1) == static_cast< typename
T::Char
>('0') )
591 if ( *(p
-1) == cDecSeparator
)
595 // Print the exponent ('E', followed by '+' or '-', followed by exactly
596 // three digits). The code in rtl_[u]str_valueOf{Float|Double} relies on
598 if( eFormat
== rtl_math_StringFormat_E
)
601 *p
++ = static_cast< typename
T::Char
>('1');
602 // maybe no nDigits if nDecPlaces < 0
603 *p
++ = static_cast< typename
T::Char
>('E');
607 *p
++ = static_cast< typename
T::Char
>('-');
610 *p
++ = static_cast< typename
T::Char
>('+');
612 *p
++ = static_cast< typename
T::Char
>(
613 nExp
/ 100 + static_cast< typename
T::Char
>('0') );
615 *p
++ = static_cast< typename
T::Char
>(
616 nExp
/ 10 + static_cast< typename
T::Char
>('0') );
617 *p
++ = static_cast< typename
T::Char
>(
618 nExp
% 10 + static_cast< typename
T::Char
>('0') );
621 if (pResultCapacity
== 0)
622 T::createString(pResult
, pBuf
, p
- pBuf
);
624 T::appendChars(pResult
, pResultCapacity
, &nResultOffset
, pBuf
,
627 if ( pBuf
!= &aBuf
[0] )
628 rtl_freeMemory(pBuf
);
633 void SAL_CALL
rtl_math_doubleToString(rtl_String
** pResult
,
634 sal_Int32
* pResultCapacity
,
635 sal_Int32 nResultOffset
, double fValue
,
636 rtl_math_StringFormat eFormat
,
637 sal_Int32 nDecPlaces
,
638 sal_Char cDecSeparator
,
639 sal_Int32
const * pGroups
,
640 sal_Char cGroupSeparator
,
641 sal_Bool bEraseTrailingDecZeros
)
644 doubleToString
< StringTraits
, StringTraits::String
>(
645 pResult
, pResultCapacity
, nResultOffset
, fValue
, eFormat
, nDecPlaces
,
646 cDecSeparator
, pGroups
, cGroupSeparator
, bEraseTrailingDecZeros
);
649 void SAL_CALL
rtl_math_doubleToUString(rtl_uString
** pResult
,
650 sal_Int32
* pResultCapacity
,
651 sal_Int32 nResultOffset
, double fValue
,
652 rtl_math_StringFormat eFormat
,
653 sal_Int32 nDecPlaces
,
654 sal_Unicode cDecSeparator
,
655 sal_Int32
const * pGroups
,
656 sal_Unicode cGroupSeparator
,
657 sal_Bool bEraseTrailingDecZeros
)
660 doubleToString
< UStringTraits
, UStringTraits::String
>(
661 pResult
, pResultCapacity
, nResultOffset
, fValue
, eFormat
, nDecPlaces
,
662 cDecSeparator
, pGroups
, cGroupSeparator
, bEraseTrailingDecZeros
);
667 // if nExp * 10 + nAdd would result in overflow
668 inline bool long10Overflow( long& nExp
, int nAdd
)
670 if ( nExp
> (LONG_MAX
/10)
671 || (nExp
== (LONG_MAX
/10) && nAdd
> (LONG_MAX
%10)) )
679 template< typename CharT
>
680 inline double stringToDouble(CharT
const * pBegin
, CharT
const * pEnd
,
681 CharT cDecSeparator
, CharT cGroupSeparator
,
682 rtl_math_ConversionStatus
* pStatus
,
683 CharT
const ** pParsedEnd
)
686 rtl_math_ConversionStatus eStatus
= rtl_math_ConversionStatus_Ok
;
688 CharT
const * p0
= pBegin
;
689 while (p0
!= pEnd
&& (*p0
== CharT(' ') || *p0
== CharT('\t')))
692 if (p0
!= pEnd
&& *p0
== CharT('-'))
700 if (p0
!= pEnd
&& *p0
== CharT('+'))
703 CharT
const * p
= p0
;
706 // #i112652# XMLSchema-2
709 if ((CharT('N') == p
[0]) && (CharT('a') == p
[1])
710 && (CharT('N') == p
[2]))
713 rtl::math::setNan( &fVal
);
716 else if ((CharT('I') == p
[0]) && (CharT('N') == p
[1])
717 && (CharT('F') == p
[2]))
721 eStatus
= rtl_math_ConversionStatus_OutOfRange
;
726 if (!bDone
) // do not recognize e.g. NaN1.23
728 // leading zeros and group separators may be safely ignored
729 while (p
!= pEnd
&& (*p
== CharT('0') || *p
== cGroupSeparator
))
732 long nValExp
= 0; // carry along exponent of mantissa
734 // integer part of mantissa
735 for (; p
!= pEnd
; ++p
)
738 if (rtl::isAsciiDigit(c
))
740 fVal
= fVal
* 10.0 + static_cast< double >( c
- CharT('0') );
743 else if (c
!= cGroupSeparator
)
747 // fraction part of mantissa
748 if (p
!= pEnd
&& *p
== cDecSeparator
)
753 while (p
!= pEnd
&& *p
== CharT('0'))
759 nValExp
= nFracExp
- 1; // no integer part => fraction exponent
760 // one decimal digit needs ld(10) ~= 3.32 bits
761 static const int nSigs
= (DBL_MANT_DIG
/ 3) + 1;
763 for (; p
!= pEnd
; ++p
)
766 if (!rtl::isAsciiDigit(c
))
769 { // further digits (more than nSigs) don't have any
771 fFrac
= fFrac
* 10.0 + static_cast<double>(c
- CharT('0'));
777 fVal
+= rtl::math::pow10Exp( fFrac
, nFracExp
);
778 else if ( nValExp
< 0 )
779 nValExp
= 0; // no digit other than 0 after decimal point
783 --nValExp
; // started with offset +1 at the first mantissa digit
786 if (p
!= p0
&& p
!= pEnd
&& (*p
== CharT('E') || *p
== CharT('e')))
788 CharT
const * const pExponent
= p
;
791 if (p
!= pEnd
&& *p
== CharT('-'))
799 if (p
!= pEnd
&& *p
== CharT('+'))
802 CharT
const * const pFirstExpDigit
= p
;
804 { // no matter what follows, zero stays zero, but carry on the
806 while (p
!= pEnd
&& rtl::isAsciiDigit(*p
))
808 if (p
== pFirstExpDigit
)
809 { // no digits in exponent, reset end of scan
815 bool bOverflow
= false;
817 for (; p
!= pEnd
; ++p
)
820 if (!rtl::isAsciiDigit(c
))
822 int i
= c
- CharT('0');
823 if ( long10Overflow( nExp
, i
) )
826 nExp
= nExp
* 10 + i
;
832 long nAllExp
= ( bOverflow
? 0 : nExp
+ nValExp
);
833 if ( nAllExp
> DBL_MAX_10_EXP
|| (bOverflow
&& !bExpSign
) )
836 eStatus
= rtl_math_ConversionStatus_OutOfRange
;
838 else if ((nAllExp
< DBL_MIN_10_EXP
) ||
839 (bOverflow
&& bExpSign
) )
842 eStatus
= rtl_math_ConversionStatus_OutOfRange
;
844 else if ( nExp
> DBL_MAX_10_EXP
|| nExp
< DBL_MIN_10_EXP
)
845 { // compensate exponents
846 fVal
= rtl::math::pow10Exp( fVal
, -nValExp
);
847 fVal
= rtl::math::pow10Exp( fVal
, nAllExp
);
850 fVal
= rtl::math::pow10Exp( fVal
, nExp
); // normal
852 else if (p
== pFirstExpDigit
)
853 { // no digits in exponent, reset end of scan
858 else if (p
- p0
== 2 && p
!= pEnd
&& p
[0] == CharT('#')
859 && p
[-1] == cDecSeparator
&& p
[-2] == CharT('1'))
861 if (pEnd
- p
>= 4 && p
[1] == CharT('I') && p
[2] == CharT('N')
862 && p
[3] == CharT('F'))
864 // "1.#INF", "+1.#INF", "-1.#INF"
867 eStatus
= rtl_math_ConversionStatus_OutOfRange
;
868 // Eat any further digits:
869 while (p
!= pEnd
&& rtl::isAsciiDigit(*p
))
872 else if (pEnd
- p
>= 4 && p
[1] == CharT('N') && p
[2] == CharT('A')
873 && p
[3] == CharT('N'))
875 // "1.#NAN", "+1.#NAN", "-1.#NAN"
877 rtl::math::setNan( &fVal
);
885 m
.md
.w32_parts
.msw
|= 0x80000000; // create negative NaN
887 bSign
= false; // don't negate again
889 // Eat any further digits:
890 while (p
!= pEnd
&& rtl::isAsciiDigit(*p
))
896 // overflow also if more than DBL_MAX_10_EXP digits without decimal
897 // separator, or 0. and more than DBL_MIN_10_EXP digits, ...
898 bool bHuge
= fVal
== HUGE_VAL
; // g++ 3.0.1 requires it this way...
900 eStatus
= rtl_math_ConversionStatus_OutOfRange
;
908 *pParsedEnd
= p
== p0
? pBegin
: p
;
915 double SAL_CALL
rtl_math_stringToDouble(sal_Char
const * pBegin
,
916 sal_Char
const * pEnd
,
917 sal_Char cDecSeparator
,
918 sal_Char cGroupSeparator
,
919 rtl_math_ConversionStatus
* pStatus
,
920 sal_Char
const ** pParsedEnd
)
923 return stringToDouble(pBegin
, pEnd
, cDecSeparator
, cGroupSeparator
, pStatus
,
927 double SAL_CALL
rtl_math_uStringToDouble(sal_Unicode
const * pBegin
,
928 sal_Unicode
const * pEnd
,
929 sal_Unicode cDecSeparator
,
930 sal_Unicode cGroupSeparator
,
931 rtl_math_ConversionStatus
* pStatus
,
932 sal_Unicode
const ** pParsedEnd
)
935 return stringToDouble(pBegin
, pEnd
, cDecSeparator
, cGroupSeparator
, pStatus
,
939 double SAL_CALL
rtl_math_round(double fValue
, int nDecPlaces
,
940 enum rtl_math_RoundingMode eMode
)
943 OSL_ASSERT(nDecPlaces
>= -20 && nDecPlaces
<= 20);
949 bool bSign
= rtl::math::isSignBitSet( fValue
);
954 if ( nDecPlaces
!= 0 )
956 // max 20 decimals, we don't have unlimited precision
957 // #38810# and no overflow on fValue*=fFac
958 if ( nDecPlaces
< -20 || 20 < nDecPlaces
|| fValue
> (DBL_MAX
/ 1e20
) )
959 return bSign
? -fValue
: fValue
;
961 fFac
= getN10Exp( nDecPlaces
);
964 //else //! uninitialized fFac, not needed
968 case rtl_math_RoundingMode_Corrected
:
970 int nExp
; // exponent for correction
972 nExp
= static_cast<int>( floor( log10( fValue
) ) );
975 int nIndex
= 15 - nExp
;
978 else if ( nIndex
<= 1 )
980 fValue
= floor( fValue
+ 0.5 + nKorrVal
[nIndex
] );
983 case rtl_math_RoundingMode_Down
:
984 fValue
= rtl::math::approxFloor( fValue
);
986 case rtl_math_RoundingMode_Up
:
987 fValue
= rtl::math::approxCeil( fValue
);
989 case rtl_math_RoundingMode_Floor
:
990 fValue
= bSign
? rtl::math::approxCeil( fValue
)
991 : rtl::math::approxFloor( fValue
);
993 case rtl_math_RoundingMode_Ceiling
:
994 fValue
= bSign
? rtl::math::approxFloor( fValue
)
995 : rtl::math::approxCeil( fValue
);
997 case rtl_math_RoundingMode_HalfDown
:
999 double f
= floor( fValue
);
1000 fValue
= ((fValue
- f
) <= 0.5) ? f
: ceil( fValue
);
1003 case rtl_math_RoundingMode_HalfUp
:
1005 double f
= floor( fValue
);
1006 fValue
= ((fValue
- f
) < 0.5) ? f
: ceil( fValue
);
1009 case rtl_math_RoundingMode_HalfEven
:
1010 #if defined FLT_ROUNDS
1012 Use fast version. FLT_ROUNDS may be defined to a function by some compilers!
1014 DBL_EPSILON is the smallest fractional number which can be represented,
1015 its reciprocal is therefore the smallest number that cannot have a
1016 fractional part. Once you add this reciprocal to `x', its fractional part
1017 is stripped off. Simply subtracting the reciprocal back out returns `x'
1018 without its fractional component.
1019 Simple, clever, and elegant - thanks to Ross Cottrell, the original author,
1020 who placed it into public domain.
1022 volatile: prevent compiler from being too smart
1024 if ( FLT_ROUNDS
== 1 )
1026 volatile double x
= fValue
+ 1.0 / DBL_EPSILON
;
1027 fValue
= x
- 1.0 / DBL_EPSILON
;
1030 #endif // FLT_ROUNDS
1032 double f
= floor( fValue
);
1033 if ( (fValue
- f
) != 0.5 )
1034 fValue
= floor( fValue
+ 0.5 );
1038 fValue
= (g
== floor( g
)) ? f
: (f
+ 1.0);
1047 if ( nDecPlaces
!= 0 )
1050 return bSign
? -fValue
: fValue
;
1053 double SAL_CALL
rtl_math_pow10Exp(double fValue
, int nExp
) SAL_THROW_EXTERN_C()
1055 return fValue
* getN10Exp( nExp
);
1058 double SAL_CALL
rtl_math_approxValue( double fValue
) SAL_THROW_EXTERN_C()
1060 if (fValue
== 0.0 || fValue
== HUGE_VAL
|| !::rtl::math::isFinite( fValue
))
1061 // We don't handle these conditions. Bail out.
1064 double fOrigValue
= fValue
;
1066 bool bSign
= ::rtl::math::isSignBitSet( fValue
);
1070 int nExp
= static_cast<int>( floor( log10( fValue
)));
1072 double fExpValue
= getN10Exp( nExp
);
1074 fValue
*= fExpValue
;
1075 // If the original value was near DBL_MIN we got an overflow. Restore and
1077 if (!rtl::math::isFinite( fValue
))
1079 fValue
= rtl_math_round( fValue
, 0, rtl_math_RoundingMode_Corrected
);
1080 fValue
/= fExpValue
;
1081 // If the original value was near DBL_MAX we got an overflow. Restore and
1083 if (!rtl::math::isFinite( fValue
))
1086 return bSign
? -fValue
: fValue
;
1089 double SAL_CALL
rtl_math_expm1( double fValue
) SAL_THROW_EXTERN_C()
1091 double fe
= exp( fValue
);
1096 return (fe
-1.0) * fValue
/ log(fe
);
1099 double SAL_CALL
rtl_math_log1p( double fValue
) SAL_THROW_EXTERN_C()
1101 // Use volatile because a compiler may be too smart "optimizing" the
1102 // condition such that in certain cases the else path was called even if
1103 // (fp==1.0) was true, where the term (fp-1.0) then resulted in 0.0 and
1104 // hence the entire expression resulted in NaN.
1105 // Happened with g++ 3.4.1 and an input value of 9.87E-18
1106 volatile double fp
= 1.0 + fValue
;
1110 return log(fp
) * fValue
/ (fp
-1.0);
1113 double SAL_CALL
rtl_math_atanh( double fValue
) SAL_THROW_EXTERN_C()
1115 return 0.5 * rtl_math_log1p( 2.0 * fValue
/ (1.0-fValue
) );
1118 /** Parent error function (erf) that calls different algorithms based on the
1119 value of x. It takes care of cases where x is negative as erf is an odd
1120 function i.e. erf(-x) = -erf(x).
1122 Kramer, W., and Blomquist, F., 2000, Algorithms with Guaranteed Error Bounds
1123 for the Error Function and the Complementary Error Function
1125 http://www.math.uni-wuppertal.de/wrswt/literatur_en.html
1127 @author Kohei Yoshida <kohei@openoffice.org>
1131 double SAL_CALL
rtl_math_erf( double x
) SAL_THROW_EXTERN_C()
1136 bool bNegative
= false;
1145 fErf
= (double) (x
*1.1283791670955125738961589031215452L);
1146 else if ( x
< 0.65 )
1147 lcl_Erf0065( x
, fErf
);
1149 fErf
= 1.0 - rtl_math_erfc( x
);
1157 /** Parent complementary error function (erfc) that calls different algorithms
1158 based on the value of x. It takes care of cases where x is negative as erfc
1159 satisfies relationship erfc(-x) = 2 - erfc(x). See the comment for Erf(x)
1160 for the source publication.
1162 @author Kohei Yoshida <kohei@openoffice.org>
1164 @see #i55735#, moved from module scaddins (#i97091#)
1167 double SAL_CALL
rtl_math_erfc( double x
) SAL_THROW_EXTERN_C()
1172 bool bNegative
= false;
1183 lcl_Erfc0600( x
, fErfc
);
1185 lcl_Erfc2654( x
, fErfc
);
1188 fErfc
= 1.0 - rtl_math_erf( x
);
1191 fErfc
= 2.0 - fErfc
;
1196 /** improved accuracy of asinh for |x| large and for x near zero
1199 double SAL_CALL
rtl_math_asinh( double fX
) SAL_THROW_EXTERN_C()
1212 return fSign
* rtl_math_log1p( fX
+ fX
*fX
/ (1.0 + sqrt( 1.0 + fX
*fX
)));
1213 else if ( fX
< 1.25e7
)
1214 return fSign
* log( fX
+ sqrt( 1.0 + fX
*fX
));
1216 return fSign
* log( 2.0*fX
);
1220 /** improved accuracy of acosh for x large and for x near 1
1223 double SAL_CALL
rtl_math_acosh( double fX
) SAL_THROW_EXTERN_C()
1225 volatile double fZ
= fX
- 1.0;
1229 ::rtl::math::setNan( &fResult
);
1232 else if ( fX
== 1.0 )
1234 else if ( fX
< 1.1 )
1235 return rtl_math_log1p( fZ
+ sqrt( fZ
*fZ
+ 2.0*fZ
));
1236 else if ( fX
< 1.25e7
)
1237 return log( fX
+ sqrt( fX
*fX
- 1.0));
1239 return log( 2.0*fX
);
1242 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */