[libc-commits] [libc] 5fedf7f - [libc] Move implementations of cosf, sinf, sincosf to src/math directory.
Siva Chandra Reddy via libc-commits
libc-commits at lists.llvm.org
Thu Apr 16 08:46:34 PDT 2020
Author: Siva Chandra Reddy
Date: 2020-04-16T08:46:10-07:00
New Revision: 5fedf7f42043f6a7d4562df2eab4a22b3346ac1a
URL: https://github.com/llvm/llvm-project/commit/5fedf7f42043f6a7d4562df2eab4a22b3346ac1a
DIFF: https://github.com/llvm/llvm-project/commit/5fedf7f42043f6a7d4562df2eab4a22b3346ac1a.diff
LOG: [libc] Move implementations of cosf, sinf, sincosf to src/math directory.
NFC intended in the implementaton. Only mechanical changes to fit the LLVM
libc implementation standard have been done.
Math testing infrastructure has been added. This infrastructure compares the
results produced by the libc with the high precision results from MPFR.
Tests making use of this infrastructure have been added for cosf, sinf and
sincosf.
Reviewers: abrachet, phosek
Differential Revision: https://reviews.llvm.org/D76825
Added:
libc/src/math/cosf.cpp
libc/src/math/cosf.h
libc/src/math/math_utils.h
libc/src/math/sincosf.cpp
libc/src/math/sincosf.h
libc/src/math/sincosf_data.cpp
libc/src/math/sincosf_utils.h
libc/src/math/sinf.cpp
libc/src/math/sinf.h
libc/test/src/math/CMakeLists.txt
libc/test/src/math/cosf_test.cpp
libc/test/src/math/float.h
libc/test/src/math/sdcomp26094.h
libc/test/src/math/sincosf_test.cpp
libc/test/src/math/sinf_test.cpp
libc/utils/MPFRWrapper/CMakeLists.txt
libc/utils/MPFRWrapper/MPFRUtils.cpp
libc/utils/MPFRWrapper/MPFRUtils.h
libc/utils/MPFRWrapper/check_mpfr.cpp
Modified:
libc/AOR_v20.02/math/test/testcases/random/float.tst
libc/config/linux/api.td
libc/lib/CMakeLists.txt
libc/src/__support/common.h.def
libc/src/math/CMakeLists.txt
libc/test/src/CMakeLists.txt
libc/utils/CMakeLists.txt
Removed:
libc/AOR_v20.02/math/cosf.c
libc/AOR_v20.02/math/sincosf.c
libc/AOR_v20.02/math/sincosf.h
libc/AOR_v20.02/math/sincosf_data.c
libc/AOR_v20.02/math/sinf.c
libc/AOR_v20.02/math/test/testcases/directed/cosf.tst
libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst
libc/AOR_v20.02/math/test/testcases/directed/sinf.tst
################################################################################
diff --git a/libc/AOR_v20.02/math/cosf.c b/libc/AOR_v20.02/math/cosf.c
deleted file mode 100644
index 1ab98a1dd60f..000000000000
--- a/libc/AOR_v20.02/math/cosf.c
+++ /dev/null
@@ -1,64 +0,0 @@
-/*
- * Single-precision cos 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 <stdint.h>
-#include <math.h>
-#include "math_config.h"
-#include "sincosf.h"
-
-/* Fast cosf implementation. Worst-case ULP is 0.5607, maximum relative
- error is 0.5303 * 2^-23. A single-step range reduction is used for
- small values. Large inputs have their range reduced using fast integer
- arithmetic. */
-float
-cosf (float y)
-{
- double x = y;
- double s;
- int n;
- const sincos_t *p = &__sincosf_table[0];
-
- if (abstop12 (y) < abstop12 (pio4))
- {
- double x2 = x * x;
-
- if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
- return 1.0f;
-
- return sinf_poly (x, x2, p, 1);
- }
- else if (likely (abstop12 (y) < abstop12 (120.0f)))
- {
- x = reduce_fast (x, p, &n);
-
- /* Setup the signs for sin and cos. */
- s = p->sign[n & 3];
-
- if (n & 2)
- p = &__sincosf_table[1];
-
- return sinf_poly (x * s, x * x, p, n ^ 1);
- }
- else if (abstop12 (y) < abstop12 (INFINITY))
- {
- uint32_t xi = asuint (y);
- int sign = xi >> 31;
-
- x = reduce_large (xi, &n);
-
- /* Setup signs for sin and cos - include original sign. */
- s = p->sign[(n + sign) & 3];
-
- if ((n + sign) & 2)
- p = &__sincosf_table[1];
-
- return sinf_poly (x * s, x * x, p, n ^ 1);
- }
- else
- return __math_invalidf (y);
-}
diff --git a/libc/AOR_v20.02/math/sincosf.c b/libc/AOR_v20.02/math/sincosf.c
deleted file mode 100644
index 819b05b21080..000000000000
--- a/libc/AOR_v20.02/math/sincosf.c
+++ /dev/null
@@ -1,80 +0,0 @@
-/*
- * Single-precision sin/cos 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 <stdint.h>
-#include <math.h>
-#include "math_config.h"
-#include "sincosf.h"
-
-/* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative
- error is 0.5303 * 2^-23. A single-step range reduction is used for
- small values. Large inputs have their range reduced using fast integer
- arithmetic. */
-void
-sincosf (float y, float *sinp, float *cosp)
-{
- double x = y;
- double s;
- int n;
- const sincos_t *p = &__sincosf_table[0];
-
- if (abstop12 (y) < abstop12 (pio4))
- {
- double x2 = x * x;
-
- if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
- {
- if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
- /* Force underflow for tiny y. */
- force_eval_float (x2);
- *sinp = y;
- *cosp = 1.0f;
- return;
- }
-
- sincosf_poly (x, x2, p, 0, sinp, cosp);
- }
- else if (abstop12 (y) < abstop12 (120.0f))
- {
- x = reduce_fast (x, p, &n);
-
- /* Setup the signs for sin and cos. */
- s = p->sign[n & 3];
-
- if (n & 2)
- p = &__sincosf_table[1];
-
- sincosf_poly (x * s, x * x, p, n, sinp, cosp);
- }
- else if (likely (abstop12 (y) < abstop12 (INFINITY)))
- {
- uint32_t xi = asuint (y);
- int sign = xi >> 31;
-
- x = reduce_large (xi, &n);
-
- /* Setup signs for sin and cos - include original sign. */
- s = p->sign[(n + sign) & 3];
-
- if ((n + sign) & 2)
- p = &__sincosf_table[1];
-
- sincosf_poly (x * s, x * x, p, n, sinp, cosp);
- }
- else
- {
- /* Return NaN if Inf or NaN for both sin and cos. */
- *sinp = *cosp = y - y;
-#if WANT_ERRNO
- /* Needed to set errno for +-Inf, the add is a hack to work
- around a gcc register allocation issue: just passing y
- affects code generation in the fast path. */
- __math_invalidf (y + y);
-#endif
- }
-}
diff --git a/libc/AOR_v20.02/math/sincosf.h b/libc/AOR_v20.02/math/sincosf.h
deleted file mode 100644
index ef40d708acde..000000000000
--- a/libc/AOR_v20.02/math/sincosf.h
+++ /dev/null
@@ -1,154 +0,0 @@
-/*
- * Header for sinf, cosf and sincosf.
- *
- * 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 <stdint.h>
-#include <math.h>
-#include "math_config.h"
-
-/* 2PI * 2^-64. */
-static const double pi63 = 0x1.921FB54442D18p-62;
-/* PI / 4. */
-static const double pio4 = 0x1.921FB54442D18p-1;
-
-/* The constants and polynomials for sine and cosine. */
-typedef struct
-{
- double sign[4]; /* Sign of sine in quadrants 0..3. */
- double hpi_inv; /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS). */
- double hpi; /* PI / 2. */
- double c0, c1, c2, c3, c4; /* Cosine polynomial. */
- double s1, s2, s3; /* Sine polynomial. */
-} sincos_t;
-
-/* Polynomial data (the cosine polynomial is negated in the 2nd entry). */
-extern const sincos_t __sincosf_table[2] HIDDEN;
-
-/* Table with 4/PI to 192 bit precision. */
-extern const uint32_t __inv_pio4[] HIDDEN;
-
-/* Top 12 bits of the float representation with the sign bit cleared. */
-static inline uint32_t
-abstop12 (float x)
-{
- return (asuint (x) >> 20) & 0x7ff;
-}
-
-/* Compute the sine and cosine of inputs X and X2 (X squared), using the
- polynomial P and store the results in SINP and COSP. N is the quadrant,
- if odd the cosine and sine polynomials are swapped. */
-static inline void
-sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
- float *cosp)
-{
- double x3, x4, x5, x6, s, c, c1, c2, s1;
-
- x4 = x2 * x2;
- x3 = x2 * x;
- c2 = p->c3 + x2 * p->c4;
- s1 = p->s2 + x2 * p->s3;
-
- /* Swap sin/cos result based on quadrant. */
- float *tmp = (n & 1 ? cosp : sinp);
- cosp = (n & 1 ? sinp : cosp);
- sinp = tmp;
-
- c1 = p->c0 + x2 * p->c1;
- x5 = x3 * x2;
- x6 = x4 * x2;
-
- s = x + x3 * p->s1;
- c = c1 + x4 * p->c2;
-
- *sinp = s + x5 * s1;
- *cosp = c + x6 * c2;
-}
-
-/* Return the sine of inputs X and X2 (X squared) using the polynomial P.
- N is the quadrant, and if odd the cosine polynomial is used. */
-static inline float
-sinf_poly (double x, double x2, const sincos_t *p, int n)
-{
- double x3, x4, x6, x7, s, c, c1, c2, s1;
-
- if ((n & 1) == 0)
- {
- x3 = x * x2;
- s1 = p->s2 + x2 * p->s3;
-
- x7 = x3 * x2;
- s = x + x3 * p->s1;
-
- return s + x7 * s1;
- }
- else
- {
- x4 = x2 * x2;
- c2 = p->c3 + x2 * p->c4;
- c1 = p->c0 + x2 * p->c1;
-
- x6 = x4 * x2;
- c = c1 + x4 * p->c2;
-
- return c + x6 * c2;
- }
-}
-
-/* Fast range reduction using single multiply-subtract. Return the modulo of
- X as a value between -PI/4 and PI/4 and store the quadrant in NP.
- The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
- is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
- the result is accurate for |X| <= 120.0. */
-static inline double
-reduce_fast (double x, const sincos_t *p, int *np)
-{
- double r;
-#if TOINT_INTRINSICS
- /* Use fast round and lround instructions when available. */
- r = x * p->hpi_inv;
- *np = converttoint (r);
- return x - roundtoint (r) * p->hpi;
-#else
- /* Use scaled float to int conversion with explicit rounding.
- hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
- This avoids inaccuracies introduced by truncating negative values. */
- r = x * p->hpi_inv;
- int n = ((int32_t)r + 0x800000) >> 24;
- *np = n;
- return x - n * p->hpi;
-#endif
-}
-
-/* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
- XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
- Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
- Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit
- multiply computes the exact 2.62-bit fixed-point modulo. Since the result
- can have at most 29 leading zeros after the binary point, the double
- precision result is accurate to 33 bits. */
-static inline double
-reduce_large (uint32_t xi, int *np)
-{
- const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
- int shift = (xi >> 23) & 7;
- uint64_t n, res0, res1, res2;
-
- xi = (xi & 0xffffff) | 0x800000;
- xi <<= shift;
-
- res0 = xi * arr[0];
- res1 = (uint64_t)xi * arr[4];
- res2 = (uint64_t)xi * arr[8];
- res0 = (res2 >> 32) | (res0 << 32);
- res0 += res1;
-
- n = (res0 + (1ULL << 61)) >> 62;
- res0 -= n << 62;
- double x = (int64_t)res0;
- *np = n;
- return x * pi63;
-}
diff --git a/libc/AOR_v20.02/math/sincosf_data.c b/libc/AOR_v20.02/math/sincosf_data.c
deleted file mode 100644
index a3d5efda6f4c..000000000000
--- a/libc/AOR_v20.02/math/sincosf_data.c
+++ /dev/null
@@ -1,64 +0,0 @@
-/*
- * Data definition for sinf, cosf and sincosf.
- *
- * 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 <stdint.h>
-#include <math.h>
-#include "math_config.h"
-#include "sincosf.h"
-
-/* The constants and polynomials for sine and cosine. The 2nd entry
- computes -cos (x) rather than cos (x) to get negation for free. */
-const sincos_t __sincosf_table[2] =
-{
- {
- { 1.0, -1.0, -1.0, 1.0 },
-#if TOINT_INTRINSICS
- 0x1.45F306DC9C883p-1,
-#else
- 0x1.45F306DC9C883p+23,
-#endif
- 0x1.921FB54442D18p0,
- 0x1p0,
- -0x1.ffffffd0c621cp-2,
- 0x1.55553e1068f19p-5,
- -0x1.6c087e89a359dp-10,
- 0x1.99343027bf8c3p-16,
- -0x1.555545995a603p-3,
- 0x1.1107605230bc4p-7,
- -0x1.994eb3774cf24p-13
- },
- {
- { 1.0, -1.0, -1.0, 1.0 },
-#if TOINT_INTRINSICS
- 0x1.45F306DC9C883p-1,
-#else
- 0x1.45F306DC9C883p+23,
-#endif
- 0x1.921FB54442D18p0,
- -0x1p0,
- 0x1.ffffffd0c621cp-2,
- -0x1.55553e1068f19p-5,
- 0x1.6c087e89a359dp-10,
- -0x1.99343027bf8c3p-16,
- -0x1.555545995a603p-3,
- 0x1.1107605230bc4p-7,
- -0x1.994eb3774cf24p-13
- }
-};
-
-/* Table with 4/PI to 192 bit precision. To avoid unaligned accesses
- only 8 new bits are added per entry, making the table 4 times larger. */
-const uint32_t __inv_pio4[24] =
-{
- 0xa2, 0xa2f9, 0xa2f983, 0xa2f9836e,
- 0xf9836e4e, 0x836e4e44, 0x6e4e4415, 0x4e441529,
- 0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1,
- 0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0,
- 0x34ddc0db, 0xddc0db62, 0xc0db6295, 0xdb629599,
- 0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041
-};
diff --git a/libc/AOR_v20.02/math/sinf.c b/libc/AOR_v20.02/math/sinf.c
deleted file mode 100644
index 644b82dd94da..000000000000
--- a/libc/AOR_v20.02/math/sinf.c
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * Single-precision sin 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 <math.h>
-#include "math_config.h"
-#include "sincosf.h"
-
-/* Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative
- error is 0.5303 * 2^-23. A single-step range reduction is used for
- small values. Large inputs have their range reduced using fast integer
- arithmetic. */
-float
-sinf (float y)
-{
- double x = y;
- double s;
- int n;
- const sincos_t *p = &__sincosf_table[0];
-
- if (abstop12 (y) < abstop12 (pio4))
- {
- s = x * x;
-
- if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
- {
- if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
- /* Force underflow for tiny y. */
- force_eval_float (s);
- return y;
- }
-
- return sinf_poly (x, s, p, 0);
- }
- else if (likely (abstop12 (y) < abstop12 (120.0f)))
- {
- x = reduce_fast (x, p, &n);
-
- /* Setup the signs for sin and cos. */
- s = p->sign[n & 3];
-
- if (n & 2)
- p = &__sincosf_table[1];
-
- return sinf_poly (x * s, x * x, p, n);
- }
- else if (abstop12 (y) < abstop12 (INFINITY))
- {
- uint32_t xi = asuint (y);
- int sign = xi >> 31;
-
- x = reduce_large (xi, &n);
-
- /* Setup signs for sin and cos - include original sign. */
- s = p->sign[(n + sign) & 3];
-
- if ((n + sign) & 2)
- p = &__sincosf_table[1];
-
- return sinf_poly (x * s, x * x, p, n);
- }
- else
- return __math_invalidf (y);
-}
diff --git a/libc/AOR_v20.02/math/test/testcases/directed/cosf.tst b/libc/AOR_v20.02/math/test/testcases/directed/cosf.tst
deleted file mode 100644
index 8c7621a4550c..000000000000
--- a/libc/AOR_v20.02/math/test/testcases/directed/cosf.tst
+++ /dev/null
@@ -1,26 +0,0 @@
-; cosf.tst - Directed test cases for SP cosine
-;
-; 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
-
-func=cosf op1=7fc00001 result=7fc00001 errno=0
-func=cosf op1=ffc00001 result=7fc00001 errno=0
-func=cosf op1=7f800001 result=7fc00001 errno=0 status=i
-func=cosf op1=ff800001 result=7fc00001 errno=0 status=i
-func=cosf op1=7f800000 result=7fc00001 errno=EDOM status=i
-func=cosf op1=ff800000 result=7fc00001 errno=EDOM status=i
-func=cosf op1=00000000 result=3f800000 errno=0
-func=cosf op1=80000000 result=3f800000 errno=0
-; SDCOMP-26094: check cosf in the cases for which the range reducer
-; returns values furthest beyond its nominal upper bound of pi/4.
-func=cosf op1=46427f1b result=3f34dc5c.565 error=0
-func=cosf op1=4647e568 result=3f34dc33.c1f error=0
-func=cosf op1=46428bac result=bf34dbf2.8e3 error=0
-func=cosf op1=4647f1f9 result=bf34dbc9.f9b error=0
-func=cosf op1=4647fe8a result=3f34db60.313 error=0
-func=cosf op1=45d8d7f1 result=bf35006a.7fd error=0
-func=cosf op1=45d371a4 result=3f350056.39b error=0
-func=cosf op1=45ce0b57 result=bf350041.f38 error=0
-func=cosf op1=45d35882 result=bf34ffec.868 error=0
-func=cosf op1=45cdf235 result=3f34ffd8.404 error=0
diff --git a/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst b/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst
deleted file mode 100644
index d22fd9802694..000000000000
--- a/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst
+++ /dev/null
@@ -1,52 +0,0 @@
-; Directed test cases for SP sincos
-;
-; 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
-
-
-func=sincosf_sinf op1=7fc00001 result=7fc00001 errno=0
-func=sincosf_sinf op1=ffc00001 result=7fc00001 errno=0
-func=sincosf_sinf op1=7f800001 result=7fc00001 errno=0 status=i
-func=sincosf_sinf op1=ff800001 result=7fc00001 errno=0 status=i
-func=sincosf_sinf op1=7f800000 result=7fc00001 errno=EDOM status=i
-func=sincosf_sinf op1=ff800000 result=7fc00001 errno=EDOM status=i
-func=sincosf_sinf op1=00000000 result=00000000 errno=0
-func=sincosf_sinf op1=80000000 result=80000000 errno=0
-func=sincosf_sinf op1=c70d39a1 result=be37fad5.7ed errno=0
-func=sincosf_sinf op1=46427f1b result=3f352d80.f9b error=0
-func=sincosf_sinf op1=4647e568 result=3f352da9.7be error=0
-func=sincosf_sinf op1=46428bac result=bf352dea.924 error=0
-func=sincosf_sinf op1=4647f1f9 result=bf352e13.146 error=0
-func=sincosf_sinf op1=4647fe8a result=3f352e7c.ac9 error=0
-func=sincosf_sinf op1=45d8d7f1 result=3f35097b.cb0 error=0
-func=sincosf_sinf op1=45d371a4 result=bf350990.102 error=0
-func=sincosf_sinf op1=45ce0b57 result=3f3509a4.554 error=0
-func=sincosf_sinf op1=45d35882 result=3f3509f9.bdb error=0
-func=sincosf_sinf op1=45cdf235 result=bf350a0e.02c error=0
-
-func=sincosf_cosf op1=7fc00001 result=7fc00001 errno=0
-func=sincosf_cosf op1=ffc00001 result=7fc00001 errno=0
-func=sincosf_cosf op1=7f800001 result=7fc00001 errno=0 status=i
-func=sincosf_cosf op1=ff800001 result=7fc00001 errno=0 status=i
-func=sincosf_cosf op1=7f800000 result=7fc00001 errno=EDOM status=i
-func=sincosf_cosf op1=ff800000 result=7fc00001 errno=EDOM status=i
-func=sincosf_cosf op1=00000000 result=3f800000 errno=0
-func=sincosf_cosf op1=80000000 result=3f800000 errno=0
-func=sincosf_cosf op1=46427f1b result=3f34dc5c.565 error=0
-func=sincosf_cosf op1=4647e568 result=3f34dc33.c1f error=0
-func=sincosf_cosf op1=46428bac result=bf34dbf2.8e3 error=0
-func=sincosf_cosf op1=4647f1f9 result=bf34dbc9.f9b error=0
-func=sincosf_cosf op1=4647fe8a result=3f34db60.313 error=0
-func=sincosf_cosf op1=45d8d7f1 result=bf35006a.7fd error=0
-func=sincosf_cosf op1=45d371a4 result=3f350056.39b error=0
-func=sincosf_cosf op1=45ce0b57 result=bf350041.f38 error=0
-func=sincosf_cosf op1=45d35882 result=bf34ffec.868 error=0
-func=sincosf_cosf op1=45cdf235 result=3f34ffd8.404 error=0
-
-; no underflow
-func=sincosf_sinf op1=17800000 result=17800000.000
-func=sincosf_cosf op1=17800000 result=3f800000.000
-; underflow
-func=sincosf_sinf op1=00400000 result=00400000.000 status=ux
-func=sincosf_cosf op1=00400000 result=3f800000.000 status=ux
diff --git a/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst b/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst
deleted file mode 100644
index 022bf1424879..000000000000
--- a/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst
+++ /dev/null
@@ -1,29 +0,0 @@
-; sinf.tst - Directed test cases for SP sine
-;
-; 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
-
-
-func=sinf op1=7fc00001 result=7fc00001 errno=0
-func=sinf op1=ffc00001 result=7fc00001 errno=0
-func=sinf op1=7f800001 result=7fc00001 errno=0 status=i
-func=sinf op1=ff800001 result=7fc00001 errno=0 status=i
-func=sinf op1=7f800000 result=7fc00001 errno=EDOM status=i
-func=sinf op1=ff800000 result=7fc00001 errno=EDOM status=i
-func=sinf op1=00000000 result=00000000 errno=0
-func=sinf op1=80000000 result=80000000 errno=0
-; Directed test for a failure I found while developing mathbench
-func=sinf op1=c70d39a1 result=be37fad5.7ed errno=0
-; SDCOMP-26094: check sinf in the cases for which the range reducer
-; returns values furthest beyond its nominal upper bound of pi/4.
-func=sinf op1=46427f1b result=3f352d80.f9b error=0
-func=sinf op1=4647e568 result=3f352da9.7be error=0
-func=sinf op1=46428bac result=bf352dea.924 error=0
-func=sinf op1=4647f1f9 result=bf352e13.146 error=0
-func=sinf op1=4647fe8a result=3f352e7c.ac9 error=0
-func=sinf op1=45d8d7f1 result=3f35097b.cb0 error=0
-func=sinf op1=45d371a4 result=bf350990.102 error=0
-func=sinf op1=45ce0b57 result=3f3509a4.554 error=0
-func=sinf op1=45d35882 result=3f3509f9.bdb error=0
-func=sinf op1=45cdf235 result=bf350a0e.02c error=0
diff --git a/libc/AOR_v20.02/math/test/testcases/random/float.tst b/libc/AOR_v20.02/math/test/testcases/random/float.tst
index c142d63cd594..aadd70a336a1 100644
--- a/libc/AOR_v20.02/math/test/testcases/random/float.tst
+++ b/libc/AOR_v20.02/math/test/testcases/random/float.tst
@@ -4,10 +4,6 @@
!! See https://llvm.org/LICENSE.txt for license information.
!! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-test sinf 10000
-test cosf 10000
-test sincosf_sinf 5000
-test sincosf_cosf 5000
test tanf 10000
test expf 10000
test exp2f 10000
diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td
index fe3064ddf227..f176caee6a4e 100644
--- a/libc/config/linux/api.td
+++ b/libc/config/linux/api.td
@@ -1,5 +1,6 @@
include "config/public_api.td"
+include "spec/gnu_ext.td"
include "spec/linux.td"
include "spec/posix.td"
include "spec/stdc.td"
@@ -123,7 +124,10 @@ def MathAPI : PublicAPI<"math.h"> {
IsNanMacro,
];
let Functions = [
+ "cosf",
"round",
+ "sincosf",
+ "sinf",
];
}
diff --git a/libc/lib/CMakeLists.txt b/libc/lib/CMakeLists.txt
index 79456c938e6e..bc245cdb481e 100644
--- a/libc/lib/CMakeLists.txt
+++ b/libc/lib/CMakeLists.txt
@@ -41,7 +41,10 @@ add_entrypoint_library(
llvmlibm
DEPENDS
# math.h entrypoints
+ libc.src.math.cosf
libc.src.math.round
+ libc.src.math.sincosf
+ libc.src.math.sinf
)
add_redirector_library(
diff --git a/libc/src/__support/common.h.def b/libc/src/__support/common.h.def
index 3abd23674c97..a1bb78d5b00f 100644
--- a/libc/src/__support/common.h.def
+++ b/libc/src/__support/common.h.def
@@ -11,6 +11,10 @@
#define LIBC_INLINE_ASM __asm__ __volatile__
+#define likely(x) __builtin_expect (!!(x), 1)
+#define unlikely(x) __builtin_expect (x, 0)
+#define UNUSED __attribute__((unused))
+
<!> Include the platform specific definitions at build time. For example, that
<!> of entrypoint macro.
%%include_file(${platform_defs})
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index e7ad9d2a73fe..5ec98ddd18dd 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -1,3 +1,23 @@
+add_header_library(
+ math_utils
+ HDRS
+ math_utils.h
+ DEPENDS
+ libc.include.errno
+ libc.include.math
+ libc.src.errno.__errno_location
+)
+
+add_object_library(
+ sincosf_utils
+ HDRS
+ sincosf_utils.h
+ SRCS
+ sincosf_data.cpp
+ DEPENDS
+ .math_utils
+)
+
add_entrypoint_object(
round
REDIRECTED
@@ -12,3 +32,39 @@ add_redirector_object(
SRC
round_redirector.cpp
)
+
+add_entrypoint_object(
+ cosf
+ SRCS
+ cosf.cpp
+ HDRS
+ cosf.h
+ DEPENDS
+ .sincosf_utils
+ libc.include.math
+ libc.src.errno.__errno_location
+)
+
+add_entrypoint_object(
+ sinf
+ SRCS
+ sinf.cpp
+ HDRS
+ sinf.h
+ DEPENDS
+ .sincosf_utils
+ libc.include.math
+ libc.src.errno.__errno_location
+)
+
+add_entrypoint_object(
+ sincosf
+ SRCS
+ sincosf.cpp
+ HDRS
+ sincosf.h
+ DEPENDS
+ .sincosf_utils
+ libc.include.math
+ libc.src.errno.__errno_location
+)
diff --git a/libc/src/math/cosf.cpp b/libc/src/math/cosf.cpp
new file mode 100644
index 000000000000..db121b2cb396
--- /dev/null
+++ b/libc/src/math/cosf.cpp
@@ -0,0 +1,64 @@
+//===-- Single-precision cos 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 "math_utils.h"
+#include "sincosf_utils.h"
+
+#include "include/math.h"
+#include "src/__support/common.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+// Fast cosf implementation. Worst-case ULP is 0.5607, maximum relative
+// error is 0.5303 * 2^-23. A single-step range reduction is used for
+// small values. Large inputs have their range reduced using fast integer
+// arithmetic.
+float LLVM_LIBC_ENTRYPOINT(cosf)(float y) {
+ double x = y;
+ double s;
+ int n;
+ const sincos_t *p = &__sincosf_table[0];
+
+ if (abstop12(y) < abstop12(pio4)) {
+ double x2 = x * x;
+
+ if (unlikely(abstop12(y) < abstop12(as_float(0x39800000))))
+ return 1.0f;
+
+ return sinf_poly(x, x2, p, 1);
+ } else if (likely(abstop12(y) < abstop12(120.0f))) {
+ x = reduce_fast(x, p, &n);
+
+ // Setup the signs for sin and cos.
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &__sincosf_table[1];
+
+ return sinf_poly(x * s, x * x, p, n ^ 1);
+ } else if (abstop12(y) < abstop12(INFINITY)) {
+ uint32_t xi = as_uint32_bits(y);
+ int sign = xi >> 31;
+
+ x = reduce_large(xi, &n);
+
+ // Setup signs for sin and cos - include original sign.
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &__sincosf_table[1];
+
+ return sinf_poly(x * s, x * x, p, n ^ 1);
+ }
+
+ return invalidf(y);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/cosf.h b/libc/src/math/cosf.h
new file mode 100644
index 000000000000..1aaabe900ba8
--- /dev/null
+++ b/libc/src/math/cosf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for cosf --------------------------*- 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_COSF_H
+#define LLVM_LIBC_SRC_MATH_COSF_H
+
+namespace __llvm_libc {
+
+float cosf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_COSF_H
diff --git a/libc/src/math/math_utils.h b/libc/src/math/math_utils.h
new file mode 100644
index 000000000000..3553673486f9
--- /dev/null
+++ b/libc/src/math/math_utils.h
@@ -0,0 +1,49 @@
+//===-- Collection of utils for implementing math functions -----*- 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_MATH_UTILS_H
+#define LLVM_LIBC_SRC_MATH_MATH_UTILS_H
+
+#include "include/errno.h"
+#include "include/math.h"
+
+#include "src/__support/common.h"
+#include "src/errno/llvmlibc_errno.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+static inline float with_errnof(float x, int err) {
+ if (math_errhandling & MATH_ERRNO)
+ llvmlibc_errno = err;
+ return x;
+}
+
+static inline uint32_t as_uint32_bits(float x) {
+ return *reinterpret_cast<uint32_t *>(&x);
+}
+
+static inline float as_float(uint32_t x) {
+ return *reinterpret_cast<float *>(&x);
+}
+
+static inline double as_double(uint64_t x) {
+ return *reinterpret_cast<double *>(&x);
+}
+
+static inline constexpr float invalidf(float x) {
+ float y = (x - x) / (x - x);
+ return isnan(x) ? y : with_errnof(y, EDOM);
+}
+
+static inline void force_eval_float(float x) { volatile float y UNUSED = x; }
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_MATH_UTILS_H
diff --git a/libc/src/math/sincosf.cpp b/libc/src/math/sincosf.cpp
new file mode 100644
index 000000000000..717feaa470d2
--- /dev/null
+++ b/libc/src/math/sincosf.cpp
@@ -0,0 +1,76 @@
+//===-- Single-precision sincos 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 "math_utils.h"
+#include "sincosf_utils.h"
+
+#include "include/math.h"
+#include "src/__support/common.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+// Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative
+// error is 0.5303 * 2^-23. A single-step range reduction is used for
+// small values. Large inputs have their range reduced using fast integer
+// arithmetic.
+void LLVM_LIBC_ENTRYPOINT(sincosf)(float y, float *sinp, float *cosp) {
+ double x = y;
+ double s;
+ int n;
+ const sincos_t *p = &__sincosf_table[0];
+
+ if (abstop12(y) < abstop12(pio4)) {
+ double x2 = x * x;
+
+ if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) {
+ if (unlikely(abstop12(y) < abstop12(as_float(0x800000))))
+ // Force underflow for tiny y.
+ force_eval_float(x2);
+ *sinp = y;
+ *cosp = 1.0f;
+ return;
+ }
+
+ sincosf_poly(x, x2, p, 0, sinp, cosp);
+ } else if (abstop12(y) < abstop12(120.0f)) {
+ x = reduce_fast(x, p, &n);
+
+ // Setup the signs for sin and cos.
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &__sincosf_table[1];
+
+ sincosf_poly(x * s, x * x, p, n, sinp, cosp);
+ } else if (likely(abstop12(y) < abstop12(INFINITY))) {
+ uint32_t xi = as_uint32_bits(y);
+ int sign = xi >> 31;
+
+ x = reduce_large(xi, &n);
+
+ // Setup signs for sin and cos - include original sign.
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &__sincosf_table[1];
+
+ sincosf_poly(x * s, x * x, p, n, sinp, cosp);
+ } else {
+ // Return NaN if Inf or NaN for both sin and cos.
+ *sinp = *cosp = y - y;
+
+ // Needed to set errno for +-Inf, the add is a hack to work
+ // around a gcc register allocation issue: just passing y
+ // affects code generation in the fast path.
+ invalidf(y + y);
+ }
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/sincosf.h b/libc/src/math/sincosf.h
new file mode 100644
index 000000000000..47ef983f4385
--- /dev/null
+++ b/libc/src/math/sincosf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for sincosf -----------------------*- 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_SINCOSF_H
+#define LLVM_LIBC_SRC_MATH_SINCOSF_H
+
+namespace __llvm_libc {
+
+void sincosf(float x, float *sinx, float *cosx);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_SINCOSF_H
diff --git a/libc/src/math/sincosf_data.cpp b/libc/src/math/sincosf_data.cpp
new file mode 100644
index 000000000000..50984570b55f
--- /dev/null
+++ b/libc/src/math/sincosf_data.cpp
@@ -0,0 +1,51 @@
+//===-- sinf/cosf data tables ---------------------------------------------===//
+//
+// 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 "math_utils.h"
+#include "sincosf_utils.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+// The constants and polynomials for sine and cosine. The 2nd entry
+// computes -cos (x) rather than cos (x) to get negation for free.
+const sincos_t __sincosf_table[2] = {
+ {{1.0, -1.0, -1.0, 1.0},
+ as_double(0x41645f306dc9c883),
+ as_double(0x3ff921fb54442d18),
+ as_double(0x3ff0000000000000),
+ as_double(0xbfdffffffd0c621c),
+ as_double(0x3fa55553e1068f19),
+ as_double(0xbf56c087e89a359d),
+ as_double(0x3ef99343027bf8c3),
+ as_double(0xbfc555545995a603),
+ as_double(0x3f81107605230bc4),
+ as_double(0xbf2994eb3774cf24)},
+ {{1.0, -1.0, -1.0, 1.0},
+ as_double(0x41645f306dc9c883),
+ as_double(0x3ff921fb54442d18),
+ as_double(0xbff0000000000000),
+ as_double(0x3fdffffffd0c621c),
+ as_double(0xbfa55553e1068f19),
+ as_double(0x3f56c087e89a359d),
+ as_double(0xbef99343027bf8c3),
+ as_double(0xbfc555545995a603),
+ as_double(0x3f81107605230bc4),
+ as_double(0xbf2994eb3774cf24)},
+};
+
+// Table with 4/PI to 192 bit precision. To avoid unaligned accesses
+// only 8 new bits are added per entry, making the table 4 times larger.
+const uint32_t __inv_pio4[24] = {
+ 0xa2, 0xa2f9, 0xa2f983, 0xa2f9836e, 0xf9836e4e, 0x836e4e44,
+ 0x6e4e4415, 0x4e441529, 0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1,
+ 0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0, 0x34ddc0db, 0xddc0db62,
+ 0xc0db6295, 0xdb629599, 0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041};
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/sincosf_utils.h b/libc/src/math/sincosf_utils.h
new file mode 100644
index 000000000000..8c54cb9c1d90
--- /dev/null
+++ b/libc/src/math/sincosf_utils.h
@@ -0,0 +1,142 @@
+//===-- Collection of utils for cosf/sinf/sincosf ---------------*- 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_SINCOSF_UTILS_H
+#define LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H
+
+#include "math_utils.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+// 2PI * 2^-64.
+static const double pi63 = as_double(0x3c1921fb54442d18);
+// PI / 4.
+static const double pio4 = as_double(0x3fe921fb54442d18);
+
+// The constants and polynomials for sine and cosine.
+typedef struct {
+ double sign[4]; // Sign of sine in quadrants 0..3.
+ double hpi_inv; // 2 / PI ( * 2^24 ).
+ double hpi; // PI / 2.
+ double c0, c1, c2, c3, c4; // Cosine polynomial.
+ double s1, s2, s3; // Sine polynomial.
+} sincos_t;
+
+// Polynomial data (the cosine polynomial is negated in the 2nd entry).
+extern const sincos_t __sincosf_table[2];
+
+// Table with 4/PI to 192 bit precision.
+extern const uint32_t __inv_pio4[];
+
+// Top 12 bits of the float representation with the sign bit cleared.
+static inline uint32_t abstop12(float x) {
+ return (as_uint32_bits(x) >> 20) & 0x7ff;
+}
+
+// Compute the sine and cosine of inputs X and X2 (X squared), using the
+// polynomial P and store the results in SINP and COSP. N is the quadrant,
+// if odd the cosine and sine polynomials are swapped.
+static inline void sincosf_poly(double x, double x2, const sincos_t *p, int n,
+ float *sinp, float *cosp) {
+ double x3, x4, x5, x6, s, c, c1, c2, s1;
+
+ x4 = x2 * x2;
+ x3 = x2 * x;
+ c2 = p->c3 + x2 * p->c4;
+ s1 = p->s2 + x2 * p->s3;
+
+ // Swap sin/cos result based on quadrant.
+ float *tmp = (n & 1 ? cosp : sinp);
+ cosp = (n & 1 ? sinp : cosp);
+ sinp = tmp;
+
+ c1 = p->c0 + x2 * p->c1;
+ x5 = x3 * x2;
+ x6 = x4 * x2;
+
+ s = x + x3 * p->s1;
+ c = c1 + x4 * p->c2;
+
+ *sinp = s + x5 * s1;
+ *cosp = c + x6 * c2;
+}
+
+// Return the sine of inputs X and X2 (X squared) using the polynomial P.
+// N is the quadrant, and if odd the cosine polynomial is used.
+static inline float sinf_poly(double x, double x2, const sincos_t *p, int n) {
+ double x3, x4, x6, x7, s, c, c1, c2, s1;
+
+ if ((n & 1) == 0) {
+ x3 = x * x2;
+ s1 = p->s2 + x2 * p->s3;
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1;
+
+ return s + x7 * s1;
+ } else {
+ x4 = x2 * x2;
+ c2 = p->c3 + x2 * p->c4;
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->c2;
+
+ return c + x6 * c2;
+ }
+}
+
+// Fast range reduction using single multiply-subtract. Return the modulo of
+// X as a value between -PI/4 and PI/4 and store the quadrant in NP.
+// The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
+// is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
+// the result is accurate for |X| <= 120.0.
+static inline double reduce_fast(double x, const sincos_t *p, int *np) {
+ double r;
+ // Use scaled float to int conversion with explicit rounding.
+ // hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
+ // This avoids inaccuracies introduced by truncating negative values.
+ r = x * p->hpi_inv;
+ int n = ((int32_t)r + 0x800000) >> 24;
+ *np = n;
+ return x - n * p->hpi;
+}
+
+// Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
+// XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
+// Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
+// Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit
+// multiply computes the exact 2.62-bit fixed-point modulo. Since the result
+// can have at most 29 leading zeros after the binary point, the double
+// precision result is accurate to 33 bits.
+static inline double reduce_large(uint32_t xi, int *np) {
+ const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
+ int shift = (xi >> 23) & 7;
+ uint64_t n, res0, res1, res2;
+
+ xi = (xi & 0xffffff) | 0x800000;
+ xi <<= shift;
+
+ res0 = xi * arr[0];
+ res1 = (uint64_t)xi * arr[4];
+ res2 = (uint64_t)xi * arr[8];
+ res0 = (res2 >> 32) | (res0 << 32);
+ res0 += res1;
+
+ n = (res0 + (1ULL << 61)) >> 62;
+ res0 -= n << 62;
+ double x = (int64_t)res0;
+ *np = n;
+ return x * pi63;
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H
diff --git a/libc/src/math/sinf.cpp b/libc/src/math/sinf.cpp
new file mode 100644
index 000000000000..634e40c50f6b
--- /dev/null
+++ b/libc/src/math/sinf.cpp
@@ -0,0 +1,68 @@
+//===-- Single-precision sin 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 "math_utils.h"
+#include "sincosf_utils.h"
+
+#include "include/math.h"
+#include "src/__support/common.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+
+// Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative
+// error is 0.5303 * 2^-23. A single-step range reduction is used for
+// small values. Large inputs have their range reduced using fast integer
+// arithmetic.
+float LLVM_LIBC_ENTRYPOINT(sinf)(float y) {
+ double x = y;
+ double s;
+ int n;
+ const sincos_t *p = &__sincosf_table[0];
+
+ if (abstop12(y) < abstop12(pio4)) {
+ s = x * x;
+
+ if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) {
+ if (unlikely(abstop12(y) < abstop12(as_float(0x800000))))
+ // Force underflow for tiny y.
+ force_eval_float(s);
+ return y;
+ }
+
+ return sinf_poly(x, s, p, 0);
+ } else if (likely(abstop12(y) < abstop12(120.0f))) {
+ x = reduce_fast(x, p, &n);
+
+ // Setup the signs for sin and cos.
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &__sincosf_table[1];
+
+ return sinf_poly(x * s, x * x, p, n);
+ } else if (abstop12(y) < abstop12(INFINITY)) {
+ uint32_t xi = as_uint32_bits(y);
+ int sign = xi >> 31;
+
+ x = reduce_large(xi, &n);
+
+ // Setup signs for sin and cos - include original sign.
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &__sincosf_table[1];
+
+ return sinf_poly(x * s, x * x, p, n);
+ }
+
+ return invalidf(y);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/sinf.h b/libc/src/math/sinf.h
new file mode 100644
index 000000000000..e63db04c51b5
--- /dev/null
+++ b/libc/src/math/sinf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for sinf --------------------------*- 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_SINF_H
+#define LLVM_LIBC_SRC_MATH_SINF_H
+
+namespace __llvm_libc {
+
+float sinf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_SINF_H
diff --git a/libc/test/src/CMakeLists.txt b/libc/test/src/CMakeLists.txt
index 209d00b6d7f9..f92a0463b359 100644
--- a/libc/test/src/CMakeLists.txt
+++ b/libc/test/src/CMakeLists.txt
@@ -1,5 +1,6 @@
add_subdirectory(assert)
add_subdirectory(errno)
+add_subdirectory(math)
add_subdirectory(signal)
add_subdirectory(stdio)
add_subdirectory(stdlib)
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
new file mode 100644
index 000000000000..4ba3ff3f5d3a
--- /dev/null
+++ b/libc/test/src/math/CMakeLists.txt
@@ -0,0 +1,80 @@
+add_libc_testsuite(libc_math_unittests)
+
+function(add_math_unittest name)
+ cmake_parse_arguments(
+ "MATH_UNITTEST"
+ "NEED_MPFR" # No optional arguments
+ "" # Single value arguments
+ "" # Multi-value arguments
+ ${ARGN}
+ )
+
+ if(MATH_UNITTEST_NEED_MPFR)
+ if(NOT LIBC_TESTS_CAN_USE_MPFR)
+ message("WARNING: Math test ${name} will be skipped as MPFR library is not available.")
+ return()
+ endif()
+ endif()
+
+ add_libc_unittest(${name} ${MATH_UNITTEST_UNPARSED_ARGUMENTS})
+ if(MATH_UNITTEST_NEED_MPFR)
+ get_fq_target_name(${name} fq_target_name)
+ target_link_libraries(${fq_target_name} PRIVATE libcMPFRWrapper -lmpfr -lgmp)
+ endif()
+endfunction(add_math_unittest)
+
+add_header_library(
+ float_utils
+ HDRS
+ float.h
+)
+
+# TODO(sivachandra): Remove the dependency on __errno_location as the tested
+# entry points depend on the already.
+add_math_unittest(
+ cosf_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ cosf_test.cpp
+ HDRS
+ sdcomp26094.h
+ DEPENDS
+ .float_utils
+ libc.src.errno.__errno_location
+ libc.src.math.cosf
+ libc.utils.CPP.standalone_cpp
+)
+
+add_math_unittest(
+ sinf_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ sinf_test.cpp
+ HDRS
+ sdcomp26094.h
+ DEPENDS
+ .float_utils
+ libc.src.errno.__errno_location
+ libc.src.math.sinf
+ libc.utils.CPP.standalone_cpp
+)
+
+add_math_unittest(
+ sincosf_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ sincosf_test.cpp
+ HDRS
+ sdcomp26094.h
+ DEPENDS
+ .float_utils
+ libc.src.errno.__errno_location
+ libc.src.math.sincosf
+ libc.utils.CPP.standalone_cpp
+)
diff --git a/libc/test/src/math/cosf_test.cpp b/libc/test/src/math/cosf_test.cpp
new file mode 100644
index 000000000000..94c66cda1b0f
--- /dev/null
+++ b/libc/test/src/math/cosf_test.cpp
@@ -0,0 +1,103 @@
+//===-- Unittests for cosf ------------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/errno/llvmlibc_errno.h"
+#include "src/math/cosf.h"
+#include "src/math/math_utils.h"
+#include "test/src/math/float.h"
+#include "test/src/math/sdcomp26094.h"
+#include "utils/CPP/Array.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <stdint.h>
+
+using __llvm_libc::as_float;
+using __llvm_libc::as_uint32_bits;
+
+using __llvm_libc::testing::FloatBits;
+using __llvm_libc::testing::sdcomp26094Values;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+// 12 additional bits of precision over the base precision of a |float|
+// value.
+static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12,
+ 3 * 0x1000 / 4};
+
+TEST(CosfTest, SpecialNumbers) {
+ llvmlibc_errno = 0;
+
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::QNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegQNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::SNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegSNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_EQ(FloatBits::One,
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::Zero))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_EQ(FloatBits::One,
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegZero))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ llvmlibc_errno = 0;
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::Inf)))));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+
+ llvmlibc_errno = 0;
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegInf)))));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+}
+
+TEST(CosfTest, 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 = as_float(v);
+ if (isnan(x) || isinf(x))
+ continue;
+ EXPECT_TRUE(mpfr::equalsCos(x, __llvm_libc::cosf(x), tolerance));
+ }
+}
+
+// For small values, cos(x) is 1.
+TEST(CosfTest, SmallValues) {
+ float x = as_float(0x17800000);
+ float result = __llvm_libc::cosf(x);
+ EXPECT_TRUE(mpfr::equalsCos(x, result, tolerance));
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(result));
+
+ x = as_float(0x00400000);
+ result = __llvm_libc::cosf(x);
+ EXPECT_TRUE(mpfr::equalsCos(x, result, tolerance));
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(result));
+}
+
+// SDCOMP-26094: check cosf in the cases for which the range reducer
+// returns values furthest beyond its nominal upper bound of pi/4.
+TEST(CosfTest, SDCOMP_26094) {
+ for (uint32_t v : sdcomp26094Values) {
+ float x = as_float(v);
+ EXPECT_TRUE(mpfr::equalsCos(x, __llvm_libc::cosf(x), tolerance));
+ }
+}
diff --git a/libc/test/src/math/float.h b/libc/test/src/math/float.h
new file mode 100644
index 000000000000..bfa15a18ce3f
--- /dev/null
+++ b/libc/test/src/math/float.h
@@ -0,0 +1,49 @@
+//===-- Single precision floating point test utils --------------*- 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_TEST_SRC_MATH_FLOAT_H
+#define LLVM_LIBC_TEST_SRC_MATH_FLOAT_H
+
+#include "src/math/math_utils.h"
+
+namespace __llvm_libc {
+namespace testing {
+
+struct FloatBits {
+ // The various NaN bit patterns here are just one of the many possible
+ // patterns. The functions isQNan and isNegQNan can help understand why.
+
+ static const uint32_t QNan = 0x7fc00000;
+ static const uint32_t NegQNan = 0xffc00000;
+
+ static const uint32_t SNan = 0x7f800001;
+ static const uint32_t NegSNan = 0xff800001;
+
+ static bool isQNan(float f) {
+ uint32_t bits = as_uint32_bits(f);
+ return ((0x7fc00000 & bits) != 0) && ((0x80000000 & bits) == 0);
+ }
+
+ static bool isNegQNan(float f) {
+ uint32_t bits = as_uint32_bits(f);
+ return 0xffc00000 & bits;
+ }
+
+ static constexpr uint32_t Zero = 0x0;
+ static constexpr uint32_t NegZero = 0x80000000;
+
+ static constexpr uint32_t Inf = 0x7f800000;
+ static constexpr uint32_t NegInf = 0xff800000;
+
+ static constexpr uint32_t One = 0x3f800000;
+};
+
+} // namespace testing
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_FLOAT_H
diff --git a/libc/test/src/math/sdcomp26094.h b/libc/test/src/math/sdcomp26094.h
new file mode 100644
index 000000000000..7ebcd60e4af2
--- /dev/null
+++ b/libc/test/src/math/sdcomp26094.h
@@ -0,0 +1,25 @@
+//===-- SDCOMP-26094 specific items -----------------------------*- 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_TEST_SRC_MATH_SDCOMP26094_H
+#define LLVM_LIBC_TEST_SRC_MATH_SDCOMP26094_H
+
+#include "utils/CPP/Array.h"
+
+namespace __llvm_libc {
+namespace testing {
+
+static constexpr __llvm_libc::cpp::Array<uint32_t, 10> sdcomp26094Values{
+ 0x46427f1b, 0x4647e568, 0x46428bac, 0x4647f1f9, 0x4647fe8a,
+ 0x45d8d7f1, 0x45d371a4, 0x45ce0b57, 0x45d35882, 0x45cdf235,
+};
+
+} // namespace testing
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SDCOMP26094_H
diff --git a/libc/test/src/math/sincosf_test.cpp b/libc/test/src/math/sincosf_test.cpp
new file mode 100644
index 000000000000..36e6b4a129a7
--- /dev/null
+++ b/libc/test/src/math/sincosf_test.cpp
@@ -0,0 +1,125 @@
+//===-- Unittests for sincosf ---------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/errno/llvmlibc_errno.h"
+#include "src/math/math_utils.h"
+#include "src/math/sincosf.h"
+#include "test/src/math/float.h"
+#include "test/src/math/sdcomp26094.h"
+#include "utils/CPP/Array.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <stdint.h>
+
+using __llvm_libc::as_float;
+using __llvm_libc::as_uint32_bits;
+
+using __llvm_libc::testing::FloatBits;
+using __llvm_libc::testing::sdcomp26094Values;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+// 12 additional bits of precision over the base precision of a |float|
+// value.
+static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12,
+ 3 * 0x1000 / 4};
+
+TEST(SinCosfTest, SpecialNumbers) {
+ llvmlibc_errno = 0;
+ float sin, cos;
+
+ __llvm_libc::sincosf(as_float(FloatBits::QNan), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ __llvm_libc::sincosf(as_float(FloatBits::NegQNan), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ __llvm_libc::sincosf(as_float(FloatBits::SNan), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ __llvm_libc::sincosf(as_float(FloatBits::NegSNan), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ __llvm_libc::sincosf(as_float(FloatBits::Zero), &sin, &cos);
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(cos));
+ EXPECT_EQ(FloatBits::Zero, as_uint32_bits(sin));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ __llvm_libc::sincosf(as_float(FloatBits::NegZero), &sin, &cos);
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(cos));
+ EXPECT_EQ(FloatBits::NegZero, as_uint32_bits(sin));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ llvmlibc_errno = 0;
+ __llvm_libc::sincosf(as_float(FloatBits::Inf), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+
+ llvmlibc_errno = 0;
+ __llvm_libc::sincosf(as_float(FloatBits::NegInf), &sin, &cos);
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos)));
+ EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin)));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+}
+
+TEST(SinCosfTest, 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 = as_float(v);
+ if (isnan(x) || isinf(x))
+ continue;
+
+ float sin, cos;
+ __llvm_libc::sincosf(x, &sin, &cos);
+ EXPECT_TRUE(mpfr::equalsCos(x, cos, tolerance));
+ EXPECT_TRUE(mpfr::equalsSin(x, sin, tolerance));
+ }
+}
+
+// For small values, cos(x) is 1 and sin(x) is x.
+TEST(SinCosfTest, SmallValues) {
+ uint32_t bits = 0x17800000;
+ float x = as_float(bits);
+ float result_cos, result_sin;
+ __llvm_libc::sincosf(x, &result_sin, &result_cos);
+ EXPECT_TRUE(mpfr::equalsCos(x, result_cos, tolerance));
+ EXPECT_TRUE(mpfr::equalsSin(x, result_sin, tolerance));
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(result_cos));
+ EXPECT_EQ(bits, as_uint32_bits(result_sin));
+
+ bits = 0x00400000;
+ x = as_float(bits);
+ __llvm_libc::sincosf(x, &result_sin, &result_cos);
+ EXPECT_TRUE(mpfr::equalsCos(x, result_cos, tolerance));
+ EXPECT_TRUE(mpfr::equalsSin(x, result_sin, tolerance));
+ EXPECT_EQ(FloatBits::One, as_uint32_bits(result_cos));
+ EXPECT_EQ(bits, as_uint32_bits(result_sin));
+}
+
+// SDCOMP-26094: check sinf in the cases for which the range reducer
+// returns values furthest beyond its nominal upper bound of pi/4.
+TEST(SinCosfTest, SDCOMP_26094) {
+ for (uint32_t v : sdcomp26094Values) {
+ float x = as_float(v);
+ float sin, cos;
+ __llvm_libc::sincosf(x, &sin, &cos);
+ EXPECT_TRUE(mpfr::equalsCos(x, cos, tolerance));
+ EXPECT_TRUE(mpfr::equalsSin(x, sin, tolerance));
+ }
+}
diff --git a/libc/test/src/math/sinf_test.cpp b/libc/test/src/math/sinf_test.cpp
new file mode 100644
index 000000000000..e4c6e818b57a
--- /dev/null
+++ b/libc/test/src/math/sinf_test.cpp
@@ -0,0 +1,110 @@
+//===-- Unittests for sinf ------------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/errno/llvmlibc_errno.h"
+#include "src/math/math_utils.h"
+#include "src/math/sinf.h"
+#include "test/src/math/float.h"
+#include "test/src/math/sdcomp26094.h"
+#include "utils/CPP/Array.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <stdint.h>
+
+using __llvm_libc::as_float;
+using __llvm_libc::as_uint32_bits;
+
+using __llvm_libc::testing::FloatBits;
+using __llvm_libc::testing::sdcomp26094Values;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+// 12 additional bits of precision over the base precision of a |float|
+// value.
+static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12,
+ 3 * 0x1000 / 4};
+
+TEST(SinfTest, SpecialNumbers) {
+ llvmlibc_errno = 0;
+
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::QNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegQNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::SNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegSNan)))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_EQ(FloatBits::Zero,
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::Zero))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ EXPECT_EQ(FloatBits::NegZero,
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegZero))));
+ EXPECT_EQ(llvmlibc_errno, 0);
+
+ llvmlibc_errno = 0;
+ EXPECT_TRUE(FloatBits::isQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::Inf)))));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+
+ llvmlibc_errno = 0;
+ EXPECT_TRUE(FloatBits::isNegQNan(
+ as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegInf)))));
+ EXPECT_EQ(llvmlibc_errno, EDOM);
+}
+
+TEST(SinfTest, 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 = as_float(v);
+ if (isnan(x) || isinf(x))
+ continue;
+ EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance));
+ }
+}
+
+TEST(SinfTest, SpecificBitPatterns) {
+ float x = as_float(0xc70d39a1);
+ EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance));
+}
+
+// For small values, sin(x) is x.
+TEST(SinfTest, SmallValues) {
+ uint32_t bits = 0x17800000;
+ float x = as_float(bits);
+ float result = __llvm_libc::sinf(x);
+ EXPECT_TRUE(mpfr::equalsSin(x, result, tolerance));
+ EXPECT_EQ(bits, as_uint32_bits(result));
+
+ bits = 0x00400000;
+ x = as_float(bits);
+ result = __llvm_libc::sinf(x);
+ EXPECT_TRUE(mpfr::equalsSin(x, result, tolerance));
+ EXPECT_EQ(bits, as_uint32_bits(result));
+}
+
+// SDCOMP-26094: check sinf in the cases for which the range reducer
+// returns values furthest beyond its nominal upper bound of pi/4.
+TEST(SinfTest, SDCOMP_26094) {
+ for (uint32_t v : sdcomp26094Values) {
+ float x = as_float(v);
+ EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance));
+ }
+}
diff --git a/libc/utils/CMakeLists.txt b/libc/utils/CMakeLists.txt
index e0aa20e19a3e..1e85d05507cc 100644
--- a/libc/utils/CMakeLists.txt
+++ b/libc/utils/CMakeLists.txt
@@ -1,5 +1,6 @@
add_subdirectory(CPP)
add_subdirectory(HdrGen)
+add_subdirectory(MPFRWrapper)
add_subdirectory(testutils)
add_subdirectory(UnitTest)
add_subdirectory(benchmarks)
diff --git a/libc/utils/MPFRWrapper/CMakeLists.txt b/libc/utils/MPFRWrapper/CMakeLists.txt
new file mode 100644
index 000000000000..aa03724aee46
--- /dev/null
+++ b/libc/utils/MPFRWrapper/CMakeLists.txt
@@ -0,0 +1,17 @@
+try_compile(
+ LIBC_TESTS_CAN_USE_MPFR
+ ${CMAKE_CURRENT_BINARY_DIR}
+ SOURCES
+ ${CMAKE_CURRENT_SOURCE_DIR}/check_mpfr.cpp
+ LINK_LIBRARIES
+ -lmpfr -lgmp
+)
+
+if(LIBC_TESTS_CAN_USE_MPFR)
+ add_library(libcMPFRWrapper
+ MPFRUtils.cpp
+ MPFRUtils.h
+ )
+else()
+ message(WARNING "Math tests using MPFR will be skipped.")
+endif()
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
new file mode 100644
index 000000000000..7bd849934fc7
--- /dev/null
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -0,0 +1,97 @@
+//===-- Utils which wrap MPFR ---------------------------------------------===//
+//
+// 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 "MPFRUtils.h"
+
+#include <iostream>
+#include <mpfr.h>
+
+namespace __llvm_libc {
+namespace testing {
+namespace mpfr {
+
+class MPFRNumber {
+ // A precision value which allows sufficiently large additional
+ // precision even compared to double precision floating point values.
+ static constexpr unsigned int mpfrPrecision = 96;
+
+ mpfr_t value;
+
+public:
+ MPFRNumber() { mpfr_init2(value, mpfrPrecision); }
+
+ explicit MPFRNumber(float x) {
+ mpfr_init2(value, mpfrPrecision);
+ mpfr_set_flt(value, x, MPFR_RNDN);
+ }
+
+ MPFRNumber(const MPFRNumber &other) {
+ mpfr_set(value, other.value, MPFR_RNDN);
+ }
+
+ ~MPFRNumber() { mpfr_clear(value); }
+
+ // Returns true if |other| is within the tolerance value |t| of this
+ // number.
+ bool isEqual(const MPFRNumber &other, const Tolerance &t) {
+ MPFRNumber tolerance(0.0);
+ uint32_t bitMask = 1 << (t.width - 1);
+ for (int exponent = -t.basePrecision; bitMask > 0; bitMask >>= 1) {
+ --exponent;
+ if (t.bits & bitMask) {
+ MPFRNumber delta;
+ mpfr_set_ui_2exp(delta.value, 1, exponent, MPFR_RNDN);
+ mpfr_add(tolerance.value, tolerance.value, delta.value, MPFR_RNDN);
+ }
+ }
+
+ MPFRNumber
diff erence;
+ if (mpfr_cmp(value, other.value) >= 0)
+ mpfr_sub(
diff erence.value, value, other.value, MPFR_RNDN);
+ else
+ mpfr_sub(
diff erence.value, other.value, value, MPFR_RNDN);
+
+ return mpfr_lessequal_p(
diff erence.value, tolerance.value);
+ }
+
+ // These functions are useful for debugging.
+ float asFloat() const { return mpfr_get_flt(value, MPFR_RNDN); }
+ double asDouble() const { return mpfr_get_d(value, MPFR_RNDN); }
+ void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
+
+public:
+ static MPFRNumber cos(float x) {
+ MPFRNumber result;
+ MPFRNumber mpfrX(x);
+ mpfr_cos(result.value, mpfrX.value, MPFR_RNDN);
+ return result;
+ }
+
+ static MPFRNumber sin(float x) {
+ MPFRNumber result;
+ MPFRNumber mpfrX(x);
+ mpfr_sin(result.value, mpfrX.value, MPFR_RNDN);
+ return result;
+ }
+};
+
+bool equalsCos(float input, float libcOutput, const Tolerance &t) {
+ MPFRNumber mpfrResult = MPFRNumber::cos(input);
+ MPFRNumber libcResult(libcOutput);
+ return mpfrResult.isEqual(libcResult, t);
+}
+
+bool equalsSin(float input, float libcOutput, const Tolerance &t) {
+ MPFRNumber mpfrResult = MPFRNumber::sin(input);
+ MPFRNumber libcResult(libcOutput);
+ return mpfrResult.isEqual(libcResult, t);
+}
+
+} // namespace mpfr
+} // namespace testing
+} // namespace __llvm_libc
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
new file mode 100644
index 000000000000..9f56ccc61fe6
--- /dev/null
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -0,0 +1,51 @@
+//===-- MPFRUtils.h ---------------------------------------------*- 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_UTILS_TESTUTILS_MPFRUTILS_H
+#define LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+namespace testing {
+namespace mpfr {
+
+struct Tolerance {
+ // Number of bits used to represent the fractional
+ // part of a value of type 'float'.
+ static constexpr unsigned int floatPrecision = 23;
+
+ // Number of bits used to represent the fractional
+ // part of a value of type 'double'.
+ static constexpr unsigned int doublePrecision = 52;
+
+ // The base precision of the number. For example, for values of
+ // type float, the base precision is the value |floatPrecision|.
+ unsigned int basePrecision;
+
+ unsigned int width; // Number of valid LSB bits in |value|.
+
+ // The bits in the tolerance value. The tolerance value will be
+ // sum(bits[width - i] * 2 ^ (- basePrecision - i)) for |i| in
+ // range [1, width].
+ uint32_t bits;
+};
+
+// Return true if |libcOutput| is within the tolerance |t| of the cos(x)
+// value as evaluated by MPFR.
+bool equalsCos(float x, float libcOutput, const Tolerance &t);
+
+// Return true if |libcOutput| is within the tolerance |t| of the sin(x)
+// value as evaluated by MPFR.
+bool equalsSin(float x, float libcOutput, const Tolerance &t);
+
+} // namespace mpfr
+} // namespace testing
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H
diff --git a/libc/utils/MPFRWrapper/check_mpfr.cpp b/libc/utils/MPFRWrapper/check_mpfr.cpp
new file mode 100644
index 000000000000..6e6282457960
--- /dev/null
+++ b/libc/utils/MPFRWrapper/check_mpfr.cpp
@@ -0,0 +1,8 @@
+#include <mpfr.h>
+
+int main() {
+ mpfr_t x;
+ mpfr_init(x);
+ mpfr_clear(x);
+ return 0;
+}
More information about the libc-commits
mailing list