[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

Author: Siva Chandra Reddy

libc/AOR_v20.02/math/cosf.c

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

libc/AOR_v20.02/math/sincosf.c

@@ -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",
   ];
 }
 

libc/AOR_v20.02/math/sincosf.h

@@ -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(

libc/AOR_v20.02/math/sincosf_data.c

@@ -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})

libc/AOR_v20.02/math/sinf.c

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

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

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

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

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

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

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

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

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

libc/config/linux/api.td

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

libc/lib/CMakeLists.txt

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

libc/src/__support/common.h.def

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

libc/src/math/CMakeLists.txt

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

libc/src/math/cosf.cpp

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

libc/src/math/cosf.h

@@ -1,5 +1,6 @@
 add_subdirectory(assert)
 add_subdirectory(errno)
+add_subdirectory(math)
 add_subdirectory(signal)
 add_subdirectory(stdio)
 add_subdirectory(stdlib)

libc/src/math/math_utils.h

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

libc/src/math/sincosf.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));
+  }
+}

libc/src/math/sincosf.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

libc/src/math/sincosf_data.cpp

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

libc/src/math/sincosf_utils.h

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

libc/src/math/sinf.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));
+  }
+}

libc/src/math/sinf.h

@@ -1,5 +1,6 @@
 add_subdirectory(CPP)
 add_subdirectory(HdrGen)
+add_subdirectory(MPFRWrapper)
 add_subdirectory(testutils)
 add_subdirectory(UnitTest)
 add_subdirectory(benchmarks)

libc/test/src/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()

libc/test/src/math/CMakeLists.txt

@@ -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 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/test/src/math/cosf_test.cpp

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

libc/test/src/math/float.h

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