// Computes product from (m+1) to (m+a)
unsigned long long compute_with_for_loop(unsigned long long a, unsigned long long m) {
if (m < 1) {
return 1;
}
if (m < 4) {
unsigned long long result = 1;
for (unsigned long long i = 1; i <= m; i++) {
result *= (a+i);
}
return result;
}
unsigned long long result = 1;
unsigned long long i = 1;
for (; i < m; i+=4) {
result *= (a+i) * (a+i+1) * (a+i+2) * (a+i+3);
}
if (i <= m) {
for (; i <= m; i++) {
result *= (a+i);
}
}
return result;
}其动机是计算(a+m)!/m! = (a+m)(a+m-1)...(m+1)
我知道tgamma()，但我想测试一下它是否更快。
如何强制编译器使用矢量化？我尝试添加-march=native -ftree-vectorize -O3，但.lst文件不包含任何矢量化。
1条答案
按热度按时间cgfeq70w1#
动机是计算（a+m）！/m！=（a+m）（a+m-1）...（m+1）我知道tgamma（），但我想测试一下它是否更快。
这取决于tgamma（）是如何计算的：）
很可能，对于大量的数字和大量的因子，你可能会使用斯特林公式
n! ~= sqrt(2*M_PI*n)(n/M_E)^n给出的阶乘近似（我认为有一个更好的近似），用它来计算n!/m!：您可以使用快速求幂算法使求幂比迭代乘法更快（您不能像在快速求幂算法中那样通过平方来压缩它们）
另一方面，伽马函数是有问题的函数，因为它在所有负整数中具有极点，因此，在不引起缓慢收敛等的情况下，不能容易地开发为大圆上的级数。恕我直言，伽马函数是在点1/2周围计算的，并且具有1/2的半径以完全精度覆盖区间[0，1]而不接近-1处的极点，使用递归定义Gamma（p + 1）= p * Gamma（p）直到p在区间[0，1]中。
但我确信有比我上面解释的更好的算法来计算它，谷歌搜索一下可能会给予你一个更好的解决方案。