Add comments, update docs

Author: wldfngrs

libc/docs/headers/math/index.rst

@@ -346,7 +346,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | sqrt      | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.7.10              | F.10.4.10                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| tan       | |check|          | |check|         |                        |                      |                        | 7.12.4.7               | F.10.1.7                   |
+| tan       | |check|          | |check|         |                        | |check|              |                        | 7.12.4.7               | F.10.1.7                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | tanh      | |check|          |                 |                        | |check|              |                        | 7.12.5.6               | F.10.2.6                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

libc/src/math/generic/sincosf16_utils.h

@@ -49,22 +49,20 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
 
 // Recall, range reduction:
 //   k = round(x * 32/pi)
-//   y = x * 32/pi - k
 //
-// The constant 0x1.45f306dc9c883p3 is 32/pi rounded to double-precision.
-// 32/pi is generated by Sollya with the following commands:
-// > display = hexadecimal;
-// > round(32/pi, D, RN);
-//
-// The precision choice of 'double' is to minimize rounding errors
-// in this initial scaling step, preserving enough bits so errors accumulated
-// while computing the subtraction: y = x * 32/pi - round(x * 32/pi)
+// The precision choice of 'double' in the following function is to minimize
+// rounding errors in this initial scaling step,
+// preserving enough bits so errors accumulated while computing the subtraction:
+// y = x * 32/pi - round(x * 32/pi)
 // are beyond the least-significant bit of single-precision used during
 // further intermediate computation.
 LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
-  double prod = x * 0x1.45f306dc9c883p3;
-  double kd = fputil::nearest_integer(prod);
+  // Generated by Sollya with:
+  // > D(32/pi);
+  constexpr double THIRTYTWO_OVER_PI = 0x1.45f306dc9c883p3;
 
+  double prod = x * THIRTYTWO_OVER_PI;
+  double kd = fputil::nearest_integer(prod);
   y = static_cast<float>(prod - kd);
 
   return static_cast<int32_t>(kd);

libc/src/math/generic/tanf16.cpp

@@ -50,9 +50,10 @@ LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
 
   // |x| <= 0x1.d1p-5
   if (LIBC_UNLIKELY(x_abs <= 0x2b44)) {
+    // |x| <= 1.398p-11
     if (LIBC_UNLIKELY(x_abs <= 0x10e6)) {
       // tan(+/-0) = +/-0
-      if (LIBC_UNLIKELY(x_abs == 0U))
+      if (LIBC_UNLIKELY(x_abs == 0))
         return x;
 
       int rounding = fputil::quick_get_round();
@@ -64,20 +65,22 @@ LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
       if ((xbits.is_pos() && rounding == FE_UPWARD) ||
           (xbits.is_neg() && rounding == FE_DOWNWARD))
         return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
-      else
-        return x;
+      return x;
     }
 
     float xsq = xf * xf;
 
+    // Degree-6 minimax odd polynomial of tan(x) generated by Sollya with:
+    // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6|], [|1, SG...|], [0, pi/32]);
     float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f,
                                     0x1.be80f6p-5f);
 
     return fputil::cast<float16>(xf * result);
   }
 
-  // tan(+/-inf)= NaN, and tan(NaN) = NaN
+  // tan(+/-inf) = NaN, and tan(NaN) = NaN
   if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+    // x = +/-inf
     if (x_abs == 0x7c00) {
       fputil::set_errno_if_required(EDOM);
       fputil::raise_except_if_required(FE_INVALID);