我知道这是一个常见的问题，但我还没有找到有用的答案。我基本上是在寻找C++中函数acos的快速近似方法，我想知道是否可以显著地超越标准方法。
但你们中的一些人可能对我的具体问题有深入的了解：我正在编写一个科学程序，需要非常快速。主算法的复杂性归结为计算以下表达式（使用不同的参数多次计算）：
sin( acos(t_1) + acos(t_2) + ... + acos(t_n) )
我有一组已知实数（双精度），表示为t_i，并且n非常小（比6还要小）。 我需要至少达到1e-10的精度。 目前我正在使用标准的C++函数sin和acos。
你认为我能通过某种方式显著提高速度吗？ 对于那些了解一些数学的人，您是否认为将正弦展开以获得一个关于t_i的代数表达式（只涉及平方根）是明智的选择？
谢谢你的回答。
sin()和acos()的简单实现，可以满足您的精度要求并且您可能想要尝试。请注意，您平台上的数学库实现很可能高度针对该平台的特定硬件功能进行了优化，并且可能也是以汇编代码编写的，以达到最大效率，因此即使精度要求从完全双精度放宽，简单的编译C代码不考虑硬件的具体情况也不太可能提供更高的性能。正如Viktor Latypov指出的那样，还值得寻找不需要调用超越函数的算法替代方案。rint()函数，你可能想使用它来代替my_rint()。在某些平台上，调用floor()函数可能会很昂贵，因为它需要动态改变机器状态。函数my_rint()、sin_core()、cos_core()和asin_core()最好是内联以获得最佳性能。你的编译器可能会在高优化级别下（例如使用-O3编译时）自动执行这些操作，或者你可以为这些函数添加适当的内联属性，例如inline或__inline，具体取决于你的工具链。http://en.wikipedia.org/wiki/Minimax_approximation_algorithm ，http://en.wikipedia.org/wiki/Remez_algorithm
对于acos()的参数化简中所用的恒等式已在注释中标注；对于sin()，我使用了Cody/Waite风格的参数化简方法，详见以下书籍：W. J. Cody, W. Waite, Software Manual for the Elementary Functions. Prentice-Hall, 1980/* 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;
}
sin( acos(t1) + acos(t2) + ... + acos(tn) )
sin( acos(x) ) and cos(acos(x))=x
因为
sin(a+b) = cos(a)sin(b)+sin(a)cos(b).
首先
sin( acos(x) ) = sqrt(1-x*x)
对 sqrt 进行 Taylor 级数展开，可以将问题简化为多项式计算。
为了澄清，以下是展开到 n=2、n=3 的级数：
sin( acos(t1) + acos(t2) ) = sin(acos(t1))cos(acos(t2)) + sin(acos(t2))cos(acos(t1) = sqrt(1-t1*t1) * t2 + sqrt(1-t2*t2) * t1
cos( acos(t2) + acos(t3) ) = cos(acos(t2)) cos(acos(t3)) - sin(acos(t2))sin(acos(t3)) = t2*t3 - sqrt(1-t2*t2)*sqrt(1-t3*t3)
sin( acos(t1) + acos(t2) + acos(t3)) =
sin(acos(t1))cos(acos(t2) + acos(t3)) + sin(acos(t2)+acos(t3) )cos(acos(t1)=
sqrt(1-t1*t1) * (t2*t3 - sqrt(1-t2*t2)*sqrt(1-t3*t3)) + (sqrt(1-t2*t2) * t3 + sqrt(1-t3*t3) * t2 ) * t1
等等。
在(-1,1)范围内计算x的平方根可以使用以下方法：
x_0 is some approximation, say, zero
x_(n+1) = 0.5 * (x_n + S/x_n) where S is the argument.
编辑：我是指“巴比伦方法”，详见维基百科的文章：计算平方根的方法。在x（0，1）范围内，进行5-6次迭代即可达到1e-10的精度。
您可以尝试创建查找表，并使用它们代替标准的c++函数，看看是否会有性能提升。
通过将内存和数据流与您的内核对齐，可以实现显著的收益。通常，这比重新创建数学函数所能带来的收益要大得多。考虑如何改进内核运算符与内存之间的访问。
使用缓冲技术可以提高内存访问速度。这取决于您的硬件平台。如果您在DSP上运行此操作，则可以将数据DMA到L2缓存中，并安排指令，以使乘法器单元完全占用。
如果您在通用CPU上，则最好使用对齐的数据，通过预取来填充缓存行。如果您有嵌套循环，则最内层循环应该前后迭代（即向前迭代，然后向后迭代），以利用缓存行等。
您还可以考虑使用多个核心并行计算的方法。如果您可以使用GPU，则可以显着提高性能（尽管精度较低）。
找出代码中耗时最长的部分。只有优化这个区域才能获得最大的收益。
处理器不喜欢分支、跳转或执行路径的变化。通常情况下，处理器必须重新加载指令流水线，这会浪费用于计算的时间。这包括函数调用。
技巧是在循环中放置更多的操作集，并减少迭代次数。
经常使用的变量应该声明为register。虽然许多SO成员已经表示编译器忽略了这个建议，但我发现事实并非如此。最坏的情况是你浪费了一些时间打字。
许多处理器在其指令流水线中有足够的空间来容纳小的for循环。这减少了重新加载指令流水线的时间。
将大型计算循环分成许多小循环。
许多处理器都有数据缓存，这是非常靠近处理器的超快速内存。处理器喜欢从处理器外部的内存中加载数据缓存一次。更多的加载需要时间，这些时间可以用来进行计算。在网上搜索“面向数据设计缓存”。
改变计算的设计，使其能够轻松适应多个处理器。许多CPU具有多个核心，可以并行执行指令。一些处理器具有足够的智能，可以自动将指令委派给它们的多个核心。
一些编译器可以优化代码以进行并行处理（查找您的编译器的编译器选项）。为并行处理设计代码将使编译器更容易进行优化。
打印出函数的汇编语言清单。更改函数的设计以匹配汇编语言或帮助编译器生成更优化的汇编语言。
如果您真的需要更高的效率，请优化汇编语言并将其作为内联汇编代码或单独模块插入。我通常更喜欢后者。
在您的情况下，取泰勒展开的前10项，分别计算并放入单独的变量中：
double term1, term2, term3, term4;
double n, n1, n2, n3, n4;
n = 1.0;
for (i = 0; i < 100; ++i)
{
n1 = n + 2;
n2 = n + 4;
n3 = n + 6;
n4 = n + 8;
term1 = 4.0/n;
term2 = 4.0/n1;
term3 = 4.0/n2;
term4 = 4.0/n3;
然后总结你所有的术语：
result = term1 - term2 + term3 - term4;
// Or try sorting by operation, if possible:
// result = term1 + term3;
// result -= term2 + term4;
n = n4 + 2;
}
让我们先考虑两个术语：
cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
或者 cos(a+b) = cos(a)*cos(b) - sqrt(1-cos(a)*cos(a))*sqrt(1-cos(b)*cos(b))
将cos移到右边
a+b = acos( cos(a)*cos(b) - sqrt(1-cos(a)*cos(a))*sqrt(1-cos(b)*cos(b)) ) ... 1
这里cos(a) = t_1，cos(b) = t_2 a = acos(t_1)，b = acos(t_2)
通过在方程（1）中代入，我们得到
acos(t_1) + acos(t_2) = acos(t_1*t_2 - sqrt(1 - t_1*t_1) * sqrt(1 - t_2*t_2))
在这里，您可以看到您已经将两个acos合并为一个。因此，您可以递归地配对所有acos并形成二叉树。最后，您将得到一个形式为sin(acos(x))的表达式，它等于sqrt(1 - x*x)。
这将改善时间复杂度。
但是，我不确定计算sqrt()的复杂度。
1e-10，你只有四个数量级的差距与double精度。我认为你可能无法在性能和准确度上得到足够的近似，但我可能是错误的。 - Jonas Schäfer