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