/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ /* * This file is part of the LibreOffice project. * * This Source Code Form is subject to the terms of the Mozilla Public * License, v. 2.0. If a copy of the MPL was not distributed with this * file, You can obtain one at http://mozilla.org/MPL/2.0/. * * This file incorporates work covered by the following license notice: * * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed * with this work for additional information regarding copyright * ownership. The ASF licenses this file to you under the Apache * License, Version 2.0 (the "License"); you may not use this file * except in compliance with the License. You may obtain a copy of * the License at http://www.apache.org/licenses/LICENSE-2.0 .
*/
// return pow(10.0,nExp) optimized for exponents in the interval [-323,308] (i.e., incl. denormals) staticdouble getN10Exp(int nExp)
{ if (nExp < minExp || nExp > maxExp) return pow(10.0, static_cast<double>(nExp)); // will return 0 or INF with IEEE 754 return n10s[nExp - minExp];
}
namespace
{ /** If value (passed as absolute value) is an integer representable as double, which we handle explicitly at some places.
*/ bool isRepresentableInteger(double fAbsValue)
{
static_assert(std::numeric_limits<double>::is_iec559
&& std::numeric_limits<double>::digits == 53);
assert(fAbsValue >= 0.0); if (fAbsValue >= 0x1p53) returnfalse;
sal_Int64 nInt = static_cast<sal_Int64>(fAbsValue); return nInt == fAbsValue;
}
/** Returns number of binary bits for fractional part of the number Expects a proper non-negative double value, not +-INF, not NAN
*/ int getBitsInFracPart(double fAbsValue)
{
assert(std::isfinite(fAbsValue) && fAbsValue >= 0.0); if (fAbsValue == 0.0) return 0; auto& rValParts = reinterpret_cast<const sal_math_Double*>(&fAbsValue)->parts; int nExponent = rValParts.exponent - 1023; if (nExponent >= 52) return 0; // All bits in fraction are in integer part of the number int nLeastSignificant = rValParts.fraction
? std::countr_zero(rValParts.fraction) + 1
: 53; // the implied leading 1 is the least significant int nFracSignificant = 53 - nLeastSignificant; int nBitsInFracPart = nFracSignificant - nExponent;
if (!bDone) // do not recognize e.g. NaN1.23
{
std::unique_ptr<char[]> bufInHeap;
std::unique_ptr<const CharT* []> bufInHeapMap;
constexpr int bufOnStackSize = 256; char bufOnStack[bufOnStackSize]; const CharT* bufOnStackMap[bufOnStackSize]; char* buf = bufOnStack; const CharT** bufmap = bufOnStackMap; int bufpos = 0; const size_t bufsize = pEnd - p + (bSign ? 2 : 1); if (bufsize > bufOnStackSize)
{
bufInHeap = std::make_unique<char[]>(bufsize);
bufInHeapMap = std::make_unique<const CharT* []>(bufsize);
buf = bufInHeap.get();
bufmap = bufInHeapMap.get();
}
if (bSign)
{
buf[0] = '-';
bufmap[0] = p; // yes, this may be the same pointer as for the next mapping
bufpos = 1;
} // Put first zero to buffer for strings like "-0" if (p != pEnd && *p == '0')
{
buf[bufpos] = '0';
bufmap[bufpos] = p;
++bufpos;
++p;
} // Leading zeros and group separators between digits may be safely // ignored. p0 < p implies that there was a leading 0 already, // consecutive group separators may not happen as *(p+1) is checked for // digit. while (p != pEnd
&& (*p == '0'
|| (*p == cGroupSeparator && p0 < p && p + 1 < pEnd
&& rtl::isAsciiDigit(*(p + 1)))))
{
++p;
}
// integer part of mantissa for (; p != pEnd; ++p)
{
CharT c = *p; if (rtl::isAsciiDigit(c))
{
buf[bufpos] = static_cast<char>(c);
bufmap[bufpos] = p;
++bufpos;
} elseif (c != cGroupSeparator)
{ break;
} elseif (p == p0 || (p + 1 == pEnd) || !rtl::isAsciiDigit(*(p + 1)))
{ // A leading or trailing (not followed by a digit) group // separator character is not a group separator. break;
}
}
// fraction part of mantissa if (p != pEnd && *p == cDecSeparator)
{
buf[bufpos] = '.';
bufmap[bufpos] = p;
++bufpos;
++p;
for (; p != pEnd; ++p)
{
CharT c = *p; if (!rtl::isAsciiDigit(c))
{ break;
}
buf[bufpos] = static_cast<char>(c);
bufmap[bufpos] = p;
++bufpos;
}
}
// Eat any further digits: while (p != pEnd && rtl::isAsciiDigit(*p))
{
++p;
}
bDone = true;
}
}
if (!bDone)
{
buf[bufpos] = '\0';
bufmap[bufpos] = p; char* pCharParseEnd;
errno = 0;
fVal = strtod_nolocale(buf, &pCharParseEnd); if (errno == ERANGE)
{ // This happens with overflow and underflow (including subnormals!) if (std::isinf(fVal))
{ // Check for the dreaded rounded to 15 digits max value // 1.79769313486232e+308 for 1.7976931348623157e+308 we wrote // everywhere, accept with or without plus sign in exponent. constchar* b = buf; if (b[0] == '-')
++b; if (((pCharParseEnd - b == 21) || (pCharParseEnd - b == 20))
&& !strncmp(b, "1.79769313486232", 16) && (b[16] == 'e' || b[16] == 'E')
&& (((pCharParseEnd - b == 21) && !strncmp(b + 17, "+308", 4))
|| ((pCharParseEnd - b == 20) && !strncmp(b + 17, "308", 3))))
{
fVal = (buf < b) ? -DBL_MAX : DBL_MAX;
} else
{
eStatus = rtl_math_ConversionStatus_OutOfRange;
}
} elseif (fVal == 0)
{
eStatus = rtl_math_ConversionStatus_OutOfRange;
} // else it's subnormal: allow it
}
p = bufmap[pCharParseEnd - buf];
bSign = false;
}
}
// overflow also if more than DBL_MAX_10_EXP digits without decimal // separator, or 0. and more than DBL_MIN_10_EXP digits, ... if (std::isinf(fVal))
eStatus = rtl_math_ConversionStatus_OutOfRange;
if (bSign)
fVal = -fVal;
if (pStatus)
*pStatus = eStatus;
if (pParsedEnd)
*pParsedEnd = p == p0 ? pBegin : p;
// Rounding to decimals between integer distance precision (gaps) does not // make sense, do not even try to multiply/divide and introduce inaccuracy. // For same reasons, do not attempt to round integers to decimals. if (nDecPlaces >= 0 && (fValue >= 0x1p52 || isRepresentableInteger(fValue))) return fOrigValue;
double fFac = 0; if (nDecPlaces != 0)
{ if (nDecPlaces > 0)
{ // Determine how many decimals are representable in the precision. // Anything greater 2^52 and 0.0 was already ruled out above. // Theoretically 0.5, 0.25, 0.125, 0.0625, 0.03125, ... const sal_math_Double* pd = reinterpret_cast<const sal_math_Double*>(&fValue); const sal_Int32 nDec = 52 - (pd->parts.exponent - 1023);
if (nDec <= 0)
{
assert(!"Shouldn't this had been caught already as large number?"); return fOrigValue;
}
if (nDec < nDecPlaces)
nDecPlaces = nDec;
}
// Avoid 1e-5 (1.0000000000000001e-05) and such inaccurate fractional // factors that later when dividing back spoil things. For negative // decimals divide first with the inverse, then multiply the rounded // value back.
fFac = getN10Exp(abs(nDecPlaces));
if (fFac == 0.0 || (nDecPlaces < 0 && !std::isfinite(fFac))) // Underflow, rounding to that many integer positions would be 0. return 0.0;
if (!std::isfinite(fFac)) // Overflow with very small values and high number of decimals. return fOrigValue;
double SAL_CALL rtl_math_approxValue(double fValue) noexcept
{ constdouble fBigInt = 0x1p41; // 2^41 -> only 11 bits left for fractional part, fine as decimal if (fValue == 0.0 || !std::isfinite(fValue) || fValue > fBigInt)
{ // We don't handle these conditions. Bail out. return fValue;
}
double fOrigValue = fValue;
bool bSign = std::signbit(fValue); if (bSign)
fValue = -fValue;
// If the value is either integer representable as double, // or only has small number of bits in fraction part, then we need not do any approximation if (isRepresentableInteger(fValue) || getBitsInFracPart(fValue) <= 11) return fOrigValue;
// If the original value was near DBL_MAX we got an overflow. Restore and // bail out. if (!std::isfinite(fValue)) return fOrigValue;
return bSign ? -fValue : fValue;
}
bool SAL_CALL rtl_math_approxEqual(double a, double b) noexcept
{ // threshold is last four bits of mantissa: staticconstdouble e48 = 0x1p-48; // threshold is half of the 15th significant decimal (see doubleToString): staticconstdouble half15thSignificand = 5E-15;
if (a == b) returntrue;
if (a == 0.0 || b == 0.0 || std::signbit(a) != std::signbit(b)) returnfalse;
constdouble d = fabs(a - b); if (!std::isfinite(d)) returnfalse; // Nan or Inf involved
a = fabs(a);
b = fabs(b); constdouble min_ab = std::min(a, b); constdouble threshold1 = min_ab * e48; constdouble threshold2 = getN10Exp(floor(log10(min_ab))) * half15thSignificand; if (d >= std::max(threshold1, threshold2)) returnfalse;
if (isRepresentableInteger(a) && isRepresentableInteger(b)) returnfalse; // special case for representable integers.
/** improved accuracy of asinh for |x| large and for x near zero @see #i97605#
*/ double SAL_CALL rtl_math_asinh(double fX) noexcept
{ if (fX == 0.0) return 0.0;
double fSign = 1.0; if (fX < 0.0)
{
fX = -fX;
fSign = -1.0;
}
if (fX < 0.125) return fSign * rtl_math_log1p(fX + fX * fX / (1.0 + sqrt(1.0 + fX * fX)));
if (fX < 1.25e7) return fSign * log(fX + sqrt(1.0 + fX * fX));
return fSign * log(2.0 * fX);
}
/** improved accuracy of acosh for x large and for x near 1 @see #i97605#
*/ double SAL_CALL rtl_math_acosh(double fX) noexcept
{ double fZ = fX - 1.0; if (fX < 1.0) return std::numeric_limits<double>::quiet_NaN(); if (fX == 1.0) return 0.0;
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.
