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