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