仅使用本机精度运算当然可以计算 [0, π] 上的余弦，任何所需的误差范围 >= 0.5 ulp。然而，目标越接近正确舍入的函数，就越需要更多的前期设计工作和运行时的计算工作。
超越函数的实现通常包括参数减少、核心近似、抵消参数减少的最终修复。在参数减少涉及减法的情况下，需要通过显式或隐式使用更高的精度来避免灾难性抵消。隐式技术可以设计为仅依赖于本机精度计算，例如通过将像 π 这样的常数拆分为未计算的总和，例如1.57079637e+0f - 4.37113883e-8f使用 IEEE-754 时binary32（单精度）。
当硬件提供融合乘加 (FMA) 运算时，通过本机精度计算实现高精度会容易得多。 OP 没有指定他们的目标平台是否提供此操作，因此我将首先展示一种非常简单的方法，仅依靠乘法和加法提供中等精度（最大误差 float映射到 IEEE-754binary32 format.
以下内容基于存档的博客文章 https://archive.ph/VyyYh作者：Colin Wallace，标题为“用切比雪夫多项式将 sin(x) 近似为 5 ULP”。它建议通过使用 sin(x)/(x*(x²-π²)) 的 x² 多项式来近似 [-π, π] 上的正弦，然后将其乘以 x*(x²-π²)。更准确地计算 a2-b2 的标准技巧是将其重写为 (a-b) * (a+b)。将 π 表示为两个浮点数 pi_high 和 pi_low 的未计算和，可以避免减法过程中发生灾难性抵消，从而将计算 x²-π² 变为((x - pi_hi) - pi_lo) * ((x + pi_hi) + pi_lo).
理想情况下，多项式核心近似应使用极小极大近似，其中min形象化了max我的错误。我在这里已经这样做了。可以使用 Maple 或 Mathematics 等各种标准工具来实现此目的，或者基于 Remez 算法创建自己的代码。
对于 [0, PI] 上的余弦计算，我们可以利用 cos (t) = sin (π/2 - t) 这一事实。将 x = (π/2 - t) 代入 x * (x - π/2) * (x + π/2) 得到 (π/2 - t) * (3π/2 - t) * (-π/2 -t）。与以前一样，常量可以分为高部分和低部分（或头部分和尾部分，使用另一种常见的习惯用法）。
/* Approximate cosine on [0, PI] with maximum error of 5.081154 ulp */
float cosine (float x)
{
const float half_pi_hi = 1.57079637e+0f; // 0x1.921fb6p+0
const float half_pi_lo = -4.37113883e-8f; // -0x1.777a5cp-25
const float three_half_pi_hi = 4.71238899e+0f; // 0x1.2d97c8p+2
const float three_half_pi_lo = -1.19248806e-8f; // -0x1.99bc5cp-27
float p, s, hpmx, thpmx, nhpmx;
/* cos(x) = sin (pi/2 - x) = sin (hpmx) */
hpmx = (half_pi_hi - x) + half_pi_lo; // pi/2 - x
thpmx = (three_half_pi_hi - x) + three_half_pi_lo; // 3*pi/2 - x
nhpmx = (-half_pi_hi - x) - half_pi_lo; // -pi/2 - x
/* P(hpmx*hpmx) ~= sin (hpmx) / (hpmx * (hpmx * hpmx - pi * pi)) */
s = hpmx * hpmx;
p = 1.32823530e-10f;// 0x1.241500p-33
p = p * s - 2.33173445e-8f; // -0x1.9096c4p-26
p = p * s + 2.52237896e-6f; // 0x1.528c48p-19
p = p * s - 1.73501656e-4f; // -0x1.6bdbfep-13
p = p * s + 6.62087509e-3f; // 0x1.b1e7dap-8
p = p * s - 1.01321183e-1f; // -0x1.9f02f6p-4
return hpmx * nhpmx * thpmx * p;
}
下面我展示了一种经典方法，它首先将参数减少为 [-π/4, π/4]，同时记录象限。然后，象限告诉我们是否需要在此主要近似区间上计算正弦或余弦的多项式近似，以及是否需要翻转最终结果的符号。此代码假设目标平台支持 IEEE-754 指定的 FMA 操作，并且它是通过标准 C 函数映射的fmaf()对于单精度。
该代码很简单，除了使用舍入模式到最近或偶数的浮点到整数转换来计算象限，该转换是通过“幻数加法”方法执行的，并与 2/ 的乘法相结合π（相当于除以 π/2）。最大误差小于1.5 ulps。
/* compute cosine on [0, PI] with maximum error of 1.429027 ulp */
float my_cosf (float a)
{
const float half_pi_hi = 1.57079637e+0f; // 0x1.921fb6p+0
const float half_pi_lo = -4.37113883e-8f; // -0x1.777a5cp-25
float c, j, r, s, sa, t;
int i;
/* subtract closest multiple of pi/2 giving reduced argument and quadrant */
j = fmaf (a, 6.36619747e-1f, 12582912.f) - 12582912.f; // 2/pi, 1.5 * 2**23
a = fmaf (j, -half_pi_hi, a);
a = fmaf (j, -half_pi_lo, a);
/* phase shift of pi/2 (one quadrant) for cosine */
i = (int)j;
i = i + 1;
sa = a * a;
/* Approximate cosine on [-PI/4,+PI/4] with maximum error of 0.87444 ulp */
c = 2.44677067e-5f; // 0x1.9a8000p-16
c = fmaf (c, sa, -1.38877297e-3f); // -0x1.6c0efap-10
c = fmaf (c, sa, 4.16666567e-2f); // 0x1.555550p-5
c = fmaf (c, sa, -5.00000000e-1f); // -0x1.000000p-1
c = fmaf (c, sa, 1.00000000e+0f); // 1.00000000p+0
/* Approximate sine on [-PI/4,+PI/4] with maximum error of 0.64196 ulp */
s = 2.86567956e-6f; // 0x1.80a000p-19
s = fmaf (s, sa, -1.98559923e-4f); // -0x1.a0690cp-13
s = fmaf (s, sa, 8.33338592e-3f); // 0x1.111182p-7
s = fmaf (s, sa, -1.66666672e-1f); // -0x1.555556p-3
t = a * sa;
s = fmaf (s, t, a);
/* select sine approximation or cosine approximation based on quadrant */
r = (i & 1) ? c : s;
/* adjust sign based on quadrant */
r = (i & 2) ? (0.0f - r) : r;
return r;
}
事实证明，在这种特殊情况下，使用 FMA 在准确性方面只提供了很小的好处。如果我将呼叫替换为fmaf(a,b,c) with ((a)*(b)+(c))，最大误差最小增加至 1.451367 ulps，即保持在 1.5 ulps 以下。
