[libc] Move implementations of cosf, sinf, sincosf to src/math directory.

[sivachandra]

Summary

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.

[sivachandra]

This patch is not ready for review yet. For one, it does not have any tests. I am only sharing to showcase the move early.

[abrachet]

libc/src/math/CMakeLists.txt

10

Will this properly propagate to targets which depend on math_utils?

libc/src/math/cosf.cpp

20–22

Maybe we could remove the two spaces after the period. This isn't common in the rest of our comments

libc/src/math/math_utils.h

2

Should this file have a comment header?

libc/src/math/sincosf.h

2

sinf -> sincosf

libc/src/math/sincosf_utils.h

34–37

We can put the data from sincosf_data and just make them constexpr perhaps.

[abrachet]

This patch is not ready for review yet. For one, it does not have any tests. I am only sharing to showcase the move early.

It's probably worth discussing how much we expect we will alter these files. Is it fine to just accept that we didn't write them and that they will just be very different from the rest of libc?

[sivachandra]

In D76825#1942971, @abrachet wrote:

It's probably worth discussing how much we expect we will alter these files. Is it fine to just accept that we didn't write them and that they will just be very different from the rest of libc?

While I do intend to make any functional changes in this move patch, we absolutely want the code to follow our implementation standards, just like the rest of the code. So, thanks for already pointing out the items missed. For now, all I have done was to run clang-format. I will fix the rest with a new upload.

[sivachandra]

In D76825#1944581, @sivachandra wrote:

In D76825#1942971, @abrachet wrote:

It's probably worth discussing how much we expect we will alter these files. Is it fine to just accept that we didn't write them and that they will just be very different from the rest of libc?

While I do intend to make any functional changes in this move patch, we absolutely want the code to follow our implementation standards, just like the rest of the code. So, thanks for already pointing out the items missed. For now, all I have done was to run clang-format. I will fix the rest with a new upload.

I mean't to say, "do NOT intend".

[PaulkaToast]

Won't deleting the files in the AOR directory cause the make check on the buildbots to fail? It might make sense to leave the files and delete them once all functions have been migrated?

[sivachandra]

In D76825#1945462, @PaulkaToast wrote:

Won't deleting the files in the AOR directory cause the make check on the buildbots to fail? It might make sense to leave the files and delete them once all functions have been migrated?

Yes, it will. I plan to move the tests also in a way that make check does not break.

[sivachandra]

Not yet ready for review. Changes in this update:

• Add test infrastructure for math functions which uses MPFR.
• Add tests for cosf and sinf

[sivachandra]

This is not yet ready for review. I have added bunch of tests but more tests need to be ported from the AOR directory. Also, the test infrastructure needs a lot of comments which I will add in my next update.

[sivachandra]

• All tests I planned to add in this round are now added.
• ULP related tests have not been ported yet, but that can be done in a different round instead of cramming up this patch with everything.
• Added comments which I planned to add to explain the various pieces of the test infrastructure.

[sivachandra]

Rebase.

[abrachet]

libc/src/__support/common.h.def

16

I think it might be more common in llvm to cast to void rather than use a macro like this. But I don't have a strong preference.

libc/src/math/cosf.cpp

60–62

Remove this else, just return invalidf(y)

libc/src/math/sincosf.cpp

14

Should we use " instead of <

libc/src/math/sincosf.h

2

Only 2 -

10

SINCOSF

libc/src/math/sincosf_utils.h

10

SINCOSF_UTILS

75

I think we can declare all of these where they are initialized

81

Not sure why they choose the name x7, isn't x5 more appropriate.

85

No need for this else

107

Declare r down here.

libc/src/math/sinf.cpp

34–36

I don't want to comment on the substance of the code too much, but I'll risk making a fool out of myself for this one. These lines logically don't make sense if we are in the first if then we will return y and s has no impact on that. So I wonder what purpose these lines serve.

64

Remove this else

libc/test/src/math/cosf_test.cpp

31

; -> .

38–40

These could become EXPECT_THAT(Succeeds(FloatBits::isQNan(as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::QNan))))), true); and no longer need the EXPECT_EQ(llvmlibc_errno, 0);

72–83

Is this meant to be here?

libc/test/src/math/sincosf_test.cpp

31

. not ;

libc/test/src/math/sinf_test.cpp

31

; -> .

libc/utils/MPFRWrapper/MPFRUtils.cpp

45

for (int exponent = -t.basePrecision; bitMask ... maybe

libc/utils/MPFRWrapper/MPFRUtils.h

39–40

This could be a namespace or at least a struct

[sivachandra]

• Address comments.
• Add tests for small input values to sinf/cosf/sincosf functions.

[sivachandra]

Thanks a lot for the detailed review.

I did not address some of the comments. The plan is to quickly absorb the AOR code into the normal LLVM-libc way and try to bring AOR developers to work with LLVM-libc. To that extant, I want to keep their structuring as intact as possible so that they do not have to go through a learning phase to understand their own code. We do want the higher-level structuring to adhere to LLVM libc's style so I tried to get that in shape with this change.

libc/src/math/CMakeLists.txt

10

That's a good point. It does not and we have to add the entrypoint deps explicitly. Not ideal but I will try to work something out with my change to propagate entrypoint deps.

libc/src/math/sinf.cpp

34–36

Neither do I know. Like I have explained, I want to keep it as is so I that we don't increase the unfamiliarity for Arm developers.

libc/test/src/math/cosf_test.cpp

38–40

I started with that but I found that so much more unreadable in this case. That said, I want to do better here. For example, if something fails, we aren't printing information about what two values/bit patterns were compared. I have some ideas which I want to start incorporating with more functions that I convert.

72–83

Yup, I was playing with tolerance and forgot to un-comment aferwards.

[sivachandra]

• Address one missed comment.

[sivachandra]

There are few hexadecimal floating literals in this patch which moved along with the AOR implementations. They trigger warnings of this form:

warning: hexadecimal floating literals are a C++17 feature [-Wc++17-extensions]

I think I can fix them by using their bit patterns as is already done in other places. Will do so in the next pass.

[abrachet]

In D76825#1979742, @sivachandra wrote:

Thanks a lot for the detailed review.

I did not address some of the comments. The plan is to quickly absorb the AOR code into the normal LLVM-libc way and try to bring AOR developers to work with LLVM-libc. To that extant, I want to keep their structuring as intact as possible so that they do not have to go through a learning phase to understand their own code. We do want the higher-level structuring to adhere to LLVM libc's style so I tried to get that in shape with this change.

Sounds good. LGTM

[phosek]

LGTM, I had the same concern as @abrachet regarding some of the readability bits (curly braces, minimizing variable scopes) but keeping the code as close to original SGTM as well. If we decide to make those changes, we could probably automate them later using a tool like clang-tidy.

libc/AOR_v20.02/math/cosf.c

-/*
- * 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);
-}

libc/AOR_v20.02/math/sincosf.h

-/*
- * 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;
-}

libc/AOR_v20.02/math/sincosf.c

-/*
- * 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
-    }
-}

libc/AOR_v20.02/math/sincosf_data.c

-/*
- * 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
-};

libc/AOR_v20.02/math/sinf.c

-/*
- * 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);
-}

libc/AOR_v20.02/math/test/testcases/directed/cosf.tst

-; 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

libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst

-; 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

libc/AOR_v20.02/math/test/testcases/directed/sinf.tst

-; 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

libc/AOR_v20.02/math/test/testcases/random/float.tst

 !! single.tst - Random test case specification for SP functions
 !!
 !! 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
 
-test sinf 10000
-test cosf 10000
-test sincosf_sinf 5000
-test sincosf_cosf 5000
 test tanf 10000
 test expf 10000
 test exp2f 10000
 test logf 10000
 test log2f 10000
 test powf 10000

libc/config/linux/api.td

 include "config/public_api.td"
 
+include "spec/gnu_ext.td"
 include "spec/linux.td"
 include "spec/posix.td"
 include "spec/stdc.td"
 
 def SizeT : TypeDecl<"size_t"> {
   let Decl = [{
     #define __need_size_t
     #include <stddef.h>
let Macros = [
     SimpleMacroDef<"INFINITY", "__builtin_inff()">,
     SimpleMacroDef<"NAN", "__builtin_nanf(\"\")">,
 
     IsFiniteMacro,
     IsInfMacro,
     IsNanMacro,
   ];
   let Functions = [
+   "cosf",
    "round",
+   "sincosf",
+   "sinf",
   ];
 }
 
 def StringAPI : PublicAPI<"string.h"> {
   let Functions = [
     "memcpy",
     "memmove",
     "memcmp",

libc/lib/CMakeLists.txt

DEPENDS
     libc.src.threads.thrd_create
     libc.src.threads.thrd_join
 )
 
 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(
   llvmlibc_redirectors
   DEPENDS
     round_redirector
 )

libc/src/__support/common.h.def

 //===-- Common internal contructs -----------------------------------------===//
 //
 // 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_SUPPORT_COMMON_H
 #define LLVM_LIBC_SUPPORT_COMMON_H
 
 #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})
 
 #endif // LLVM_LIBC_SUPPORT_COMMON_H

libc/src/math/CMakeLists.txt

+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
   SRCS
     round.cpp
   HDRS
     round.h
 )
 
 add_redirector_object(
   round_redirector
   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
+)

libc/src/math/cosf.h

+//===-- 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

libc/src/math/cosf.cpp

+//===-- 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

libc/src/math/math_utils.h

+//===-- 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

libc/src/math/sincosf.h

+//===-- 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

libc/src/math/sincosf.cpp

+//===-- 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

libc/src/math/sincosf_data.cpp

+//===-- 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

libc/src/math/sincosf_utils.h

+//===-- 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

libc/src/math/sinf.h

+//===-- 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

libc/src/math/sinf.cpp

+//===-- 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

libc/test/src/CMakeLists.txt

 add_subdirectory(assert)
 add_subdirectory(errno)
+add_subdirectory(math)
 add_subdirectory(signal)
 add_subdirectory(stdio)
 add_subdirectory(stdlib)
 add_subdirectory(string)
 add_subdirectory(sys)
 add_subdirectory(threads)

libc/test/src/math/CMakeLists.txt

+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
+)

libc/test/src/math/cosf_test.cpp

+//===-- 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));
+  }
+}

libc/test/src/math/float.h

+//===-- 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

libc/test/src/math/sdcomp26094.h

+//===-- 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

libc/test/src/math/sincosf_test.cpp

+//===-- 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));
+  }
+}

libc/test/src/math/sinf_test.cpp

+//===-- 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));
+  }
+}

libc/utils/CMakeLists.txt

 add_subdirectory(CPP)
 add_subdirectory(HdrGen)
+add_subdirectory(MPFRWrapper)
 add_subdirectory(testutils)
 add_subdirectory(UnitTest)
 add_subdirectory(benchmarks)

libc/utils/MPFRWrapper/CMakeLists.txt

+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()

libc/utils/MPFRWrapper/MPFRUtils.h

+//===-- 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

libc/utils/MPFRWrapper/MPFRUtils.cpp

+//===-- 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 difference;
+    if (mpfr_cmp(value, other.value) >= 0)
+      mpfr_sub(difference.value, value, other.value, MPFR_RNDN);
+    else
+      mpfr_sub(difference.value, other.value, value, MPFR_RNDN);
+
+    return mpfr_lessequal_p(difference.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

libc/utils/MPFRWrapper/check_mpfr.cpp

+#include <mpfr.h>
+
+int main() {
+  mpfr_t x;
+  mpfr_init(x);
+  mpfr_clear(x);
+  return 0;
+}