C: 使用大量 sin() 来提高函数的性能答案

【问题标题】：C: Improving performance of function with heavy sin() usageC: 使用大量 sin() 来提高函数的性能
【发布时间】：2013-12-31 01:22:46
【问题描述】：

我有一个 C 函数，它根据经过的时间计算 4 个正弦的值。使用 gprof，我发现这个函数使用了 100%（准确地说是 100.7%）的 CPU 时间。

void
update_sines(void)
{
    clock_gettime(CLOCK_MONOTONIC, &spec);
    s = spec.tv_sec;
    ms = spec.tv_nsec * 0.0000001;
    etime = concatenate((long)s, ms);

    int k;
    for (k = 0; k < 799; ++k)
    {
        double A1 = 145 * sin((RAND1 * k + etime) * 0.00333) + RAND5;           // Amplitude
        double A2 = 100 * sin((RAND2 * k + etime) * 0.00333) + RAND4;           // Amplitude
        double A3 = 168 * sin((RAND3 * k + etime) * 0.00333) + RAND3;           // Amplitude
        double A4 = 136 * sin((RAND4 * k + etime) * 0.00333) + RAND2;           // Amplitude

        double B1 = 3 + RAND1 + (sin((RAND5 * k) * etime) * 0.00216);           // Period
        double B2 = 3 + RAND2 + (sin((RAND4 * k) * etime) * 0.002);         // Period
        double B3 = 3 + RAND3 + (sin((RAND3 * k) * etime) * 0.00245);           // Period
        double B4 = 3 + RAND4 + (sin((RAND2 * k) * etime) * 0.002);         // Period

        double x = k;                                   // Current x

        double C1 = 0.6 * etime;                            // X axis move
        double C2 = 0.9 * etime;                            // X axis move
        double C3 = 1.2 * etime;                            // X axis move
        double C4 = 0.8 * etime + 200;                          // X axis move

        double D1 = RAND1 + sin(RAND1 * x * 0.00166) * 4;               // Y axis move
        double D2 = RAND2 + sin(RAND2 * x * 0.002) * 4;                 // Y axis move
        double D3 = RAND3 + cos(RAND3 * x * 0.0025) * 4;                // Y axis move
        double D4 = RAND4 + sin(RAND4 * x * 0.002) * 4;                 // Y axis move

        sine1[k] = A1 * sin((B1 * x + C1) * 0.0025) + D1;
        sine2[k] = A2 * sin((B2 * x + C2) * 0.00333) + D2 + 100;
        sine3[k] = A3 * cos((B3 * x + C3) * 0.002) + D3 + 50;
        sine4[k] = A4 * sin((B4 * x + C4) * 0.00333) + D4 + 100;
    }

}

这是 gprof 的输出：

Flat profile:

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls  Ts/call  Ts/call  name    
100.07      0.04     0.04

我目前使用这个获得大约 30-31 fps 的帧速率。现在我想有一种更有效的方法来做到这一点。

正如您所注意到的，我已经将所有除法更改为乘法，但这对性能影响很小。

我怎样才能提高这个数学繁重的函数的性能？

【问题讨论】：

考虑到您在 4 中执行此操作，它可能非常适合转换为 SIMD（例如 SSE 内在函数），特别是如果您不介意使用浮点精度。
在我看来，A、B 和 C 值是完全独立的。为什么不生成 4 个线程，每个 sineX 一个线程？
您将/ 300 替换为* 0.00333 似乎表明您对准确性的期望极低。通过使用sinf（写得很好的sinf 比sin 快）和单精度变量和常量，使整个计算成为单精度。它仍然比0.00333 更准确，是1/300 的精确近似值。
出于兴趣，该函数的目的到底是什么？从维度 POV 来看，sin * sin(sin + sin) + sin 已经够奇怪了，但这完全取决于当前系统时间？那有什么作用？
我不知道为什么我的问题被搁置了？我的代码中有一个瓶颈，因为我经常调用 sin 函数，有什么解决方案可以解决这个问题？我不明白这是多么广泛。但如果是的话，我很抱歉。

标签： c performance math optimization trigonometry

【解决方案1】：

除了其他答案中给出的所有其他建议之外，这是一个纯粹的算法优化。

在大多数情况下，您正在计算sin(k * a + b) 形式的东西，其中a 和b 是常量，k 是循环变量。如果您还要计算cos(k * a + b)，那么您可以使用二维rotation matrix 来形成递归关系（以矩阵形式）：

|cos(k*a + b)| = |cos(a)  -sin(a)| * |cos((k-1)*a + b)|
|sin(k*a + b)|   |sin(a)   cos(a)|   |sin((k-1)*a + b)|

换句话说，您可以根据前一次迭代的值来计算当前迭代的值。因此，您只需要对k == 0 进行完整的三角计算，但其余的可以通过此重复计算（一旦您计算了cos(a) 和sin(a)，它们是常量）。因此，您消除了 75% 的三角函数调用（不清楚是否可以为最后一组三角函数调用相同的技巧）。

【讨论】：

哦，天哪，我在我头上。好的，让我看看我是否明白。你基本上是说，如果我知道 k == 0 的值，我可以算出下一个值，我可以算出 k == 1 的值，而无需再次计算整个方程？
@ReX357：当然。公平地说，对于k == 0，您需要cos 和sin，但是一旦有了这些，您就可以为k == 1 推导出cos 和sin，等等。我太累了，无法在这里进行推导，但请看一下我链接到的表中的前两个方程。将k * a + b 替换为alpha，将-1 * a 替换为beta，然后完成！
@ReX357：事实上，一种更直观的思考方式是，如果你同时拥有sin 和cos 并将它们绘制为二维坐标，你最终会绘制出一个圆圈。要将一个点稍微绕一个圆移动（这是您从迭代到迭代所做的事情），您需要应用旋转。 M 和 N 是适当的 2D rotation matrix 的两个元素。
好的，让我看看我是否做对了。当 k == 0 时，我会同时获得 cos 和 sin 的 A1 的值，并将它们存储在 a1_sin 和 a1_cos 中。现在当 k == 1 时，我可以通过用 M * a1_sin + N * a1_cos 的值交换 sin 来简单地找到 A1 的值？
@ReX357：是的，它类似于a1_sin = M * a1_sin_old + N * a1_cos_old; a1_cos = -N * a1_sin_old + M * a1_cos_old;（我可能把减号放在了错误的地方；从我链接到的旋转矩阵的维基页面上完成它）。

【解决方案2】：

如果您不需要所有这些精度，请为您需要的 sin() 值创建查找，因此如果 1 度就足够了，请使用 double sin_lookup[360], etc.. 如果浮点精度足够，则可能使用 float sin_lookup[360]。

此外，正如 cmets 中所述，在某个点上，根据 Keith，“您还可以考虑在查找值之间使用线性插值，这应该会在相当小的性能成本”

编辑：还可以考虑将硬编码的 A1、A2、A3、A4 模式更改为大小为 [4] 的数组，并从 0 循环到 3 - 应该允许在许多平台上进行矢量化并允许并行化，而无需管理线程

EDIT2：一些代码和结果

（用 C++ 编码只是为了便于在精度之间进行比较，计算在 C 中是相同的）

class simple_trig
{
public:
        simple_trig(size_t prec) : precision(prec)
        {
                static const double PI=3.141592653589793;
                const double dprec=(double)prec;
                const double quotient=(2.0*PI)/dprec;
                rev_quotient=dprec/(2.0*PI);
                values.reserve(prec);

                for (int i=0; i < precision; ++i)
                {
                        values[i]=::sin(quotient*(double)i);
                }
        }

        double sin(double x) const
        {
                double cvt=x*rev_quotient;
                int index=(int)cvt;
                double delta=cvt-(double)index;
                int lookup1=index%precision;
                int lookup2=(index+1)%precision;
                return values[lookup1]*(1.0-delta)+values[lookup2]*delta;
        }

        double cos(double x) const
        {
                double cvt=x*rev_quotient;
                int index=(int)cvt;
                double delta=cvt-(double)index;
                int lookup1=(index+precision/4)%precision;
                int lookup2=(index+precision/4+1)%precision;
                return values[lookup1]*(1.0-delta)+values[lookup2]*delta;
        }

private:
        const size_t precision;
        double rev_quotient;
        std::vector<double> values;
};

示例低为 100，中为 1000，高为 10,000

X=0 Sin=0 Sin Low=0 Sin Med=0 Sin High=0
X=0 Cos=1 Cos Low=1 Cos Med=1 Cos High=1
X=0.5 Sin=0.479426 Sin Low=0.479389 Sin Med=0.479423 Sin High=0.479426
X=0.5 Cos=0.877583 Cos Low=0.877512 Cos Med=0.877578 Cos High=0.877583
X=1.33333 Sin=0.971938 Sin Low=0.971607 Sin Med=0.971935 Sin High=0.971938
X=1.33333 Cos=0.235238 Cos Low=0.235162 Cos Med=0.235237 Cos High=0.235238
X=2.25 Sin=0.778073 Sin Low=0.777834 Sin Med=0.778072 Sin High=0.778073
X=2.25 Cos=-0.628174 Cos Low=-0.627986 Cos Med=-0.628173 Cos High=-0.628174
X=3.2 Sin=-0.0583741 Sin Low=-0.0583689 Sin Med=-0.0583739 Sin High=-0.0583741
X=3.2 Cos=-0.998295 Cos Low=-0.998166 Cos Med=-0.998291 Cos High=-0.998295
X=4.16667 Sin=-0.854753 Sin Low=-0.854387 Sin Med=-0.854751 Sin High=-0.854753
X=4.16667 Cos=-0.519036 Cos Low=-0.518818 Cos Med=-0.519034 Cos High=-0.519036
X=5.14286 Sin=-0.90877 Sin Low=-0.908542 Sin Med=-0.908766 Sin High=-0.90877
X=5.14286 Cos=0.417296 Cos Low=0.417195 Cos Med=0.417294 Cos High=0.417296
X=6.125 Sin=-0.157526 Sin Low=-0.157449 Sin Med=-0.157526 Sin High=-0.157526
X=6.125 Cos=0.987515 Cos Low=0.987028 Cos Med=0.987512 Cos High=0.987515
X=7.11111 Sin=0.73653 Sin Low=0.736316 Sin Med=0.736527 Sin High=0.73653
X=7.11111 Cos=0.676405 Cos Low=0.676213 Cos Med=0.676403 Cos High=0.676405
X=8.1 Sin=0.96989 Sin Low=0.969741 Sin Med=0.969887 Sin High=0.96989
X=8.1 Cos=-0.243544 Cos Low=-0.24351 Cos Med=-0.243544 Cos High=-0.243544
X=9.09091 Sin=0.327701 Sin Low=0.327558 Sin Med=0.3277 Sin High=0.327701
X=9.09091 Cos=-0.944782 Cos Low=-0.944381 Cos Med=-0.944779 Cos High=-0.944782
X=10.0833 Sin=-0.611975 Sin Low=-0.611673 Sin Med=-0.611973 Sin High=-0.611975
X=10.0833 Cos=-0.790877 Cos Low=-0.790488 Cos Med=-0.790875 Cos High=-0.790877

【讨论】：

这假设调用sin() 是瓶颈（这似乎很可能，但值得检查）。您还可以考虑在查找值之间使用线性插值，这应该以相当小的性能成本为您提供更好的准确性（合理的连续函数而不是阶跃函数）。
@KeithThompson 好主意，这取决于需要多少精度和分配多少空间，可能不需要度数，甚至 1/10 度，但除此之外它很可能会得到回报
并且取决于应用程序（我不假装知道 OP 在做什么），使您的伪sin 函数或多或少是连续的而不是阶跃函数可能是一个非常好的事情，不仅仅是获得几乎正确的值。
@KeithThompson 已经有一段时间了，但是否也有某个点可以使用 x 代替 sin(x)
不是确定点，但是是的，lim(n -> 0) sin(x)/x = 1（我希望我的 ersatz 数学符号足够清楚）。

【解决方案3】：

在我看来，sine1、sine2、sine3 和 sine4 数组是完全独立的。所以你基本上是为 4 个没有依赖关系的不同数组运行一个 for 循环。

产生 4 个线程，每个线程 1 个，因此您有 4 个 for 循环同时运行。在多核机器上，这应该会显着加快您的功能。事实上，它应该是完美的 4 倍加速 (+- ...)。

【讨论】：

【解决方案4】：

实际上将线程的使用（考虑一下 OpenMP）和表的使用结合起来是一个好主意。如果可能，请使用 float 而不是 double，并且根据平台，您还可以使用 simd 指令，但后者会使线程的使用变得不必要。

干杯

【讨论】：

OpenMP 在这里可能会有所帮助，但我担心每次调用此函数时都会产生设置开销。据推测，单次调用此函数的运行时间非常短。
@OliCharlesworth 你是对的，但我的建议是线程；) OpenMP 只是在使用线程被认为太难的情况下。
@OliCharlesworth - 它在 4 个变量集上循环 800 次
@GlennTeitelbaum：确实，但我希望最多需要几毫秒？将其加快 4 倍（或更多，如果您开始加入此处建议的其他一些优化），并且线程启动时间将很重要。我想你可以准备好一个线程池，但现在事情开始变得复杂了;)

【解决方案5】：

这是使用已接受答案中建议的旋转矩阵的 C++ sn-p。

   float a = 0.343;
   float b = 2.3232;
   float sina{};
   float cosa{};
   sincosf(a, &sina, &cosa);
   float resSin{};
   float resCos{};

   for (int k = 0; k < 5; k++) {
     if (k == 0) {
       sincosf(b, &resSin, &resCos);
     } else {
       float newResCos, newResSin;
       newResCos = cosa * resCos - sina * resSin;
       newResSin = sina * resCos + cosa * resSin;
       resCos = newResCos;
       resSin = newResSin;
     }
   }

【讨论】：