[libc] Add x86_64 implementations of double precision cos, sin and tan.

[sivachandra]

Summary

The implementations use the x86_64 FPU instructions. These instructions
are extremely slow compared to a polynomial based software
implementation. Also, their accuracy falls drastically once the input
goes beyond 2PI. To improve both the speed and accuracy, we will be
taking the following approach going forward:

1. As a follow up to this CL, we will implement a range reduction algorithm

which will expand the accuracy to the entire double precision range.

2. After that, we will replace the HW instructions with a polynomial

implementation to improve the run time.

After step 2, the implementations will be accurate, performant and target
architecture independent.

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
 
     # math.h entrypoints
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
+    libc.src.math.cos
     libc.src.math.cosf
     libc.src.math.expf
     libc.src.math.exp2f
     libc.src.math.fabs
     libc.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.remquo
     libc.src.math.remquol
     libc.src.math.rint
     libc.src.math.rintf
     libc.src.math.rintl
     libc.src.math.round
     libc.src.math.roundf
     libc.src.math.roundl
+    libc.src.math.sin
     libc.src.math.sincosf
     libc.src.math.sinf
     libc.src.math.sqrt
     libc.src.math.sqrtf
     libc.src.math.sqrtl
+    libc.src.math.tan
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 )
 
 if(LLVM_LIBC_FULL_BUILD)
   list(APPEND TARGET_LIBC_ENTRYPOINTS
     # assert.h entrypoints

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"logb", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logbf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"logbl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"modf", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoublePtr>]>,
           FunctionSpec<"modff", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatPtr>]>,
           FunctionSpec<"modfl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoublePtr>]>,
 
+          FunctionSpec<"cos", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"cosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          FunctionSpec<"sin", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"sinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          FunctionSpec<"tan", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
 
           FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
           FunctionSpec<"remainderf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"remainder", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"remainderl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 

libc/src/math/CMakeLists.txt

 add_math_entrypoint_object(ceil)
 add_math_entrypoint_object(ceilf)
 add_math_entrypoint_object(ceill)
 
 add_math_entrypoint_object(copysign)
 add_math_entrypoint_object(copysignf)
 add_math_entrypoint_object(copysignl)
 
+add_math_entrypoint_object(cos)
 add_math_entrypoint_object(cosf)
 
 add_math_entrypoint_object(expf)
 
 add_math_entrypoint_object(exp2f)
 
 add_math_entrypoint_object(fabs)
 add_math_entrypoint_object(fabsf)

 add_math_entrypoint_object(rintl)
 
 add_math_entrypoint_object(round)
 add_math_entrypoint_object(roundf)
 add_math_entrypoint_object(roundl)
 
 add_math_entrypoint_object(sincosf)
 
+add_math_entrypoint_object(sin)
 add_math_entrypoint_object(sinf)
 
 add_math_entrypoint_object(sqrt)
 add_math_entrypoint_object(sqrtf)
 add_math_entrypoint_object(sqrtl)
 
+add_math_entrypoint_object(tan)
+
 add_math_entrypoint_object(trunc)
 add_math_entrypoint_object(truncf)
 add_math_entrypoint_object(truncl)

libc/src/math/cos.h

+//===-- Implementation header for cos ---------------------------*- 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_COS_H
+#define LLVM_LIBC_SRC_MATH_COS_H
+
+namespace __llvm_libc {
+
+double cos(double x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_COS_H

libc/src/math/sin.h

+//===-- Implementation header for sin ---------------------------*- 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_SIN_H
+#define LLVM_LIBC_SRC_MATH_SIN_H
+
+namespace __llvm_libc {
+
+double sin(double x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_SIN_H

libc/src/math/tan.h

+//===-- Implementation header for tan ---------------------------*- 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_TAN_H
+#define LLVM_LIBC_SRC_MATH_TAN_H
+
+namespace __llvm_libc {
+
+double tan(double x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_TAN_H

libc/src/math/x86_64/CMakeLists.txt

+add_entrypoint_object(
+  cos
+  SRCS
+    cos.cpp
+  HDRS
+    ../cos.h
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  sin
+  SRCS
+    sin.cpp
+  HDRS
+    ../sin.h
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  tan
+  SRCS
+    tan.cpp
+  HDRS
+    ../tan.h
+  COMPILE_OPTIONS
+    -O2
+)

libc/src/math/x86_64/cos.cpp

+//===-- Implementation of the cos function for x86_64 ---------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/cos.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(double, cos, (double x)) {
+  double result;
+  __asm__ __volatile__("fcos" : "=t"(result) : "f"(x));
+  return result;
+}
+
+} // namespace __llvm_libc

libc/src/math/x86_64/sin.cpp

+//===-- Implementation of the sin function for x86_64 ---------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/sin.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(double, sin, (double x)) {
+  double result;
+  __asm__ __volatile__("fsin" : "=t"(result) : "f"(x));
+  return result;
+}
+
+} // namespace __llvm_libc

libc/src/math/x86_64/tan.cpp

+//===-- Implementation of the tan function for x86_64 ---------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/tan.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(double, tan, (double x)) {
+  double result;
+  double one;
+  // The fptan instruction pushes the number 1 on to the FP stack after
+  // computing tan. So, we read out the one before popping the actual result.
+  __asm__ __volatile__("fptan" : "=t"(one) : "f"(x));
+  __asm__ __volatile__("fstpl %0" : "=m"(result));
+  return result;
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   DEPENDS
     libc.include.errno
     libc.src.math.cosf
     libc.utils.CPP.standalone_cpp
     libc.utils.FPUtil.fputil
 )
 
 add_fp_unittest(
+  cos_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    cos_test.cpp
+  DEPENDS
+    libc.src.math.cos
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
   sinf_test
   NEED_MPFR
   SUITE
     libc_math_unittests
   SRCS
     sinf_test.cpp
   HDRS
     sdcomp26094.h
   DEPENDS
     libc.include.errno
     libc.src.math.sinf
     libc.utils.CPP.standalone_cpp
     libc.utils.FPUtil.fputil
 )
 
 add_fp_unittest(
+  sin_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    sin_test.cpp
+  DEPENDS
+    libc.src.math.sin
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
   sincosf_test
   NEED_MPFR
   SUITE
     libc_math_unittests
   SRCS
     sincosf_test.cpp
   HDRS
     sdcomp26094.h
add_fp_unittest(
   SRCS
     fma_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.fma
     libc.utils.FPUtil.fputil
 )
 
+add_fp_unittest(
+  tan_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    tan_test.cpp
+  DEPENDS
+    libc.src.math.tan
+    libc.utils.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/cos_test.cpp

+//===-- Unittests for cos -------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/cos.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <math.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(double)
+
+TEST(LlvmLibccosTest, Range) {
+  static constexpr double _2pi = 6.283185307179586;
+  constexpr UIntType count = 10000000;
+  constexpr UIntType step = UIntType(-1) / count;
+  for (UIntType i = 0, v = 0; i <= count; ++i, v += step) {
+    double x = double(FPBits(v));
+    // TODO: Expand the range of testing after range reduction is implemented.
+    if (isnan(x) || isinf(x) || x > _2pi || x < -_2pi)
+      continue;
+
+    ASSERT_MPFR_MATCH(mpfr::Operation::Cos, x, __llvm_libc::cos(x), 1.0);
+  }
+}

libc/test/src/math/sin_test.cpp

+//===-- Unittests for sin -------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/sin.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <math.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(double)
+
+TEST(LlvmLibcSinTest, Range) {
+  static constexpr double _2pi = 6.283185307179586;
+  constexpr UIntType count = 10000000;
+  constexpr UIntType step = UIntType(-1) / count;
+  for (UIntType i = 0, v = 0; i <= count; ++i, v += step) {
+    double x = double(FPBits(v));
+    // TODO: Expand the range of testing after range reduction is implemented.
+    if (isnan(x) || isinf(x) || x > _2pi || x < -_2pi)
+      continue;
+
+    ASSERT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sin(x), 1.0);
+  }
+}

libc/test/src/math/tan_test.cpp

+//===-- Unittests for tan -------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/tan.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <math.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(double)
+
+TEST(LlvmLibctanTest, Range) {
+  static constexpr double _2pi = 6.283185307179586;
+  constexpr UIntType count = 10000000;
+  constexpr UIntType step = UIntType(-1) / count;
+  for (UIntType i = 0, v = 0; i <= count; ++i, v += step) {
+    double x = double(FPBits(v));
+    // TODO: Expand the range of testing after range reduction is implemented.
+    if (isnan(x) || isinf(x) || x > _2pi || x < -_2pi)
+      continue;
+
+    ASSERT_MPFR_MATCH(mpfr::Operation::Tan, x, __llvm_libc::tan(x), 1.0);
+  }
+}

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   Ceil,
   Cos,
   Exp,
   Exp2,
   Floor,
   Round,
   Sin,
   Sqrt,
+  Tan,
   Trunc,
   EndUnaryOperationsSingleOutput,
 
   // Operations which take a single floating point nubmer as input
   // but produce two outputs. The first ouput is a floating point
   // number of the same type as the input. The second output is of type
   // 'int'.
   BeginUnaryOperationsTwoOutputs,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber sqrt() const {
     MPFRNumber result;
     mpfr_sqrt(result.value, value, MPFR_RNDN);
     return result;
   }
 
+  MPFRNumber tan() const {
+    MPFRNumber result;
+    mpfr_tan(result.value, value, MPFR_RNDN);
+    return result;
+  }
+
   MPFRNumber trunc() const {
     MPFRNumber result;
     mpfr_trunc(result.value, value);
     return result;
   }
 
   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
     MPFRNumber result(*this);
unaryOperation(Operation op, InputType input) {
   case Operation::Floor:
     return mpfrInput.floor();
   case Operation::Round:
     return mpfrInput.round();
   case Operation::Sin:
     return mpfrInput.sin();
   case Operation::Sqrt:
     return mpfrInput.sqrt();
+  case Operation::Tan:
+    return mpfrInput.tan();
   case Operation::Trunc:
     return mpfrInput.trunc();
   default:
     __builtin_unreachable();
   }
 }
 
 template <typename InputType>