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