c++ - 加速我的球体体积代码(嵌套 while 循环)

Question

我正在尝试加快计算球体体积的代码(参见下面的代码)。这个球体的体积是通过计算小体积段 dv 并将它们相加成一个体积 vol 产生的。

实际上，这段代码只是在我将计算应用于具有对称属性的其他球体(如对象)之前进行完整性检查，因此我应该能够通过在小体积上计算并乘以结束来提高代码速度结果。

将 while (phid<=(360.0/adstep)) 和 while (thetad<=(180.0/adstep)) 中的 360 和 180 分别替换为 180 和 90，然后将所需的计算四分之一，这意味着您可以简单地将最终结果相乘体积增加 4.0。

如果我将 phi 设置为 并将 theta 保留为 180，从而将计算减半，则此方法有效。但是当我将 theta 设置为 90 时它不喜欢它。

输出:

Phi 360, Theta 180
Actual Volume       Calculated Volume   % Difference
4.18879020478639053 4.18878971565348923 0.00001167718922403

Phi 180, Theta 180
4.18879020478639053 4.18878971565618219 0.00001167712493440

Phi 180, Theta 90
4.18879020478639053 4.18586538829648180 0.06987363946500515

您可以在上面看到前两个计算几乎相同(我假设差异是由于精度误差)，而最后一个给出明显不同的结果。嵌套循环会导致问题吗？

任何帮助将不胜感激，因为我在我的研究(谷歌和堆栈溢出)中没有发现任何东西来描述我遇到的问题。

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

int main()
{
double thetar, phir, maxr, vol, dv, vol2, arstep, adstep, rad, rad3, thetad, phid, ntheta, nphi;
cout << fixed << setprecision(17); // Set output precision to defined number of decimal places. Note Double has up to 15 decimal place accuracy

vol=0.0;            // Initialise volume and set at zero
adstep=0.1;     // Steps to rotate angles in degrees
arstep=(adstep/180.0)*M_PI;     // Angle steps in radians
phid=1.0;           // Phi in degrees starting at adstep
maxr = 1.0;         // Radius of the sphere

// Loop to calculate volume

while (phid<=(360.0/adstep))            // Loop over Phi divided by adstep. This scales the loop to the desired number of calculations.
{
    phir=((phid*adstep)/180.0)*M_PI;        // Phi in radians
    thetad=1.0;                             // Theta in degrees, reset to initial adstep value
    while (thetad<=(180.0/adstep))          // Loop over Theta divided by adstep. Like Phi loop, this scales the loop to the desired number of calculations
    {
        thetar=((thetad*adstep)/180.0)*M_PI;    // Convert theta degrees to radians
        dv = ((maxr*maxr*maxr) * sin(thetar) * arstep * arstep) / 3.0;      // Volume of current segment
        vol += dv;      // Summing all the dv value together to generate a global volume
        thetad+=1.0;    // Increase theta (degrees) by a single step
    }
    phid+=1.0;      // Increase phi (degrees) by a single step
}

vol = vol*1.0; // Volume compensated for any reduction in phi and theta

rad3 = (3.0*vol)/(4.0*M_PI);    // volume equivalent radius^3
rad = pow(rad3,(1.0/3.0));      // volume equivalent radius
vol2 = (4.0/3.0)*M_PI*(maxr*maxr*maxr);     // Calculated volume of a sphere given initial maxr

// Diagnostic output
cout << vol2 << " " << vol << " " << ((vol2-vol)/vol)*100.0 << endl;

}

编辑:将 phid 和 thetad 的起始值更正为 1.0

编辑 2:我只想为 future 的观众更新，使用 Kahan 求和算法 (https://en.wikipedia.org/wiki/Kahan_summation_algorithm) 实际上已经消除了我所有的精度误差，这是由于将一个小数加到一个大数上造成的。还有其他方法，但这是最简单的方法之一，可以完成我需要的工作。对于后代来说，这是取自关于该主题的维基百科页面的示例伪代码:

function KahanSum(input)
 var sum = 0.0
 var c = 0.0                  // A running compensation for lost low-order bits.
 for i = 1 to input.length do
    var y = input[i] - c     // So far, so good: c is zero.
    var t = sum + y          // Alas, sum is big, y small, so low-order digits of y are lost.
    c = (t - sum) - y // (t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
    sum = t           // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
    // Next time around, the lost low part will be added to y in a fresh attempt.
 return sum

Answer 1

就速度而言，我怀疑(没有分析它)大量时间浪费在弧度和度数之间的转换，以及计算所有这些 sin 上。 AFAICT，thetar 在外循环的每次迭代中循环遍历相同的值，因此在主循环之前预先计算一次 sin(thetar) 可能会更有效，并在您的内部循环中进行简单的查找。

至于数值稳定性，随着 vol 变得比 dv 越来越大，您将开始失去越来越多的精度。如果您可以将所有 dv 存储在一个数组中，然后使用分而治之的方法而不是线性传递对其求和，原则上您会获得更好的结果。然后我再次计算(仅)6 480 000 总迭代，所以我认为 double 累加器(持有 15-17 个有效的 base-10 数字)实际上可以处理丢失 6- 7 位数没有太大问题。