diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -110,6 +110,7 @@ libc.src.math.ceil libc.src.math.ceilf libc.src.math.ceill + libc.src.math.coshf libc.src.math.cosf libc.src.math.expf libc.src.math.exp2f @@ -185,6 +186,7 @@ libc.src.math.roundf libc.src.math.roundl libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -129,6 +129,7 @@ libc.src.math.ceil libc.src.math.ceilf libc.src.math.ceill + libc.src.math.coshf libc.src.math.cosf libc.src.math.expf libc.src.math.exp2f @@ -204,6 +205,7 @@ libc.src.math.roundf libc.src.math.roundl libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -135,6 +135,7 @@ libc.src.math.ceilf libc.src.math.ceill libc.src.math.cos + libc.src.math.coshf libc.src.math.cosf libc.src.math.expf libc.src.math.exp2f @@ -211,6 +212,7 @@ libc.src.math.roundl libc.src.math.sin libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -113,6 +113,7 @@ libc.src.math.ceill libc.src.math.cos libc.src.math.cosf + libc.src.math.coshf libc.src.math.expf libc.src.math.exp2f libc.src.math.expm1f @@ -189,6 +190,7 @@ libc.src.math.sin libc.src.math.sincosf libc.src.math.sinf + libc.src.math.sinhf libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -470,6 +470,9 @@ FunctionSpec<"nextafterf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"nextafter", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"nextafterl", RetValSpec, [ArgSpec, ArgSpec]>, + + FunctionSpec<"coshf", RetValSpec, [ArgSpec]>, + FunctionSpec<"sinhf", RetValSpec, [ArgSpec]>, ] >; diff --git a/libc/src/__support/FPUtil/FPBits.h b/libc/src/__support/FPUtil/FPBits.h --- a/libc/src/__support/FPUtil/FPBits.h +++ b/libc/src/__support/FPUtil/FPBits.h @@ -57,31 +57,33 @@ bits |= mantVal; } - UIntType get_mantissa() const { return bits & FloatProp::MANTISSA_MASK; } + inline UIntType get_mantissa() const { + return bits & FloatProp::MANTISSA_MASK; + } // The function return mantissa with implicit bit set for normal values. constexpr UIntType get_explicit_mantissa() { return (FloatProp::MANTISSA_MASK + 1) | (FloatProp::MANTISSA_MASK & bits); } - void set_unbiased_exponent(UIntType expVal) { + inline void set_unbiased_exponent(UIntType expVal) { expVal = (expVal << (FloatProp::MANTISSA_WIDTH)) & FloatProp::EXPONENT_MASK; bits &= ~(FloatProp::EXPONENT_MASK); bits |= expVal; } - uint16_t get_unbiased_exponent() const { + inline uint16_t get_unbiased_exponent() const { return uint16_t((bits & FloatProp::EXPONENT_MASK) >> (FloatProp::MANTISSA_WIDTH)); } - void set_sign(bool signVal) { + inline void set_sign(bool signVal) { bits |= FloatProp::SIGN_MASK; if (!signVal) bits -= FloatProp::SIGN_MASK; } - bool get_sign() const { return (bits & FloatProp::SIGN_MASK) != 0; } + inline bool get_sign() const { return (bits & FloatProp::SIGN_MASK) != 0; } static_assert(sizeof(T) == sizeof(UIntType), "Data type and integral representation have different sizes."); @@ -150,8 +152,8 @@ static constexpr FPBits neg_zero() { return zero(true); } - static constexpr FPBits inf() { - FPBits bits; + static constexpr FPBits inf(bool sign = false) { + FPBits bits(sign ? FloatProp::SIGN_MASK : UIntType(0)); bits.set_unbiased_exponent(MAX_EXPONENT); return bits; } diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -76,6 +76,7 @@ add_math_entrypoint_object(cos) add_math_entrypoint_object(cosf) +add_math_entrypoint_object(coshf) add_math_entrypoint_object(expf) @@ -181,6 +182,7 @@ add_math_entrypoint_object(sin) add_math_entrypoint_object(sinf) +add_math_entrypoint_object(sinhf) add_math_entrypoint_object(sqrt) add_math_entrypoint_object(sqrtf) diff --git a/libc/src/math/coshf.h b/libc/src/math/coshf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/coshf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for coshf -------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_COSHF_H +#define LLVM_LIBC_SRC_MATH_COSHF_H + +namespace __llvm_libc { + +float coshf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_COSHF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -475,6 +475,7 @@ expf.cpp HDRS ../expf.h + expxf.h DEPENDS .common_constants libc.src.__support.FPUtil.fputil @@ -491,6 +492,7 @@ HDRS ../exp2f.h DEPENDS + .common_constants libc.src.__support.FPUtil.fputil libc.src.__support.FPUtil.polyeval libc.include.math @@ -504,6 +506,7 @@ expm1f.cpp HDRS ../expm1f.h + expxf.h DEPENDS .common_constants libc.src.__support.FPUtil.fputil @@ -1115,4 +1118,35 @@ libc.src.__support.FPUtil.generic.fmod COMPILE_OPTIONS -O3 -) \ No newline at end of file +) + +add_entrypoint_object( + coshf + SRCS + coshf.cpp + HDRS + ../coshf.h + DEPENDS + .common_constants + libc.src.__support.FPUtil.fputil + libc.src.__support.FPUtil.polyeval + libc.include.math + COMPILE_OPTIONS + -O3 +) + +add_entrypoint_object( + sinhf + SRCS + sinhf.cpp + HDRS + ../sinhf.h + DEPENDS + .common_constants + libc.src.__support.FPUtil.fputil + libc.src.__support.FPUtil.polyeval + libc.include.math + COMPILE_OPTIONS + -O3 +) + diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h --- a/libc/src/math/generic/common_constants.h +++ b/libc/src/math/generic/common_constants.h @@ -17,19 +17,11 @@ // Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127. extern const double LOG_F[128]; -// Lookup table for exp(m) with m = -104, ..., 89. -// -104 = floor(log(single precision's min denormal)) -// 89 = ceil(log(single precision's max normal)) -// Table is generated with Sollya as follow: -// > display = hexadecimal; -// > for i from -104 to 89 do { D(exp(i)); }; -extern const double EXP_M1[195]; - -// Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127. -// Table is generated with Sollya as follow: -// > display = hexadecimal; -// > for i from 0 to 127 do { D(exp(i / 128)); }; -extern const double EXP_M2[128]; +static constexpr int EXP_bits_p = 5; +static constexpr int EXP_num_p = 1 << EXP_bits_p; + +// > for i from 0 to 31 do { D(2^(i / 32)) - 1; }; +extern const double EXP_M[EXP_num_p]; } // namespace __llvm_libc diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp --- a/libc/src/math/generic/common_constants.cpp +++ b/libc/src/math/generic/common_constants.cpp @@ -102,128 +102,18 @@ 0x1.58cadb5cd7989p-1, 0x1.5ad404c359f2cp-1, 0x1.5cdb1dc6c1764p-1, 0x1.5ee02a9241675p-1, 0x1.60e32f44788d8p-1}; -// Lookup table for exp(m) with m = -104, ..., 89. -// -104 = floor(log(single precision's min denormal)) -// 89 = ceil(log(single precision's max normal)) -// Table is generated with Sollya as follow: -// > display = hexadecimal; -// > for i from -104 to 89 do { D(exp(i)); }; -const double EXP_M1[195] = { - 0x1.f1e6b68529e33p-151, 0x1.525be4e4e601dp-149, 0x1.cbe0a45f75eb1p-148, - 0x1.3884e838aea68p-146, 0x1.a8c1f14e2af5dp-145, 0x1.20a717e64a9bdp-143, - 0x1.8851d84118908p-142, 0x1.0a9bdfb02d240p-140, 0x1.6a5bea046b42ep-139, - 0x1.ec7f3b269efa8p-138, 0x1.4eafb87eab0f2p-136, 0x1.c6e2d05bbc000p-135, - 0x1.35208867c2683p-133, 0x1.a425b317eeacdp-132, 0x1.1d8508fa8246ap-130, - 0x1.840fbc08fdc8ap-129, 0x1.07b7112bc1ffep-127, 0x1.666d0dad2961dp-126, - 0x1.e726c3f64d0fep-125, 0x1.4b0dc07cabf98p-123, 0x1.c1f2daf3b6a46p-122, - 0x1.31c5957a47de2p-120, 0x1.9f96445648b9fp-119, 0x1.1a6baeadb4fd1p-117, - 0x1.7fd974d372e45p-116, 0x1.04da4d1452919p-114, 0x1.62891f06b3450p-113, - 0x1.e1dd273aa8a4ap-112, 0x1.4775e0840bfddp-110, 0x1.bd109d9d94bdap-109, - 0x1.2e73f53fba844p-107, 0x1.9b138170d6bfep-106, 0x1.175af0cf60ec5p-104, - 0x1.7baee1bffa80bp-103, 0x1.02057d1245cebp-101, 0x1.5eafffb34ba31p-100, - 0x1.dca23bae16424p-99, 0x1.43e7fc88b8056p-97, 0x1.b83bf23a9a9ebp-96, - 0x1.2b2b8dd05b318p-94, 0x1.969d47321e4ccp-93, 0x1.1452b7723aed2p-91, - 0x1.778fe2497184cp-90, 0x1.fe7116182e9ccp-89, 0x1.5ae191a99585ap-87, - 0x1.d775d87da854dp-86, 0x1.4063f8cc8bb98p-84, 0x1.b374b315f87c1p-83, - 0x1.27ec458c65e3cp-81, 0x1.923372c67a074p-80, 0x1.1152eaeb73c08p-78, - 0x1.737c5645114b5p-77, 0x1.f8e6c24b5592ep-76, 0x1.571db733a9d61p-74, - 0x1.d257d547e083fp-73, 0x1.3ce9b9de78f85p-71, 0x1.aebabae3a41b5p-70, - 0x1.24b6031b49bdap-68, 0x1.8dd5e1bb09d7ep-67, 0x1.0e5b73d1ff53dp-65, - 0x1.6f741de1748ecp-64, 0x1.f36bd37f42f3ep-63, 0x1.536452ee2f75cp-61, - 0x1.cd480a1b74820p-60, 0x1.39792499b1a24p-58, 0x1.aa0de4bf35b38p-57, - 0x1.2188ad6ae3303p-55, 0x1.898471fca6055p-54, 0x1.0b6c3afdde064p-52, - 0x1.6b7719a59f0e0p-51, 0x1.ee001eed62aa0p-50, 0x1.4fb547c775da8p-48, - 0x1.c8464f7616468p-47, 0x1.36121e24d3bbap-45, 0x1.a56e0c2ac7f75p-44, - 0x1.1e642baeb84a0p-42, 0x1.853f01d6d53bap-41, 0x1.0885298767e9ap-39, - 0x1.67852a7007e42p-38, 0x1.e8a37a45fc32ep-37, 0x1.4c1078fe9228ap-35, - 0x1.c3527e433fab1p-34, 0x1.32b48bf117da2p-32, 0x1.a0db0d0ddb3ecp-31, - 0x1.1b48655f37267p-29, 0x1.81056ff2c5772p-28, 0x1.05a628c699fa1p-26, - 0x1.639e3175a689dp-25, 0x1.e355bbaee85cbp-24, 0x1.4875ca227ec38p-22, - 0x1.be6c6fdb01612p-21, 0x1.2f6053b981d98p-19, 0x1.9c54c3b43bc8bp-18, - 0x1.18354238f6764p-16, 0x1.7cd79b5647c9bp-15, 0x1.02cf22526545ap-13, - 0x1.5fc21041027adp-12, 0x1.de16b9c24a98fp-11, 0x1.44e51f113d4d6p-9, - 0x1.b993fe00d5376p-8, 0x1.2c155b8213cf4p-6, 0x1.97db0ccceb0afp-5, - 0x1.152aaa3bf81ccp-3, 0x1.78b56362cef38p-2, 0x1.0000000000000p+0, - 0x1.5bf0a8b145769p+1, 0x1.d8e64b8d4ddaep+2, 0x1.415e5bf6fb106p+4, - 0x1.b4c902e273a58p+5, 0x1.28d389970338fp+7, 0x1.936dc5690c08fp+8, - 0x1.122885aaeddaap+10, 0x1.749ea7d470c6ep+11, 0x1.fa7157c470f82p+12, - 0x1.5829dcf950560p+14, 0x1.d3c4488ee4f7fp+15, 0x1.3de1654d37c9ap+17, - 0x1.b00b5916ac955p+18, 0x1.259ac48bf05d7p+20, 0x1.8f0ccafad2a87p+21, - 0x1.0f2ebd0a80020p+23, 0x1.709348c0ea4f9p+24, 0x1.f4f22091940bdp+25, - 0x1.546d8f9ed26e1p+27, 0x1.ceb088b68e804p+28, 0x1.3a6e1fd9eecfdp+30, - 0x1.ab5adb9c43600p+31, 0x1.226af33b1fdc1p+33, 0x1.8ab7fb5475fb7p+34, - 0x1.0c3d3920962c9p+36, 0x1.6c932696a6b5dp+37, 0x1.ef822f7f6731dp+38, - 0x1.50bba3796379ap+40, 0x1.c9aae4631c056p+41, 0x1.370470aec28edp+43, - 0x1.a6b765d8cdf6dp+44, 0x1.1f43fcc4b662cp+46, 0x1.866f34a725782p+47, - 0x1.0953e2f3a1ef7p+49, 0x1.689e221bc8d5bp+50, 0x1.ea215a1d20d76p+51, - 0x1.4d13fbb1a001ap+53, 0x1.c4b334617cc67p+54, 0x1.33a43d282a519p+56, - 0x1.a220d397972ebp+57, 0x1.1c25c88df6862p+59, 0x1.8232558201159p+60, - 0x1.0672a3c9eb871p+62, 0x1.64b41c6d37832p+63, 0x1.e4cf766fe49bep+64, - 0x1.49767bc0483e3p+66, 0x1.bfc951eb8bb76p+67, 0x1.304d6aeca254bp+69, - 0x1.9d97010884251p+70, 0x1.19103e4080b45p+72, 0x1.7e013cd114461p+73, - 0x1.03996528e074cp+75, 0x1.60d4f6fdac731p+76, 0x1.df8c5af17ba3bp+77, - 0x1.45e3076d61699p+79, 0x1.baed16a6e0da7p+80, 0x1.2cffdfebde1a1p+82, - 0x1.9919cabefcb69p+83, 0x1.160345c9953e3p+85, 0x1.79dbc9dc53c66p+86, - 0x1.00c810d464097p+88, 0x1.5d009394c5c27p+89, 0x1.da57de8f107a8p+90, - 0x1.425982cf597cdp+92, 0x1.b61e5ca3a5e31p+93, 0x1.29bb825dfcf87p+95, - 0x1.94a90db0d6fe2p+96, 0x1.12fec759586fdp+98, 0x1.75c1dc469e3afp+99, - 0x1.fbfd219c43b04p+100, 0x1.5936d44e1a146p+102, 0x1.d531d8a7ee79cp+103, - 0x1.3ed9d24a2d51bp+105, 0x1.b15cfe5b6e17bp+106, 0x1.268038c2c0e00p+108, - 0x1.9044a73545d48p+109, 0x1.1002ab6218b38p+111, 0x1.71b3540cbf921p+112, - 0x1.f6799ea9c414ap+113, 0x1.55779b984f3ebp+115, 0x1.d01a210c44aa4p+116, - 0x1.3b63da8e91210p+118, 0x1.aca8d6b0116b8p+119, 0x1.234de9e0c74e9p+121, - 0x1.8bec7503ca477p+122, 0x1.0d0eda9796b90p+124, 0x1.6db0118477245p+125, - 0x1.f1056dc7bf22dp+126, 0x1.51c2cc3433801p+128, 0x1.cb108ffbec164p+129, -}; - -// Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127. -// Table is generated with Sollya as follow: -// > display = hexadecimal; -// > for i from 0 to 127 do { D(exp(i / 128)); }; -const double EXP_M2[128] = { - 0x1.0000000000000p0, 0x1.0202015600446p0, 0x1.04080ab55de39p0, - 0x1.06122436410ddp0, 0x1.08205601127edp0, 0x1.0a32a84e9c1f6p0, - 0x1.0c49236829e8cp0, 0x1.0e63cfa7ab09dp0, 0x1.1082b577d34edp0, - 0x1.12a5dd543ccc5p0, 0x1.14cd4fc989cd6p0, 0x1.16f9157587069p0, - 0x1.192937074e0cdp0, 0x1.1b5dbd3f68122p0, 0x1.1d96b0eff0e79p0, - 0x1.1fd41afcba45ep0, 0x1.2216045b6f5cdp0, 0x1.245c7613b8a9bp0, - 0x1.26a7793f60164p0, 0x1.28f7170a755fdp0, 0x1.2b4b58b372c79p0, - 0x1.2da4478b620c7p0, 0x1.3001ecf601af7p0, 0x1.32645269ea829p0, - 0x1.34cb8170b5835p0, 0x1.373783a722012p0, 0x1.39a862bd3c106p0, - 0x1.3c1e2876834aap0, 0x1.3e98deaa11dccp0, 0x1.41188f42c3e32p0, - 0x1.439d443f5f159p0, 0x1.462707b2bac21p0, 0x1.48b5e3c3e8186p0, - 0x1.4b49e2ae5ac67p0, 0x1.4de30ec211e60p0, 0x1.50817263c13cdp0, - 0x1.5325180cfacf7p0, 0x1.55ce0a4c58c7cp0, 0x1.587c53c5a7af0p0, - 0x1.5b2fff3210fd9p0, 0x1.5de9176045ff5p0, 0x1.60a7a734ab0e8p0, - 0x1.636bb9a983258p0, 0x1.663559cf1bc7cp0, 0x1.690492cbf9433p0, - 0x1.6bd96fdd034a2p0, 0x1.6eb3fc55b1e76p0, 0x1.719443a03acb9p0, - 0x1.747a513dbef6ap0, 0x1.776630c678bc1p0, 0x1.7a57ede9ea23ep0, - 0x1.7d4f946f0ba8dp0, 0x1.804d30347b546p0, 0x1.8350cd30ac390p0, - 0x1.865a7772164c5p0, 0x1.896a3b1f66a0ep0, 0x1.8c802477b0010p0, - 0x1.8f9c3fd29beafp0, 0x1.92be99a09bf00p0, 0x1.95e73e6b1b75ep0, - 0x1.99163ad4b1dccp0, 0x1.9c4b9b995509bp0, 0x1.9f876d8e8c566p0, - 0x1.a2c9bda3a3e78p0, 0x1.a61298e1e069cp0, 0x1.a9620c6cb3374p0, - 0x1.acb82581eee54p0, 0x1.b014f179fc3b8p0, 0x1.b3787dc80f95fp0, - 0x1.b6e2d7fa5eb18p0, 0x1.ba540dba56e56p0, 0x1.bdcc2cccd3c85p0, - 0x1.c14b431256446p0, 0x1.c4d15e873c193p0, 0x1.c85e8d43f7cd0p0, - 0x1.cbf2dd7d490f2p0, 0x1.cf8e5d84758a9p0, 0x1.d3311bc7822b4p0, - 0x1.d6db26d16cd67p0, 0x1.da8c8d4a66969p0, 0x1.de455df80e3c0p0, - 0x1.e205a7bdab73ep0, 0x1.e5cd799c6a54ep0, 0x1.e99ce2b397649p0, - 0x1.ed73f240dc142p0, 0x1.f152b7a07bb76p0, 0x1.f539424d90f5ep0, - 0x1.f927a1e24bb76p0, 0x1.fd1de6182f8c9p0, 0x1.008e0f64294abp1, - 0x1.02912df5ce72ap1, 0x1.049856cd84339p1, 0x1.06a39207f0a09p1, - 0x1.08b2e7d2035cfp1, 0x1.0ac6606916501p1, 0x1.0cde041b0e9aep1, - 0x1.0ef9db467dcf8p1, 0x1.1119ee5ac36b6p1, 0x1.133e45d82e952p1, - 0x1.1566ea50201d7p1, 0x1.1793e4652cc50p1, 0x1.19c53ccb3fc6bp1, - 0x1.1bfafc47bda73p1, 0x1.1e352bb1a74adp1, 0x1.2073d3f1bd518p1, - 0x1.22b6fe02a3b9cp1, 0x1.24feb2f105cb8p1, 0x1.274afbdbba4a6p1, - 0x1.299be1f3e7f1cp1, 0x1.2bf16e7d2a38cp1, 0x1.2e4baacdb6614p1, - 0x1.30aaa04e80d05p1, 0x1.330e587b62b28p1, 0x1.3576dce33feadp1, - 0x1.37e437282d4eep1, 0x1.3a5670ff972edp1, 0x1.3ccd9432682b4p1, - 0x1.3f49aa9d30590p1, 0x1.41cabe304cb34p1, 0x1.4450d8f00edd4p1, - 0x1.46dc04f4e5338p1, 0x1.496c4c6b832dap1, 0x1.4c01b9950a111p1, - 0x1.4e9c56c731f5dp1, 0x1.513c2e6c731d7p1, 0x1.53e14b042f9cap1, - 0x1.568bb722dd593p1, 0x1.593b7d72305bbp1, -}; - +// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27] +// printf("%.13a,\n", d[i]); +const double EXP_M[EXP_num_p] = { + -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2, + -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3, + -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3, + -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4, + -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5, + -0x1.5f134923757f3p-6, 0x0.0000000000000p+0, 0x1.66c34c5615d0fp-6, + 0x1.6ab0d9f3121ecp-5, 0x1.1301d0125b50ap-4, 0x1.72b83c7d517aep-4, + 0x1.d4873168b9aa8p-4, 0x1.1c3d373ab11c3p-3, 0x1.4f4efa8fef709p-3, + 0x1.837f0518db8a9p-3, 0x1.b8d39b9d54e55p-3, 0x1.ef5326091a112p-3, + 0x1.13821818624b4p-2, 0x1.2ff6b54d8a89cp-2, 0x1.4d0ad5a753e07p-2, + 0x1.6ac1f752150a5p-2, 0x1.891fac0e95613p-2}; } // namespace __llvm_libc diff --git a/libc/src/math/generic/coshf.cpp b/libc/src/math/generic/coshf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/coshf.cpp @@ -0,0 +1,46 @@ +//===-- Single-precision cosh function ------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/math/coshf.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/generic/expxf.h" + +namespace __llvm_libc { + +LLVM_LIBC_FUNCTION(float, coshf, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + xbits.set_sign(false); + x = xbits.get_val(); + + uint32_t x_u = xbits.uintval(); + + // |x| <= 2^-26 + if (unlikely(x_u <= 0x3280'0000U)) { + return 1.0f + x; + } + + // When |x| >= 90, or x is inf or nan + if (unlikely(x_u >= 0x42b4'0000U)) { + if (xbits.is_inf_or_nan()) + return x + FPBits::inf().get_val(); + + int rounding = fputil::get_round(); + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); + + errno = ERANGE; + + return x + FPBits::inf().get_val(); + } + + double ep = exp_eval(x) + exp_eval(-x); + return 0.5 * ep + 1.0; +} + +} // namespace __llvm_libc diff --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp --- a/libc/src/math/generic/exp2f.cpp +++ b/libc/src/math/generic/exp2f.cpp @@ -7,6 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/exp2f.h" +#include "common_constants.h" #include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FMA.h" @@ -18,34 +19,14 @@ namespace __llvm_libc { -// Lookup table for 2^(m * 2^(-6)) with m = 0, ..., 63. -// Table is generated with Sollya as follow: -// > display = hexadecimal; -// > for i from 0 to 63 do { D(2^(i / 64)); }; -static constexpr double EXP_M[64] = { - 0x1.0000000000000p0, 0x1.02c9a3e778061p0, 0x1.059b0d3158574p0, - 0x1.0874518759bc8p0, 0x1.0b5586cf9890fp0, 0x1.0e3ec32d3d1a2p0, - 0x1.11301d0125b51p0, 0x1.1429aaea92de0p0, 0x1.172b83c7d517bp0, - 0x1.1a35beb6fcb75p0, 0x1.1d4873168b9aap0, 0x1.2063b88628cd6p0, - 0x1.2387a6e756238p0, 0x1.26b4565e27cddp0, 0x1.29e9df51fdee1p0, - 0x1.2d285a6e4030bp0, 0x1.306fe0a31b715p0, 0x1.33c08b26416ffp0, - 0x1.371a7373aa9cbp0, 0x1.3a7db34e59ff7p0, 0x1.3dea64c123422p0, - 0x1.4160a21f72e2ap0, 0x1.44e086061892dp0, 0x1.486a2b5c13cd0p0, - 0x1.4bfdad5362a27p0, 0x1.4f9b2769d2ca7p0, 0x1.5342b569d4f82p0, - 0x1.56f4736b527dap0, 0x1.5ab07dd485429p0, 0x1.5e76f15ad2148p0, - 0x1.6247eb03a5585p0, 0x1.6623882552225p0, 0x1.6a09e667f3bcdp0, - 0x1.6dfb23c651a2fp0, 0x1.71f75e8ec5f74p0, 0x1.75feb564267c9p0, - 0x1.7a11473eb0187p0, 0x1.7e2f336cf4e62p0, 0x1.82589994cce13p0, - 0x1.868d99b4492edp0, 0x1.8ace5422aa0dbp0, 0x1.8f1ae99157736p0, - 0x1.93737b0cdc5e5p0, 0x1.97d829fde4e50p0, 0x1.9c49182a3f090p0, - 0x1.a0c667b5de565p0, 0x1.a5503b23e255dp0, 0x1.a9e6b5579fdbfp0, - 0x1.ae89f995ad3adp0, 0x1.b33a2b84f15fbp0, 0x1.b7f76f2fb5e47p0, - 0x1.bcc1e904bc1d2p0, 0x1.c199bdd85529cp0, 0x1.c67f12e57d14bp0, - 0x1.cb720dcef9069p0, 0x1.d072d4a07897cp0, 0x1.d5818dcfba487p0, - 0x1.da9e603db3285p0, 0x1.dfc97337b9b5fp0, 0x1.e502ee78b3ff6p0, - 0x1.ea4afa2a490dap0, 0x1.efa1bee615a27p0, 0x1.f50765b6e4540p0, - 0x1.fa7c1819e90d8p0, -}; +constexpr double l2h = 0.6931471787393093109130859375; +constexpr double l2l = 1.8206359985041461839581765680755E-9; +constexpr double mlp = double(EXP_num_p); +constexpr double mld = (1.0 / mlp); + +constexpr uint32_t exval1 = 0x3b42'9d37U; +constexpr uint32_t exval2 = 0xbcf3'a937U; +constexpr uint32_t exval_mask = exval1 & exval2; LLVM_LIBC_FUNCTION(float, exp2f, (float x)) { using FPBits = typename fputil::FPBits; @@ -54,36 +35,6 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // Exceptional values. - switch (x_u) { - case 0x3b42'9d37U: // x = 0x1.853a6ep-9f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.00870ap+0f; - break; - case 0x3c02'a9adU: // x = 0x1.05535ap-7f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.016b46p+0f; - break; - case 0x3ca6'6e26U: { // x = 0x1.4cdc4cp-6f - int round_mode = fputil::get_round(); - if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD) - return 0x1.03a16ap+0f; - return 0x1.03a168p+0f; - } - case 0x3d92'a282U: // x = 0x1.254504p-4f - if (fputil::get_round() == FE_UPWARD) - return 0x1.0d0688p+0f; - return 0x1.0d0686p+0f; - case 0xbcf3'a937U: // x = -0x1.e7526ep-6f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.f58d62p-1f; - break; - case 0xb8d3'd026U: // x = -0x1.a7a04cp-14f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.fff6d2p-1f; - break; - } - // // When |x| >= 128, |x| < 2^-25, or x is nan if (unlikely(x_abs >= 0x4300'0000U || x_abs <= 0x3280'0000U)) { // |x| < 2^-25 @@ -101,7 +52,7 @@ errno = ERANGE; } // x is +inf or nan - return x + static_cast(FPBits::inf()); + return x + FPBits::inf().get_val(); } // x < -150 if (x_u >= 0xc316'0000U) { @@ -112,54 +63,46 @@ if (xbits.is_nan()) return x; if (fputil::get_round() == FE_UPWARD) - return static_cast(FPBits(FPBits::MIN_SUBNORMAL)); + return FPBits(FPBits::MIN_SUBNORMAL).get_val(); if (x != 0.0f) errno = ERANGE; return 0.0f; } } - // For -150 <= x < 128, to compute 2^x, we perform the following range - // reduction: find hi, mid, lo such that: - // x = hi + mid + lo, in which - // hi is an integer, - // mid * 2^6 is an integer - // -2^(-7) <= lo < 2^-7. - // In particular, - // hi + mid = round(x * 2^6) * 2^(-6). - // Then, - // 2^(x) = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo. - // Multiply by 2^hi is simply adding hi to the exponent field. We store - // exp(mid) in the lookup tables EXP_M. exp(lo) is computed using a degree-4 - // minimax polynomial generated by Sollya. - // x_hi = round(hi + mid). - // The default rounding mode for float-to-int conversion in C++ is - // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5 - // before conversion. - int x_hi = - static_cast(x * 0x1.0p+6f + (xbits.get_sign() ? -0.5f : 0.5f)); - // For 2-complement integers, arithmetic right shift is the same as dividing - // by a power of 2 and then round down (toward negative infinity). - int e_hi = (x_hi >> 6) + - static_cast(fputil::FloatProperties::EXPONENT_BIAS); - fputil::FPBits exp_hi( - static_cast(e_hi) - << fputil::FloatProperties::MANTISSA_WIDTH); - // mid = x_hi & 0x0000'003fU; - double exp_hi_mid = static_cast(exp_hi) * EXP_M[x_hi & 0x3f]; - // Subtract (hi + mid) from x to get lo. - x -= static_cast(x_hi) * 0x1.0p-6f; - double xd = static_cast(x); - // Degree-4 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > Q = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-7, 2^-7]); - // > Q; - double exp_lo = - fputil::polyeval(xd, 0x1p0, 0x1.62e42fefa2417p-1, 0x1.ebfbdff82f809p-3, - 0x1.c6b0b92131c47p-5, 0x1.3b2ab6fb568a3p-7); - double result = exp_hi_mid * exp_lo; - return static_cast(result); + if (unlikely(x_u & exval_mask) == exval_mask) { + if (unlikely(x_u == exval1)) { // x = 0x1.853a6ep-9f + if (fputil::get_round() == FE_TONEAREST) + return 0x1.00870ap+0f; + } else if (unlikely(x_u == exval2)) { // x = -0x1.e7526ep-6f + if (fputil::get_round() == FE_TONEAREST) + return 0x1.f58d62p-1f; + } + } + + double xdbl = x; + int ps = xdbl * mlp + 16384.5; + double psx = mld * (ps - 16384.0); + double mult_f; + { + ps += 1 << (EXP_bits_p - 1); + int dg = (ps >> EXP_bits_p) - (16384 >> EXP_bits_p); + fputil::FPBits bs; + bs.set_unbiased_exponent(fputil::FPBits::EXPONENT_BIAS + dg); + mult_f = bs.get_val(); + } + + double ml = EXP_M[ps & (EXP_num_p - 1)]; + double dx = xdbl - psx; + // N[Table[Ln[2]^n/n!,{n,1,6}],30] + double pe = dx * fputil::polyeval(dx, l2l, 0.240226506959100712333551263163, + 0.0555041086648215799531422637686, + 0.00961812910762847716197907157366, + 0.00133335581464284434234122219880, + 0.000154035303933816099544370973327) + + l2h * dx; + + return mult_f * (ml * pe + pe + ml + 1.0); } } // namespace __llvm_libc diff --git a/libc/src/math/generic/expf.cpp b/libc/src/math/generic/expf.cpp --- a/libc/src/math/generic/expf.cpp +++ b/libc/src/math/generic/expf.cpp @@ -7,15 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/expf.h" -#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. -#include "src/__support/FPUtil/BasicOperations.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FMA.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/common.h" - -#include +#include "src/math/generic/expxf.h" namespace __llvm_libc { @@ -26,11 +18,6 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // Exceptional values - if (unlikely(x_u == 0xc236'bd8cU)) { // x = -0x1.6d7b18p+5f - return 0x1.108a58p-66f - x * 0x1.0p-95f; - } - // When |x| >= 89, |x| < 2^-25, or x is nan if (unlikely(x_abs >= 0x42b2'0000U || x_abs <= 0x3280'0000U)) { // |x| < 2^-25 @@ -41,66 +28,35 @@ // When x < log(2^-150) or nan if (xbits.uintval() >= 0xc2cf'f1b5U) { // exp(-Inf) = 0 - if (xbits.is_inf()) + if (xbits.is_inf()) { return 0.0f; + } // exp(nan) = nan if (xbits.is_nan()) return x; - if (fputil::get_round() == FE_UPWARD) - return static_cast(FPBits(FPBits::MIN_SUBNORMAL)); + + if (unlikely(fputil::get_round() == FE_UPWARD)) + return FPBits(FPBits::MIN_SUBNORMAL).get_val(); errno = ERANGE; return 0.0f; } + // x >= 89 or nan if (!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000)) { // x is finite if (xbits.uintval() < 0x7f80'0000U) { int rounding = fputil::get_round(); - if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) - return static_cast(FPBits(FPBits::MAX_NORMAL)); + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); errno = ERANGE; } // x is +inf or nan - return x + static_cast(FPBits::inf()); + return x + FPBits::inf().get_val(); } } - // For -104 < x < 89, to compute exp(x), we perform the following range - // reduction: find hi, mid, lo such that: - // x = hi + mid + lo, in which - // hi is an integer, - // mid * 2^7 is an integer - // -2^(-8) <= lo < 2^-8. - // In particular, - // hi + mid = round(x * 2^7) * 2^(-7). - // Then, - // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo). - // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 - // respectively. exp(lo) is computed using a degree-4 minimax polynomial - // generated by Sollya. - // x_hi = (hi + mid) * 2^7 = round(x * 2^7). - // The default rounding mode for float-to-int conversion in C++ is - // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5 - // before conversion. - int x_hi = static_cast(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f)); - // Subtract (hi + mid) from x to get lo. - x -= static_cast(x_hi) * 0x1.0p-7f; - double xd = static_cast(x); - x_hi += 104 << 7; - // hi = x_hi >> 7 - double exp_hi = EXP_M1[x_hi >> 7]; - // mid * 2^7 = x_hi & 0x0000'007fU; - double exp_mid = EXP_M2[x_hi & 0x7f]; - // Degree-4 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]); - // > Q; - double exp_lo = - fputil::polyeval(xd, 0x1p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1, - 0x1.555566668e5e7p-3, 0x1.55555555ef243p-5); - return static_cast(exp_hi * exp_mid * exp_lo); + return exp_eval(x); } } // namespace __llvm_libc diff --git a/libc/src/math/generic/expm1f.cpp b/libc/src/math/generic/expm1f.cpp --- a/libc/src/math/generic/expm1f.cpp +++ b/libc/src/math/generic/expm1f.cpp @@ -7,15 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/expm1f.h" -#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. -#include "src/__support/FPUtil/BasicOperations.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FMA.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/common.h" - -#include +#include "src/math/generic/expxf.h" namespace __llvm_libc { @@ -26,132 +18,44 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // Exceptional value - if (unlikely(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f - int round_mode = fputil::get_round(); - if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD) - return 0x1.8dbe64p-3f; - return 0x1.8dbe62p-3f; - } - -#if !defined(LIBC_TARGET_HAS_FMA) - if (unlikely(x_u == 0xbdc1'c6cbU)) { // x = -0x1.838d96p-4f - int round_mode = fputil::get_round(); - if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD) - return -0x1.71c884p-4f; - return -0x1.71c882p-4f; - } -#endif // LIBC_TARGET_HAS_FMA + // When |x| >= 89, |x| < 2^-25, or x is nan or x < -18 + if (unlikely(x_abs >= 0x42b2'0000U || x_abs <= 0x3280'0000U || + x_u >= 0xc190'0000U)) { + // |x| < 2^-25 + if (xbits.get_unbiased_exponent() <= 101) { + return unlikely(x_abs == 0) ? x : (x + 0.5 * x * x); + } - // When |x| > 25*log(2), or nan - if (unlikely(x_abs >= 0x418a'a123U)) { - // x < log(2^-25) - if (xbits.get_sign()) { + // When x < log(2^-150) or nan + if (xbits.uintval() >= 0xc180'0000U) { // exp(-Inf) = 0 - if (xbits.is_inf()) + if (xbits.is_inf()) { return -1.0f; + } + // exp(nan) = nan if (xbits.is_nan()) return x; - int round_mode = fputil::get_round(); - if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) - return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f - return -1.0f; - } else { - // x >= 89 or nan - if (xbits.uintval() >= 0x42b2'0000) { - if (xbits.uintval() < 0x7f80'0000U) { - int rounding = fputil::get_round(); - if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) - return static_cast(FPBits(FPBits::MAX_NORMAL)); - errno = ERANGE; - } - return x + static_cast(FPBits::inf()); - } + return -1.0f + opt_barrier(FPBits(FPBits::MIN_NORMAL).get_val()); } - } - // |x| < 2^-4 - if (x_abs < 0x3d80'0000U) { - // |x| < 2^-25 - if (x_abs < 0x3300'0000U) { - // x = -0.0f - if (unlikely(xbits.uintval() == 0x8000'0000U)) - return x; - // When |x| < 2^-25, the relative error of the approximation e^x - 1 ~ x - // is: - // |(e^x - 1) - x| / |e^x - 1| < |x^2| / |x| - // = |x| - // < 2^-25 - // < epsilon(1)/2. - // So the correctly rounded values of expm1(x) are: - // = x + eps(x) if rounding mode = FE_UPWARD, - // or (rounding mode = FE_TOWARDZERO and x is - // negative), - // = x otherwise. - // To simplify the rounding decision and make it more efficient, we use - // fma(x, x, x) ~ x + x^2 instead. - // Note: to use the formula x + x^2 to decide the correct rounding, we - // do need fma(x, x, x) to prevent underflow caused by x*x when |x| < - // 2^-76. For targets without FMA instructions, we simply use double for - // intermediate results as it is more efficient than using an emulated - // version of FMA. -#if defined(LIBC_TARGET_HAS_FMA) - return fputil::fma(x, x, x); -#else - double xd = x; - return static_cast(fputil::multiply_add(xd, xd, xd)); -#endif // LIBC_TARGET_HAS_FMA - } + // x >= 89 or nan + if (!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000)) { + // x is finite + if (xbits.uintval() < 0x7f80'0000U) { + int rounding = fputil::get_round(); + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); - // 2^-25 <= |x| < 2^-4 - double xd = static_cast(x); - double xsq = xd * xd; - // Degree-8 minimax polynomial generated by Sollya with: - // > display = hexadecimal; - // > P = fpminimax((expm1(x) - x)/x^2, 6, [|D...|], [-2^-4, 2^-4]); - double r = - fputil::polyeval(xd, 0x1p-1, 0x1.55555555557ddp-3, 0x1.55555555552fap-5, - 0x1.111110fcd58b7p-7, 0x1.6c16c1717660bp-10, - 0x1.a0241f0006d62p-13, 0x1.a01e3f8d3c06p-16); - return static_cast(fputil::multiply_add(r, xsq, xd)); + errno = ERANGE; + } + // x is +inf or nan + return x + FPBits::inf().get_val(); + } } - // For -18 < x < 89, to compute expm1(x), we perform the following range - // reduction: find hi, mid, lo such that: - // x = hi + mid + lo, in which - // hi is an integer, - // mid * 2^7 is an integer - // -2^(-8) <= lo < 2^-8. - // In particular, - // hi + mid = round(x * 2^7) * 2^(-7). - // Then, - // expm1(x) = exp(hi + mid + lo) - 1 = exp(hi) * exp(mid) * exp(lo) - 1. - // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 - // respectively. exp(lo) is computed using a degree-4 minimax polynomial - // generated by Sollya. - - // x_hi = hi + mid. - int x_hi = static_cast(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f)); - // Subtract (hi + mid) from x to get lo. - x -= static_cast(x_hi) * 0x1.0p-7f; - double xd = static_cast(x); - x_hi += 104 << 7; - // hi = x_hi >> 7 - double exp_hi = EXP_M1[x_hi >> 7]; - // lo = x_hi & 0x0000'007fU; - double exp_mid = EXP_M2[x_hi & 0x7f]; - double exp_hi_mid = exp_hi * exp_mid; - // Degree-4 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]); - // > Q; - double exp_lo = - fputil::polyeval(xd, 0x1.0p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1, - 0x1.555566668e5e7p-3, 0x1.55555555ef243p-5); - return static_cast(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0)); + return exp_eval(x); } } // namespace __llvm_libc diff --git a/libc/src/math/generic/expxf.h b/libc/src/math/generic/expxf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/expxf.h @@ -0,0 +1,83 @@ +//===-- Single-precision general exp function -----------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H +#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H + +#include "common_constants.h" // Lookup tables EXP_M +#include "math_utils.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/common.h" +#include + +#include + +namespace __llvm_libc { + +// The algorithm represents exp(x) as +// exp(x) = 2^(ln(2) * i) * 2^(ln(2) * j / NUM_P )) * exp(dx) +// where i integer value, j integer in range [-NUM_P/2, NUM_P/2). +// 2^(ln(2) * j / NUM_P )) is a table values: 1.0 + EXP_M +// exp(dx) calculates by taylor expansion. + +// Inversion of ln(2). Multiplication by EXP_num_p due to sampling by 1 / +// EXP_num_p Precise value of the constant is not needed. +static constexpr double LN2_INV = 1.44269504088896340735 * EXP_num_p; + +// LN2_HIGH + LN2_LOW = ln(2) with precision higher than double(ln(2)) +static constexpr double LN2_HIGH = 0x1.62e42fefa0000p-1 / EXP_num_p; +static constexpr double LN2_LOW = 0x1.cf79abc9e3b3ap-40 / EXP_num_p; + +// Velaue to add to make integer digits always positive. +// Needed for correct rounding to nearest value in both positive and negative +// range. +static constexpr int ADD_TO_POS = (1 << 14); + +// Max possible exp abs value +// 110 is abs of maximum value, which can be passed to "exp_eval" function +static_assert(110 * EXP_num_p < ADD_TO_POS, "Incorrect ADD_TO_POS value."); + +template static inline double exp_eval(double x) { + int64_t ps = fputil::multiply_add(LN2_INV, x, double(ADD_TO_POS) + 0.5); + // Negative sign due to to multiply_add optimization + double psx = double(ADD_TO_POS) - ps; + ps += 1 << (EXP_bits_p - 1); + double mult_e1; + int64_t dg = (ps >> EXP_bits_p) - (ADD_TO_POS >> EXP_bits_p); + { + fputil::FPBits bs; + bs.set_unbiased_exponent(fputil::FPBits::EXPONENT_BIAS + dg); + mult_e1 = bs.get_val(); + } + + double ml = EXP_M[ps & (EXP_num_p - 1)]; + double dx = fputil::multiply_add(psx, LN2_LOW, + fputil::multiply_add(psx, LN2_HIGH, x)); + + // Taylor series coefficients + // double c_x1 = fputil::multiply_add(dx, 1.0 / 2.0, 1.0); + // double c_x3 = fputil::multiply_add(dx, 1.0 / 24.0, 1.0 / 6.0); + // double c_x5 = fputil::multiply_add(dx, 1.0 / 720.0, 1.0 / 120.0); + // double pe = dx * fputil::polyeval(dx * dx, c_x1, c_x3, c_x5); + + double pe = dx * fputil::polyeval(dx, 1.0, 1.0 / 2.0, 1.0 / 6.0, 1.0 / 24.0, + 1.0 / 120.0, 1.0 / 720.0); + + double result = fputil::multiply_add(ml, pe, pe) + ml; + if constexpr (!EXPM1) { + return fputil::multiply_add(mult_e1, result, mult_e1); + } else { + return fputil::multiply_add(mult_e1, result, mult_e1 - 1.0); + } +} + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/sinhf.cpp @@ -0,0 +1,68 @@ +//===-- Single-precision sinh function ------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/math/sinhf.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/generic/expxf.h" + +namespace __llvm_libc { + +LLVM_LIBC_FUNCTION(float, sinhf, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + bool sign = xbits.get_sign(); + uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK; + + // |x| <= 2^-26 + if (unlikely(x_abs <= 0x3280'0000U)) { + return unlikely(x_abs == 0) ? x : (x + 0.25 * x * x * x); + } + + // When |x| >= 90, or x is inf or nan + if (unlikely(x_abs >= 0x42b4'0000U)) { + if (xbits.is_nan()) + return x + 1.0f; // sNaN to qNaN + signal + + if (xbits.is_inf()) + return x; + + int rounding = fputil::get_round(); + if (sign) { + if (unlikely(rounding == FE_UPWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK) + .get_val(); + } else { + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); + } + + errno = ERANGE; + + return x + FPBits::inf(sign).get_val(); + } + + // |x| <= 0.078125 + if (unlikely(x_abs <= 0x3da0'0000U)) { + // |x| = 0.0005589424981735646724700927734375 + if (unlikely(x_abs == 0x3a12'85ffU)) { + if (fputil::get_round() == FE_TONEAREST) + return x; + } + + double xdbl = x; + double x2 = xdbl * xdbl; + // Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x); + double pe = fputil::polyeval(x2, 0.0, 0x1.5555555556583p-3, + 0x1.111110d239f1fp-7, 0x1.a02b5a284013cp-13); + return fputil::multiply_add(xdbl, pe, xdbl); + } + + return 0.5 * (exp_eval(x) - exp_eval(-x)); +} + +} // namespace __llvm_libc diff --git a/libc/src/math/sinhf.h b/libc/src/math/sinhf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/sinhf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for sinhf -------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_SINHF_H +#define LLVM_LIBC_SRC_MATH_SINHF_H + +namespace __llvm_libc { + +float sinhf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINHF_H diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1299,6 +1299,38 @@ libc.src.__support.FPUtil.fputil ) +add_fp_unittest( + coshf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + coshf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + libc.include.errno + libc.src.math.coshf + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fputil +) + +add_fp_unittest( + sinhf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + sinhf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + libc.include.errno + libc.src.math.sinhf + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fputil +) + add_subdirectory(generic) add_subdirectory(exhaustive) add_subdirectory(differential_testing) diff --git a/libc/test/src/math/coshf_test.cpp b/libc/test/src/math/coshf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/coshf_test.cpp @@ -0,0 +1,78 @@ +//===-- Unittests for coshf -----------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/__support/CPP/Array.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/coshf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" +#include + +#include +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +DECLARE_SPECIAL_CONSTANTS(float) + +TEST(LlvmLibcCoshfTest, SpecialNumbers) { + errno = 0; + + EXPECT_FP_EQ(aNaN, __llvm_libc::coshf(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(inf, __llvm_libc::coshf(inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(inf, __llvm_libc::coshf(neg_inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(1.0f, __llvm_libc::coshf(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(1.0f, __llvm_libc::coshf(-0.0f)); + EXPECT_MATH_ERRNO(0); +} + +TEST(LlvmLibcCoshfTest, Overflow) { + errno = 0; + EXPECT_FP_EQ(inf, __llvm_libc::coshf(float(FPBits(0x7f7fffffU)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::coshf(float(FPBits(0x42cffff8U)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::coshf(float(FPBits(0x42d00008U)))); + EXPECT_MATH_ERRNO(ERANGE); +} + +TEST(LlvmLibcCoshfTest, InFloatRange) { + constexpr uint32_t COUNT = 1000000; + constexpr uint32_t STEP = UINT32_MAX / COUNT; + for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { + float x = float(FPBits(v)); + if (isnan(x) || isinf(x)) + continue; + ASSERT_MPFR_MATCH(mpfr::Operation::Cosh, x, __llvm_libc::coshf(x), 1.0); + } +} + +TEST(LlvmLibcCoshfTest, SmallValues) { + float x = float(FPBits(0x17800000U)); + float result = __llvm_libc::coshf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Cosh, x, result, 1.0); + EXPECT_FP_EQ(1.0f, result); + + x = float(FPBits(0x0040000U)); + result = __llvm_libc::coshf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Cosh, x, result, 1.0); + EXPECT_FP_EQ(1.0f, result); +} diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -197,3 +197,38 @@ libc.src.__support.FPUtil.fputil libc.src.__support.FPUtil.generic.fmod ) + +add_fp_unittest( + coshf_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + coshf_test.cpp + DEPENDS + .exhaustive_test + libc.include.math + libc.src.math.coshf + libc.src.__support.FPUtil.fputil + LINK_LIBRARIES + -lpthread +) + +add_fp_unittest( + sinhf_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + sinhf_test.cpp + DEPENDS + .exhaustive_test + libc.include.math + libc.src.math.sinhf + libc.src.__support.FPUtil.fputil + LINK_LIBRARIES + -lpthread +) + diff --git a/libc/test/src/math/exhaustive/coshf_test.cpp b/libc/test/src/math/exhaustive/coshf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exhaustive/coshf_test.cpp @@ -0,0 +1,60 @@ +//===-- Exhaustive test for expf ------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "exhaustive_test.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/coshf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" + +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcCoshfExhaustiveTest : public LlvmLibcExhaustiveTest { + bool check(uint32_t start, uint32_t stop, + mpfr::RoundingMode rounding) override { + mpfr::ForceRoundingMode r(rounding); + uint32_t bits = start; + bool result = true; + do { + FPBits xbits(bits); + float x = float(xbits); + result &= EXPECT_MPFR_MATCH(mpfr::Operation::Cosh, x, + __llvm_libc::coshf(x), 0.5, rounding); + } while (bits++ < stop); + return result; + } +}; + +static const int NUM_THREADS = std::thread::hardware_concurrency(); + +// Range: [0, 90]; +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x42b4'0000U; + +TEST_F(LlvmLibcCoshfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcCoshfExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcCoshfExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcCoshfExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} \ No newline at end of file diff --git a/libc/test/src/math/exhaustive/sinhf_test.cpp b/libc/test/src/math/exhaustive/sinhf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exhaustive/sinhf_test.cpp @@ -0,0 +1,83 @@ +//===-- Exhaustive test for expf ------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "exhaustive_test.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/sinhf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" + +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcSinhfExhaustiveTest : public LlvmLibcExhaustiveTest { + bool check(uint32_t start, uint32_t stop, + mpfr::RoundingMode rounding) override { + mpfr::ForceRoundingMode r(rounding); + uint32_t bits = start; + bool result = true; + do { + FPBits xbits(bits); + float x = float(xbits); + result &= EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, + __llvm_libc::sinhf(x), 0.5, rounding); + } while (bits++ < stop); + return result; + } +}; + +static const int NUM_THREADS = std::thread::hardware_concurrency(); + +// Range: [0, 90]; +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x42b4'0000U; + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} + +// Range: [-90, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xc2b4'0000U; + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundUp) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundDown) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} diff --git a/libc/test/src/math/exp2f_test.cpp b/libc/test/src/math/exp2f_test.cpp --- a/libc/test/src/math/exp2f_test.cpp +++ b/libc/test/src/math/exp2f_test.cpp @@ -54,18 +54,18 @@ TEST(LlvmLibcExp2fTest, TrickyInputs) { constexpr int N = 12; constexpr uint32_t INPUTS[N] = { - 0x3b429d37U, /*0x1.853a6ep-9f*/ - 0x3c02a9adU, /*0x1.05535ap-7f*/ - 0x3ca66e26U, /*0x1.4cdc4cp-6f*/ - 0x3d92a282U, /*0x1.254504p-4f*/ - 0x42fa0001U, /*0x1.f40002p+6f*/ - 0x42ffffffU, /*0x1.fffffep+6f*/ - 0xb8d3d026U, /*-0x1.a7a04cp-14f*/ - 0xbcf3a937U, /*-0x1.e7526ep-6f*/ - 0xc2fa0001U, /*-0x1.f40002p+6f*/ - 0xc2fc0000U, /*-0x1.f8p+6f*/ - 0xc2fc0001U, /*-0x1.f80002p+6f*/ - 0xc3150000U, /*-0x1.2ap+7f*/ + 0x3b429d37U, // 0x1.853a6ep-9f + 0x3c02a9adU, // 0x1.05535ap-7f + 0x3ca66e26U, // 0x1.4cdc4cp-6f + 0x3d92a282U, // 0x1.254504p-4f + 0x42fa0001U, // 0x1.f40002p+6f + 0x42ffffffU, // 0x1.fffffep+6f + 0xb8d3d026U, // -0x1.a7a04cp-14f + 0xbcf3a937U, // -0x1.e7526ep-6f + 0xc2fa0001U, // -0x1.f40002p+6f + 0xc2fc0000U, // -0x1.f8p+6f + 0xc2fc0001U, // -0x1.f80002p+6f + 0xc3150000U, // -0x1.2ap+7f }; for (int i = 0; i < N; ++i) { errno = 0; diff --git a/libc/test/src/math/expf_test.cpp b/libc/test/src/math/expf_test.cpp --- a/libc/test/src/math/expf_test.cpp +++ b/libc/test/src/math/expf_test.cpp @@ -73,6 +73,22 @@ float x; errno = 0; + x = float(FPBits(0xbbf0edf1U)); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, __llvm_libc::expf(x), + 0.5); + EXPECT_MATH_ERRNO(0); + + errno = 0; + x = float(FPBits(0xbbf0'edf1U)); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, __llvm_libc::expf(x), + 0.5); + EXPECT_MATH_ERRNO(0); + + x = float(FPBits(0xc169'12cdU)); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, __llvm_libc::expf(x), + 0.5); + EXPECT_MATH_ERRNO(0); + x = float(FPBits(0x42affff8U)); ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, __llvm_libc::expf(x), 0.5); diff --git a/libc/test/src/math/expm1f_test.cpp b/libc/test/src/math/expm1f_test.cpp --- a/libc/test/src/math/expm1f_test.cpp +++ b/libc/test/src/math/expm1f_test.cpp @@ -68,6 +68,17 @@ float x; errno = 0; + x = -120.0; + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x, + __llvm_libc::expm1f(x), 0.5); + EXPECT_MATH_ERRNO(0); + + errno = 0; + x = float(FPBits(0x989e18e7U)); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x, + __llvm_libc::expm1f(x), 0.5); + EXPECT_MATH_ERRNO(0); + x = float(FPBits(0x42affff8U)); ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 0.5); diff --git a/libc/test/src/math/sinhf_test.cpp b/libc/test/src/math/sinhf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/sinhf_test.cpp @@ -0,0 +1,79 @@ +//===-- Unittests for sinhf -----------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/__support/CPP/Array.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/sinhf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" +#include + +#include +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +DECLARE_SPECIAL_CONSTANTS(float) + +TEST(LlvmLibcSinhfTest, SpecialNumbers) { + errno = 0; + + EXPECT_FP_EQ(aNaN, __llvm_libc::sinhf(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(0.0f, __llvm_libc::sinhf(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(-0.0f, __llvm_libc::sinhf(-0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(neg_inf, __llvm_libc::sinhf(neg_inf)); + EXPECT_MATH_ERRNO(0); +} + +TEST(LlvmLibcSinhfTest, InFloatRange) { + constexpr uint32_t COUNT = 1000000; + constexpr uint32_t STEP = UINT32_MAX / COUNT; + for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { + float x = float(FPBits(v)); + if (isnan(x) || isinf(x)) + continue; + ASSERT_MPFR_MATCH(mpfr::Operation::Sinh, x, __llvm_libc::sinhf(x), 1.0); + } +} + +// For small values, sinh(x) is x. +TEST(LlvmLibcSinhfTest, SmallValues) { + float x = float(FPBits(uint32_t(0x17800000))); + float result = __llvm_libc::sinhf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 1.0); + EXPECT_FP_EQ(x, result); + + x = float(FPBits(uint32_t(0x00400000))); + result = __llvm_libc::sinhf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 1.0); + EXPECT_FP_EQ(x, result); +} + +TEST(LlvmLibcExpfTest, Overflow) { + errno = 0; + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x7f7fffffU)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42cffff8U)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42d00008U)))); + EXPECT_MATH_ERRNO(ERANGE); +} \ No newline at end of file diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -26,6 +26,7 @@ Abs, Ceil, Cos, + Cosh, Exp, Exp2, Expm1, @@ -39,6 +40,7 @@ ModPIOver4, Round, Sin, + Sinh, Sqrt, Tan, Trunc, diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -223,6 +223,12 @@ return result; } + MPFRNumber cosh() const { + MPFRNumber result(*this); + mpfr_cosh(result.value, value, mpfr_rounding); + return result; + } + MPFRNumber exp() const { MPFRNumber result(*this); mpfr_exp(result.value, value, mpfr_rounding); @@ -360,6 +366,12 @@ return result; } + MPFRNumber sinh() const { + MPFRNumber result(*this); + mpfr_sinh(result.value, value, mpfr_rounding); + return result; + } + MPFRNumber sqrt() const { MPFRNumber result(*this); mpfr_sqrt(result.value, value, mpfr_rounding); @@ -510,6 +522,8 @@ return mpfrInput.ceil(); case Operation::Cos: return mpfrInput.cos(); + case Operation::Cosh: + return mpfrInput.cosh(); case Operation::Exp: return mpfrInput.exp(); case Operation::Exp2: @@ -536,6 +550,8 @@ return mpfrInput.round(); case Operation::Sin: return mpfrInput.sin(); + case Operation::Sinh: + return mpfrInput.sinh(); case Operation::Sqrt: return mpfrInput.sqrt(); case Operation::Tan: