?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
tdf#150621 Changed Korean word counting to use words
[LibreOffice.git] / sal / rtl / math.cxx
blob766bde5d92b95302c6bf8ca97f80c5d66fe1cc56
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/character.hxx>
24 #include <rtl/math.hxx>
25
26 #include <algorithm>
27 #include <bit>
28 #include <cassert>
29 #include <cfenv>
30 #include <cmath>
31 #include <float.h>
32 #include <limits>
33 #include <memory>
34 #include <stdlib.h>
35
36 #include "strtmpl.hxx"
37
38 #include <dtoa.h>
39
40 constexpr int minExp = -323, maxExp = 308;
41 constexpr double n10s[] = {
42 1e-323, 1e-322, 1e-321, 1e-320, 1e-319, 1e-318, 1e-317, 1e-316, 1e-315, 1e-314, 1e-313, 1e-312,
43 1e-311, 1e-310, 1e-309, 1e-308, 1e-307, 1e-306, 1e-305, 1e-304, 1e-303, 1e-302, 1e-301, 1e-300,
44 1e-299, 1e-298, 1e-297, 1e-296, 1e-295, 1e-294, 1e-293, 1e-292, 1e-291, 1e-290, 1e-289, 1e-288,
45 1e-287, 1e-286, 1e-285, 1e-284, 1e-283, 1e-282, 1e-281, 1e-280, 1e-279, 1e-278, 1e-277, 1e-276,
46 1e-275, 1e-274, 1e-273, 1e-272, 1e-271, 1e-270, 1e-269, 1e-268, 1e-267, 1e-266, 1e-265, 1e-264,
47 1e-263, 1e-262, 1e-261, 1e-260, 1e-259, 1e-258, 1e-257, 1e-256, 1e-255, 1e-254, 1e-253, 1e-252,
48 1e-251, 1e-250, 1e-249, 1e-248, 1e-247, 1e-246, 1e-245, 1e-244, 1e-243, 1e-242, 1e-241, 1e-240,
49 1e-239, 1e-238, 1e-237, 1e-236, 1e-235, 1e-234, 1e-233, 1e-232, 1e-231, 1e-230, 1e-229, 1e-228,
50 1e-227, 1e-226, 1e-225, 1e-224, 1e-223, 1e-222, 1e-221, 1e-220, 1e-219, 1e-218, 1e-217, 1e-216,
51 1e-215, 1e-214, 1e-213, 1e-212, 1e-211, 1e-210, 1e-209, 1e-208, 1e-207, 1e-206, 1e-205, 1e-204,
52 1e-203, 1e-202, 1e-201, 1e-200, 1e-199, 1e-198, 1e-197, 1e-196, 1e-195, 1e-194, 1e-193, 1e-192,
53 1e-191, 1e-190, 1e-189, 1e-188, 1e-187, 1e-186, 1e-185, 1e-184, 1e-183, 1e-182, 1e-181, 1e-180,
54 1e-179, 1e-178, 1e-177, 1e-176, 1e-175, 1e-174, 1e-173, 1e-172, 1e-171, 1e-170, 1e-169, 1e-168,
55 1e-167, 1e-166, 1e-165, 1e-164, 1e-163, 1e-162, 1e-161, 1e-160, 1e-159, 1e-158, 1e-157, 1e-156,
56 1e-155, 1e-154, 1e-153, 1e-152, 1e-151, 1e-150, 1e-149, 1e-148, 1e-147, 1e-146, 1e-145, 1e-144,
57 1e-143, 1e-142, 1e-141, 1e-140, 1e-139, 1e-138, 1e-137, 1e-136, 1e-135, 1e-134, 1e-133, 1e-132,
58 1e-131, 1e-130, 1e-129, 1e-128, 1e-127, 1e-126, 1e-125, 1e-124, 1e-123, 1e-122, 1e-121, 1e-120,
59 1e-119, 1e-118, 1e-117, 1e-116, 1e-115, 1e-114, 1e-113, 1e-112, 1e-111, 1e-110, 1e-109, 1e-108,
60 1e-107, 1e-106, 1e-105, 1e-104, 1e-103, 1e-102, 1e-101, 1e-100, 1e-99, 1e-98, 1e-97, 1e-96,
61 1e-95, 1e-94, 1e-93, 1e-92, 1e-91, 1e-90, 1e-89, 1e-88, 1e-87, 1e-86, 1e-85, 1e-84,
62 1e-83, 1e-82, 1e-81, 1e-80, 1e-79, 1e-78, 1e-77, 1e-76, 1e-75, 1e-74, 1e-73, 1e-72,
63 1e-71, 1e-70, 1e-69, 1e-68, 1e-67, 1e-66, 1e-65, 1e-64, 1e-63, 1e-62, 1e-61, 1e-60,
64 1e-59, 1e-58, 1e-57, 1e-56, 1e-55, 1e-54, 1e-53, 1e-52, 1e-51, 1e-50, 1e-49, 1e-48,
65 1e-47, 1e-46, 1e-45, 1e-44, 1e-43, 1e-42, 1e-41, 1e-40, 1e-39, 1e-38, 1e-37, 1e-36,
66 1e-35, 1e-34, 1e-33, 1e-32, 1e-31, 1e-30, 1e-29, 1e-28, 1e-27, 1e-26, 1e-25, 1e-24,
67 1e-23, 1e-22, 1e-21, 1e-20, 1e-19, 1e-18, 1e-17, 1e-16, 1e-15, 1e-14, 1e-13, 1e-12,
68 1e-11, 1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0,
69 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12,
70 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22, 1e23, 1e24,
71 1e25, 1e26, 1e27, 1e28, 1e29, 1e30, 1e31, 1e32, 1e33, 1e34, 1e35, 1e36,
72 1e37, 1e38, 1e39, 1e40, 1e41, 1e42, 1e43, 1e44, 1e45, 1e46, 1e47, 1e48,
73 1e49, 1e50, 1e51, 1e52, 1e53, 1e54, 1e55, 1e56, 1e57, 1e58, 1e59, 1e60,
74 1e61, 1e62, 1e63, 1e64, 1e65, 1e66, 1e67, 1e68, 1e69, 1e70, 1e71, 1e72,
75 1e73, 1e74, 1e75, 1e76, 1e77, 1e78, 1e79, 1e80, 1e81, 1e82, 1e83, 1e84,
76 1e85, 1e86, 1e87, 1e88, 1e89, 1e90, 1e91, 1e92, 1e93, 1e94, 1e95, 1e96,
77 1e97, 1e98, 1e99, 1e100, 1e101, 1e102, 1e103, 1e104, 1e105, 1e106, 1e107, 1e108,
78 1e109, 1e110, 1e111, 1e112, 1e113, 1e114, 1e115, 1e116, 1e117, 1e118, 1e119, 1e120,
79 1e121, 1e122, 1e123, 1e124, 1e125, 1e126, 1e127, 1e128, 1e129, 1e130, 1e131, 1e132,
80 1e133, 1e134, 1e135, 1e136, 1e137, 1e138, 1e139, 1e140, 1e141, 1e142, 1e143, 1e144,
81 1e145, 1e146, 1e147, 1e148, 1e149, 1e150, 1e151, 1e152, 1e153, 1e154, 1e155, 1e156,
82 1e157, 1e158, 1e159, 1e160, 1e161, 1e162, 1e163, 1e164, 1e165, 1e166, 1e167, 1e168,
83 1e169, 1e170, 1e171, 1e172, 1e173, 1e174, 1e175, 1e176, 1e177, 1e178, 1e179, 1e180,
84 1e181, 1e182, 1e183, 1e184, 1e185, 1e186, 1e187, 1e188, 1e189, 1e190, 1e191, 1e192,
85 1e193, 1e194, 1e195, 1e196, 1e197, 1e198, 1e199, 1e200, 1e201, 1e202, 1e203, 1e204,
86 1e205, 1e206, 1e207, 1e208, 1e209, 1e210, 1e211, 1e212, 1e213, 1e214, 1e215, 1e216,
87 1e217, 1e218, 1e219, 1e220, 1e221, 1e222, 1e223, 1e224, 1e225, 1e226, 1e227, 1e228,
88 1e229, 1e230, 1e231, 1e232, 1e233, 1e234, 1e235, 1e236, 1e237, 1e238, 1e239, 1e240,
89 1e241, 1e242, 1e243, 1e244, 1e245, 1e246, 1e247, 1e248, 1e249, 1e250, 1e251, 1e252,
90 1e253, 1e254, 1e255, 1e256, 1e257, 1e258, 1e259, 1e260, 1e261, 1e262, 1e263, 1e264,
91 1e265, 1e266, 1e267, 1e268, 1e269, 1e270, 1e271, 1e272, 1e273, 1e274, 1e275, 1e276,
92 1e277, 1e278, 1e279, 1e280, 1e281, 1e282, 1e283, 1e284, 1e285, 1e286, 1e287, 1e288,
93 1e289, 1e290, 1e291, 1e292, 1e293, 1e294, 1e295, 1e296, 1e297, 1e298, 1e299, 1e300,
94 1e301, 1e302, 1e303, 1e304, 1e305, 1e306, 1e307, 1e308,
95 };
96 static_assert(SAL_N_ELEMENTS(n10s) == maxExp - minExp + 1);
97
98 // return pow(10.0,nExp) optimized for exponents in the interval [-323,308] (i.e., incl. denormals)
99 static double getN10Exp(int nExp)
100 {
101 if (nExp < minExp || nExp > maxExp)
102 return pow(10.0, static_cast<double>(nExp)); // will return 0 or INF with IEEE 754
103 return n10s[nExp - minExp];
104 }
105
106 namespace
107 {
108 /** If value (passed as absolute value) is an integer representable as double,
109 which we handle explicitly at some places.
110 */
111 bool isRepresentableInteger(double fAbsValue)
112 {
113 static_assert(std::numeric_limits<double>::is_iec559
114 && std::numeric_limits<double>::digits == 53);
115 assert(fAbsValue >= 0.0);
116 if (fAbsValue >= 0x1p53)
117 return false;
118 sal_Int64 nInt = static_cast<sal_Int64>(fAbsValue);
119 return nInt == fAbsValue;
120 }
121
122 /** Returns number of binary bits for fractional part of the number
123 Expects a proper non-negative double value, not +-INF, not NAN
124 */
125 int getBitsInFracPart(double fAbsValue)
126 {
127 assert(std::isfinite(fAbsValue) && fAbsValue >= 0.0);
128 if (fAbsValue == 0.0)
129 return 0;
130 auto& rValParts = reinterpret_cast<const sal_math_Double*>(&fAbsValue)->parts;
131 int nExponent = rValParts.exponent - 1023;
132 if (nExponent >= 52)
133 return 0; // All bits in fraction are in integer part of the number
134 int nLeastSignificant = rValParts.fraction
135 ? std::countr_zero(rValParts.fraction) + 1
136 : 53; // the implied leading 1 is the least significant
137 int nFracSignificant = 53 - nLeastSignificant;
138 int nBitsInFracPart = nFracSignificant - nExponent;
139
140 return std::max(nBitsInFracPart, 0);
141 }
142 }
143
144 void SAL_CALL rtl_math_doubleToString(rtl_String** pResult, sal_Int32* pResultCapacity,
145 sal_Int32 nResultOffset, double fValue,
146 rtl_math_StringFormat eFormat, sal_Int32 nDecPlaces,
147 char cDecSeparator, sal_Int32 const* pGroups,
148 char cGroupSeparator, sal_Bool bEraseTrailingDecZeros)
149 SAL_THROW_EXTERN_C()
150 {
151 rtl::str::doubleToString(pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
152 cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
153 }
154
155 void SAL_CALL rtl_math_doubleToUString(rtl_uString** pResult, sal_Int32* pResultCapacity,
156 sal_Int32 nResultOffset, double fValue,
157 rtl_math_StringFormat eFormat, sal_Int32 nDecPlaces,
158 sal_Unicode cDecSeparator, sal_Int32 const* pGroups,
159 sal_Unicode cGroupSeparator, sal_Bool bEraseTrailingDecZeros)
160 SAL_THROW_EXTERN_C()
161 {
162 rtl::str::doubleToString(pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
163 cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
164 }
165
166 namespace
167 {
168 template <typename CharT>
169 double stringToDouble(CharT const* pBegin, CharT const* pEnd, CharT cDecSeparator,
170 CharT cGroupSeparator, rtl_math_ConversionStatus* pStatus,
171 CharT const** pParsedEnd)
172 {
173 double fVal = 0.0;
174 rtl_math_ConversionStatus eStatus = rtl_math_ConversionStatus_Ok;
175
176 CharT const* p0 = pBegin;
177 while (p0 != pEnd && (*p0 == ' ' || *p0 == '\t'))
178 {
179 ++p0;
180 }
181
182 bool bSign;
183 bool explicitSign = false;
184 if (p0 != pEnd && *p0 == '-')
185 {
186 bSign = true;
187 explicitSign = true;
188 ++p0;
189 }
190 else
191 {
192 bSign = false;
193 if (p0 != pEnd && *p0 == '+')
194 {
195 explicitSign = true;
196 ++p0;
197 }
198 }
199
200 CharT const* p = p0;
201 bool bDone = false;
202
203 // #i112652# XMLSchema-2
204 if ((pEnd - p) >= 3)
205 {
206 if (!explicitSign && ('N' == p[0]) && ('a' == p[1]) && ('N' == p[2]))
207 {
208 p += 3;
209 fVal = std::numeric_limits<double>::quiet_NaN();
210 bDone = true;
211 }
212 else if (('I' == p[0]) && ('N' == p[1]) && ('F' == p[2]))
213 {
214 p += 3;
215 fVal = HUGE_VAL;
216 eStatus = rtl_math_ConversionStatus_OutOfRange;
217 bDone = true;
218 }
219 }
220
221 if (!bDone) // do not recognize e.g. NaN1.23
222 {
223 std::unique_ptr<char[]> bufInHeap;
224 std::unique_ptr<const CharT* []> bufInHeapMap;
225 constexpr int bufOnStackSize = 256;
226 char bufOnStack[bufOnStackSize];
227 const CharT* bufOnStackMap[bufOnStackSize];
228 char* buf = bufOnStack;
229 const CharT** bufmap = bufOnStackMap;
230 int bufpos = 0;
231 const size_t bufsize = pEnd - p + (bSign ? 2 : 1);
232 if (bufsize > bufOnStackSize)
233 {
234 bufInHeap = std::make_unique<char[]>(bufsize);
235 bufInHeapMap = std::make_unique<const CharT* []>(bufsize);
236 buf = bufInHeap.get();
237 bufmap = bufInHeapMap.get();
238 }
239
240 if (bSign)
241 {
242 buf[0] = '-';
243 bufmap[0] = p; // yes, this may be the same pointer as for the next mapping
244 bufpos = 1;
245 }
246 // Put first zero to buffer for strings like "-0"
247 if (p != pEnd && *p == '0')
248 {
249 buf[bufpos] = '0';
250 bufmap[bufpos] = p;
251 ++bufpos;
252 ++p;
253 }
254 // Leading zeros and group separators between digits may be safely
255 // ignored. p0 < p implies that there was a leading 0 already,
256 // consecutive group separators may not happen as *(p+1) is checked for
257 // digit.
258 while (p != pEnd
259 && (*p == '0'
260 || (*p == cGroupSeparator && p0 < p && p + 1 < pEnd
261 && rtl::isAsciiDigit(*(p + 1)))))
262 {
263 ++p;
264 }
265
266 // integer part of mantissa
267 for (; p != pEnd; ++p)
268 {
269 CharT c = *p;
270 if (rtl::isAsciiDigit(c))
271 {
272 buf[bufpos] = static_cast<char>(c);
273 bufmap[bufpos] = p;
274 ++bufpos;
275 }
276 else if (c != cGroupSeparator)
277 {
278 break;
279 }
280 else if (p == p0 || (p + 1 == pEnd) || !rtl::isAsciiDigit(*(p + 1)))
281 {
282 // A leading or trailing (not followed by a digit) group
283 // separator character is not a group separator.
284 break;
285 }
286 }
287
288 // fraction part of mantissa
289 if (p != pEnd && *p == cDecSeparator)
290 {
291 buf[bufpos] = '.';
292 bufmap[bufpos] = p;
293 ++bufpos;
294 ++p;
295
296 for (; p != pEnd; ++p)
297 {
298 CharT c = *p;
299 if (!rtl::isAsciiDigit(c))
300 {
301 break;
302 }
303 buf[bufpos] = static_cast<char>(c);
304 bufmap[bufpos] = p;
305 ++bufpos;
306 }
307 }
308
309 // Exponent
310 if (p != p0 && p != pEnd && (*p == 'E' || *p == 'e'))
311 {
312 buf[bufpos] = 'E';
313 bufmap[bufpos] = p;
314 ++bufpos;
315 ++p;
316 if (p != pEnd && *p == '-')
317 {
318 buf[bufpos] = '-';
319 bufmap[bufpos] = p;
320 ++bufpos;
321 ++p;
322 }
323 else if (p != pEnd && *p == '+')
324 ++p;
325
326 for (; p != pEnd; ++p)
327 {
328 CharT c = *p;
329 if (!rtl::isAsciiDigit(c))
330 break;
331
332 buf[bufpos] = static_cast<char>(c);
333 bufmap[bufpos] = p;
334 ++bufpos;
335 }
336 }
337 else if (p - p0 == 2 && p != pEnd && p[0] == '#' && p[-1] == cDecSeparator && p[-2] == '1')
338 {
339 if (pEnd - p >= 4 && p[1] == 'I' && p[2] == 'N' && p[3] == 'F')
340 {
341 // "1.#INF", "+1.#INF", "-1.#INF"
342 p += 4;
343 fVal = HUGE_VAL;
344 eStatus = rtl_math_ConversionStatus_OutOfRange;
345 // Eat any further digits:
346 while (p != pEnd && rtl::isAsciiDigit(*p))
347 ++p;
348 bDone = true;
349 }
350 else if (pEnd - p >= 4 && p[1] == 'N' && p[2] == 'A' && p[3] == 'N')
351 {
352 // "1.#NAN", "+1.#NAN", "-1.#NAN"
353 p += 4;
354 fVal = std::copysign(std::numeric_limits<double>::quiet_NaN(), bSign ? -1.0 : 1.0);
355 bSign = false; // don't negate again
356
357 // Eat any further digits:
358 while (p != pEnd && rtl::isAsciiDigit(*p))
359 {
360 ++p;
361 }
362 bDone = true;
363 }
364 }
365
366 if (!bDone)
367 {
368 buf[bufpos] = '\0';
369 bufmap[bufpos] = p;
370 char* pCharParseEnd;
371 errno = 0;
372 fVal = strtod_nolocale(buf, &pCharParseEnd);
373 if (errno == ERANGE)
374 {
375 // Check for the dreaded rounded to 15 digits max value
376 // 1.79769313486232e+308 for 1.7976931348623157e+308 we wrote
377 // everywhere, accept with or without plus sign in exponent.
378 const char* b = buf;
379 if (b[0] == '-')
380 ++b;
381 if (((pCharParseEnd - b == 21) || (pCharParseEnd - b == 20))
382 && !strncmp(b, "1.79769313486232", 16) && (b[16] == 'e' || b[16] == 'E')
383 && (((pCharParseEnd - b == 21) && !strncmp(b + 17, "+308", 4))
384 || ((pCharParseEnd - b == 20) && !strncmp(b + 17, "308", 3))))
385 {
386 fVal = (buf < b) ? -DBL_MAX : DBL_MAX;
387 }
388 else
389 {
390 eStatus = rtl_math_ConversionStatus_OutOfRange;
391 }
392 }
393 p = bufmap[pCharParseEnd - buf];
394 bSign = false;
395 }
396 }
397
398 // overflow also if more than DBL_MAX_10_EXP digits without decimal
399 // separator, or 0. and more than DBL_MIN_10_EXP digits, ...
400 if (std::isinf(fVal))
401 eStatus = rtl_math_ConversionStatus_OutOfRange;
402
403 if (bSign)
404 fVal = -fVal;
405
406 if (pStatus)
407 *pStatus = eStatus;
408
409 if (pParsedEnd)
410 *pParsedEnd = p == p0 ? pBegin : p;
411
412 return fVal;
413 }
414 }
415
416 double SAL_CALL rtl_math_stringToDouble(char const* pBegin, char const* pEnd, char cDecSeparator,
417 char cGroupSeparator, rtl_math_ConversionStatus* pStatus,
418 char const** pParsedEnd) SAL_THROW_EXTERN_C()
419 {
420 return stringToDouble(reinterpret_cast<unsigned char const*>(pBegin),
421 reinterpret_cast<unsigned char const*>(pEnd),
422 static_cast<unsigned char>(cDecSeparator),
423 static_cast<unsigned char>(cGroupSeparator), pStatus,
424 reinterpret_cast<unsigned char const**>(pParsedEnd));
425 }
426
427 double SAL_CALL rtl_math_uStringToDouble(sal_Unicode const* pBegin, sal_Unicode const* pEnd,
428 sal_Unicode cDecSeparator, sal_Unicode cGroupSeparator,
429 rtl_math_ConversionStatus* pStatus,
430 sal_Unicode const** pParsedEnd) SAL_THROW_EXTERN_C()
431 {
432 return stringToDouble(pBegin, pEnd, cDecSeparator, cGroupSeparator, pStatus, pParsedEnd);
433 }
434
435 double SAL_CALL rtl_math_round(double fValue, int nDecPlaces, enum rtl_math_RoundingMode eMode)
436 SAL_THROW_EXTERN_C()
437 {
438 if (!std::isfinite(fValue))
439 return fValue;
440
441 if (fValue == 0.0)
442 return fValue;
443
444 const double fOrigValue = fValue;
445
446 // sign adjustment
447 bool bSign = std::signbit(fValue);
448 if (bSign)
449 fValue = -fValue;
450
451 // Rounding to decimals between integer distance precision (gaps) does not
452 // make sense, do not even try to multiply/divide and introduce inaccuracy.
453 // For same reasons, do not attempt to round integers to decimals.
454 if (nDecPlaces >= 0 && (fValue >= 0x1p52 || isRepresentableInteger(fValue)))
455 return fOrigValue;
456
457 double fFac = 0;
458 if (nDecPlaces != 0)
459 {
460 if (nDecPlaces > 0)
461 {
462 // Determine how many decimals are representable in the precision.
463 // Anything greater 2^52 and 0.0 was already ruled out above.
464 // Theoretically 0.5, 0.25, 0.125, 0.0625, 0.03125, ...
465 const sal_math_Double* pd = reinterpret_cast<const sal_math_Double*>(&fValue);
466 const sal_Int32 nDec = 52 - (pd->parts.exponent - 1023);
467
468 if (nDec <= 0)
469 {
470 assert(!"Shouldn't this had been caught already as large number?");
471 return fOrigValue;
472 }
473
474 if (nDec < nDecPlaces)
475 nDecPlaces = nDec;
476 }
477
478 // Avoid 1e-5 (1.0000000000000001e-05) and such inaccurate fractional
479 // factors that later when dividing back spoil things. For negative
480 // decimals divide first with the inverse, then multiply the rounded
481 // value back.
482 fFac = getN10Exp(abs(nDecPlaces));
483
484 if (fFac == 0.0 || (nDecPlaces < 0 && !std::isfinite(fFac)))
485 // Underflow, rounding to that many integer positions would be 0.
486 return 0.0;
487
488 if (!std::isfinite(fFac))
489 // Overflow with very small values and high number of decimals.
490 return fOrigValue;
491
492 if (nDecPlaces < 0)
493 fValue /= fFac;
494 else
495 fValue *= fFac;
496
497 if (!std::isfinite(fValue))
498 return fOrigValue;
499 }
500
501 // Round only if not already in distance precision gaps of integers, where
502 // for [2^52,2^53) adding 0.5 would even yield the next representable
503 // integer.
504 if (fValue < 0x1p52)
505 {
506 switch (eMode)
507 {
508 case rtl_math_RoundingMode_Corrected:
509 fValue = rtl::math::approxFloor(fValue + 0.5);
510 break;
511 case rtl_math_RoundingMode_Down:
512 fValue = rtl::math::approxFloor(fValue);
513 break;
514 case rtl_math_RoundingMode_Up:
515 fValue = rtl::math::approxCeil(fValue);
516 break;
517 case rtl_math_RoundingMode_Floor:
518 fValue = bSign ? rtl::math::approxCeil(fValue) : rtl::math::approxFloor(fValue);
519 break;
520 case rtl_math_RoundingMode_Ceiling:
521 fValue = bSign ? rtl::math::approxFloor(fValue) : rtl::math::approxCeil(fValue);
522 break;
523 case rtl_math_RoundingMode_HalfDown:
524 {
525 double f = floor(fValue);
526 fValue = ((fValue - f) <= 0.5) ? f : ceil(fValue);
527 }
528 break;
529 case rtl_math_RoundingMode_HalfUp:
530 {
531 double f = floor(fValue);
532 fValue = ((fValue - f) < 0.5) ? f : ceil(fValue);
533 }
534 break;
535 case rtl_math_RoundingMode_HalfEven:
536 if (const int oldMode = std::fegetround(); std::fesetround(FE_TONEAREST) == 0)
537 {
538 fValue = std::nearbyint(fValue);
539 std::fesetround(oldMode);
540 }
541 else
542 {
543 double f = floor(fValue);
544 if ((fValue - f) != 0.5)
545 {
546 fValue = floor(fValue + 0.5);
547 }
548 else
549 {
550 double g = f / 2.0;
551 fValue = (g == floor(g)) ? f : (f + 1.0);
552 }
553 }
554 break;
555 default:
556 OSL_ASSERT(false);
557 break;
558 }
559 }
560
561 if (nDecPlaces != 0)
562 {
563 if (nDecPlaces < 0)
564 fValue *= fFac;
565 else
566 fValue /= fFac;
567 }
568
569 if (!std::isfinite(fValue))
570 return fOrigValue;
571
572 return bSign ? -fValue : fValue;
573 }
574
575 double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THROW_EXTERN_C()
576 {
577 return fValue * getN10Exp(nExp);
578 }
579
580 double SAL_CALL rtl_math_approxValue(double fValue) SAL_THROW_EXTERN_C()
581 {
582 const double fBigInt = 0x1p41; // 2^41 -> only 11 bits left for fractional part, fine as decimal
583 if (fValue == 0.0 || !std::isfinite(fValue) || fValue > fBigInt)
584 {
585 // We don't handle these conditions. Bail out.
586 return fValue;
587 }
588
589 double fOrigValue = fValue;
590
591 bool bSign = std::signbit(fValue);
592 if (bSign)
593 fValue = -fValue;
594
595 // If the value is either integer representable as double,
596 // or only has small number of bits in fraction part, then we need not do any approximation
597 if (isRepresentableInteger(fValue) || getBitsInFracPart(fValue) <= 11)
598 return fOrigValue;
599
600 int nExp = static_cast<int>(floor(log10(fValue)));
601 nExp = 14 - nExp;
602 double fExpValue = getN10Exp(abs(nExp));
603
604 if (nExp < 0)
605 fValue /= fExpValue;
606 else
607 fValue *= fExpValue;
608
609 // If the original value was near DBL_MIN we got an overflow. Restore and
610 // bail out.
611 if (!std::isfinite(fValue))
612 return fOrigValue;
613
614 fValue = std::round(fValue);
615
616 if (nExp < 0)
617 fValue *= fExpValue;
618 else
619 fValue /= fExpValue;
620
621 // If the original value was near DBL_MAX we got an overflow. Restore and
622 // bail out.
623 if (!std::isfinite(fValue))
624 return fOrigValue;
625
626 return bSign ? -fValue : fValue;
627 }
628
629 bool SAL_CALL rtl_math_approxEqual(double a, double b) SAL_THROW_EXTERN_C()
630 {
631 static const double e48 = 0x1p-48;
632
633 if (a == b)
634 return true;
635
636 if (a == 0.0 || b == 0.0 || std::signbit(a) != std::signbit(b))
637 return false;
638
639 const double d = fabs(a - b);
640 if (!std::isfinite(d))
641 return false; // Nan or Inf involved
642
643 a = fabs(a);
644 if (d >= (a * e48))
645 return false;
646 b = fabs(b);
647 if (d >= (b * e48))
648 return false;
649
650 if (isRepresentableInteger(a) && isRepresentableInteger(b))
651 return false; // special case for representable integers.
652
653 return true;
654 }
655
656 double SAL_CALL rtl_math_expm1(double fValue) SAL_THROW_EXTERN_C() { return expm1(fValue); }
657
658 double SAL_CALL rtl_math_log1p(double fValue) SAL_THROW_EXTERN_C()
659 {
660 #ifdef __APPLE__
661 if (fValue == -0.0)
662 return fValue; // macOS 10.8 libc returns 0.0 for -0.0
663 #endif
664
665 return log1p(fValue);
666 }
667
668 double SAL_CALL rtl_math_atanh(double fValue) SAL_THROW_EXTERN_C() { return ::atanh(fValue); }
669
670 /** Parent error function (erf) */
671 double SAL_CALL rtl_math_erf(double x) SAL_THROW_EXTERN_C() { return erf(x); }
672
673 /** Parent complementary error function (erfc) */
674 double SAL_CALL rtl_math_erfc(double x) SAL_THROW_EXTERN_C() { return erfc(x); }
675
676 /** improved accuracy of asinh for |x| large and for x near zero
677 @see #i97605#
678 */
679 double SAL_CALL rtl_math_asinh(double fX) SAL_THROW_EXTERN_C()
680 {
681 if (fX == 0.0)
682 return 0.0;
683
684 double fSign = 1.0;
685 if (fX < 0.0)
686 {
687 fX = -fX;
688 fSign = -1.0;
689 }
690
691 if (fX < 0.125)
692 return fSign * rtl_math_log1p(fX + fX * fX / (1.0 + sqrt(1.0 + fX * fX)));
693
694 if (fX < 1.25e7)
695 return fSign * log(fX + sqrt(1.0 + fX * fX));
696
697 return fSign * log(2.0 * fX);
698 }
699
700 /** improved accuracy of acosh for x large and for x near 1
701 @see #i97605#
702 */
703 double SAL_CALL rtl_math_acosh(double fX) SAL_THROW_EXTERN_C()
704 {
705 volatile double fZ = fX - 1.0;
706 if (fX < 1.0)
707 return std::numeric_limits<double>::quiet_NaN();
708 if (fX == 1.0)
709 return 0.0;
710
711 if (fX < 1.1)
712 return rtl_math_log1p(fZ + sqrt(fZ * fZ + 2.0 * fZ));
713
714 if (fX < 1.25e7)
715 return log(fX + sqrt(fX * fX - 1.0));
716
717 return log(2.0 * fX);
718 }
719
720 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
FREE OFFICE SUITE