Bug Summary

File:sal/rtl/source/math.cxx
Location:line 622, column 13
Description:Dereference of null pointer

Annotated Source Code

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