下面的代码提供了sin() 和acos() 的简单实现，它们应该满足您的准确性要求并且您可能想尝试一下。请注意，您平台上的数学库实现很可能针对该平台的特定硬件功能进行了高度调整，并且可能还以汇编形式编码以实现最高效率，因此不太可能提供不满足硬件细节的简单编译 C 代码更高的性能，即使精度要求从全双精度有所放宽。正如 Viktor Latypov 指出的那样，寻找不需要昂贵调用超越数学函数的算法替代方案也可能是值得的。
在下面的代码中，我尝试使用简单、可移植的结构。如果您的编译器支持rint() 函数[由C99 和C++11 指定]，您可能希望使用它而不是my_rint()。在某些平台上，对floor() 的调用可能会很昂贵，因为它需要动态更改机器状态。函数 my_rint()、sin_core()、cos_core() 和 asin_core() 需要内联以获得最佳性能。您的编译器可能会在高优化级别自动执行此操作（例如，使用 -O3 编译时），或者您可以为这些函数添加适当的内联属性，例如inline 或 __inline 取决于您的工具链。
对您的平台一无所知，我选择了简单的多项式近似，使用 Estrin 方案和 Horner 方案对其进行评估。有关这些评估方案的描述，请参见 Wikipedia：
http://en.wikipedia.org/wiki/Estrin%27s_scheme ,
http://en.wikipedia.org/wiki/Horner_scheme
近似值本身属于极大极小值类型，并且是使用 Remez 算法为此答案自定义生成的：
http://en.wikipedia.org/wiki/Minimax_approximation_algorithm ,
http://en.wikipedia.org/wiki/Remez_algorithm
acos() 的参数归约中使用的标识在 cmets 中注明，对于 sin()，我使用了 Cody/Waite 风格的参数归约，如下书中所述：
W。 J. Cody, W. Waite，基本功能软件手册。普伦蒂斯·霍尔，1980 年
cmets 中提到的误差范围是近似值，未经严格测试或证明。
/* not quite rint(), i.e. results not properly rounded to nearest-or-even */
double my_rint (double x)
{
double t = floor (fabs(x) + 0.5);
return (x < 0.0) ? -t : t;
}
/* minimax approximation to cos on [-pi/4, pi/4] with rel. err. ~= 7.5e-13 */
double cos_core (double x)
{
double x8, x4, x2;
x2 = x * x;
x4 = x2 * x2;
x8 = x4 * x4;
/* evaluate polynomial using Estrin's scheme */
return (-2.7236370439787708e-7 * x2 + 2.4799852696610628e-5) * x8 +
(-1.3888885054799695e-3 * x2 + 4.1666666636943683e-2) * x4 +
(-4.9999999999963024e-1 * x2 + 1.0000000000000000e+0);
}
/* minimax approximation to sin on [-pi/4, pi/4] with rel. err. ~= 5.5e-12 */
double sin_core (double x)
{
double x4, x2, t;
x2 = x * x;
x4 = x2 * x2;
/* evaluate polynomial using a mix of Estrin's and Horner's scheme */
return ((2.7181216275479732e-6 * x2 - 1.9839312269456257e-4) * x4 +
(8.3333293048425631e-3 * x2 - 1.6666666640797048e-1)) * x2 * x + x;
}
/* minimax approximation to arcsin on [0, 0.5625] with rel. err. ~= 1.5e-11 */
double asin_core (double x)
{
double x8, x4, x2;
x2 = x * x;
x4 = x2 * x2;
x8 = x4 * x4;
/* evaluate polynomial using a mix of Estrin's and Horner's scheme */
return (((4.5334220547132049e-2 * x2 - 1.1226216762576600e-2) * x4 +
(2.6334281471361822e-2 * x2 + 2.0596336163223834e-2)) * x8 +
(3.0582043602875735e-2 * x2 + 4.4630538556294605e-2) * x4 +
(7.5000364034134126e-2 * x2 + 1.6666666300567365e-1)) * x2 * x + x;
}
/* relative error < 7e-12 on [-50000, 50000] */
double my_sin (double x)
{
double q, t;
int quadrant;
/* Cody-Waite style argument reduction */
q = my_rint (x * 6.3661977236758138e-1);
quadrant = (int)q;
t = x - q * 1.5707963267923333e+00;
t = t - q * 2.5633441515945189e-12;
if (quadrant & 1) {
t = cos_core(t);
} else {
t = sin_core(t);
}
return (quadrant & 2) ? -t : t;
}
/* relative error < 2e-11 on [-1, 1] */
double my_acos (double x)
{
double xa, t;
xa = fabs (x);
/* arcsin(x) = pi/2 - 2 * arcsin (sqrt ((1-x) / 2))
* arccos(x) = pi/2 - arcsin(x)
* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2))
*/
if (xa > 0.5625) {
t = 2.0 * asin_core (sqrt (0.5 * (1.0 - xa)));
} else {
t = 1.5707963267948966 - asin_core (xa);
}
/* arccos (-x) = pi - arccos(x) */
return (x < 0.0) ? (3.1415926535897932 - t) : t;
}