我知道这是一个反复出现的问题，但我还没有真正找到有用的答案。我基本上是在寻找acosC++ 中函数的快速近似值，我想知道我是否可以显着击败标准函数。
但是你们中的一些人可能对我的具体问题有见解：我正在编写一个我需要非常快的科学程序。主要算法的复杂性归结为计算以下表达式（多次使用不同的参数）：
sin( acos(t_1) + acos(t_2) + ... + acos(t_n) )
其中t_i是已知的实数（双）数，并且n非常小（例如小于 6）。我需要至少 1e-10 的精度。我目前正在使用标准sin和acosC++ 函数。
你认为我能以某种方式显着提高速度吗？对于那些了解一些数学的人，您认为扩展该正弦以获得代数表达式t_i（仅涉及平方根）是否明智？
谢谢你的回答。
下面的代码提供了简单的实现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
近似值本身属于 minimax 类型，并且是使用 Remez 算法为该答案自定义生成的：
http://en.wikipedia.org/wiki/Minimax_approximation_algorithm，http://en.wikipedia.org/wiki/Remez_algorithm _ _
参数缩减中使用的恒等式在acos()注释中注明，因为sin()我使用了 Cody/Waite 风格的参数缩减，如下书中所述：
WJ Cody, W. Waite，基本功能软件手册。普伦蒂斯大厅，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) )
boils down to the calculation of
sin( acos(x) ) and cos(acos(x))=x
because
sin(a+b) = cos(a)sin(b)+sin(a)cos(b).
The first thing is
sin( acos(x) ) = sqrt(1-x*x)
Taylor series expansion for the sqrt reduces the problem to polynomial calculations.
To clarify, here's the expansion to 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
and so on.
The sqrt() for x in (-1,1) can be computed using
x_0 is some approximation, say, zero
x_(n+1) = 0.5 * (x_n + S/x_n) where S is the argument.
EDIT: I mean the "Babylonian method", see Wikipedia's article for details. You will need not more than 5-6 iterations to achieve 1e-10 with x in (0,1).
正如 Jonas Wielicki 在评论中提到的那样，您无法做出太多精确的权衡。
您最好的选择是尝试使用函数的处理器内在函数（如果您的编译器还没有这样做）并使用一些数学来减少必要的计算量。
同样非常重要的是将所有内容保持在 CPU 友好的格式中，确保很少有缓存未命中等。
如果您正在计算大量函数，例如acos迁移到 GPU 是您的选择吗？
您可以尝试创建查找表，并使用它们而不是标准的 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 带到 RHS
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() 的复杂性。