我怎样才能优化这个计算？(x ^ a + y ^ a + z ^ a)^(1/a)

Question

正如标题所示.我需要做这样的大量计算:

re = (x^a + y^a + z^a)^(1/a).

其中{ x,y,z }> = 0.更具体地,a是正浮点常数,x,y,z是浮点数.这^是一个取幂运算符.

目前,我不想使用SIMD,但希望有其他一些技巧来加速它.

static void heavy_load(void) {
  static struct xyz_t {
    float x,y,z;
  };
  struct xyz_t xyzs[10000];
  float re[10000] = {.0f};
  const float a = 0.2;

  /* here fill xyzs using some random positive floating point values */

  for (i = 0; i < 10000; ++ i) {
    /* here is what we need to optimize */
    re[i] = pow((pow(xyzs[i].x, a) + pow(xyzs[i].y, a) + pow(xyzs[i].z, a)), 1.0/a);
  }
}

Answer 1

首先要做的是通过分解其中一个术语来摆脱其中一个取幂.以下代码使用

(x^a + y^a + z^a)^(1/a) = x * ((1 + (y/x)^a + (z/x)^a)^(1/a))

考虑到三者中最大的一个将更安全一点,也许更准确.

另一个优化是利用a为0.1或0.2的事实.使用Chebychev多项式近似来近似x ^ a.下面的代码只有x ^ 0.1的近似值; x ^ 0.2就是那个平方.最后,由于1/a是一个小整数(5或10),因此可以用少量乘法替换最终取幂.

要查看函数中发生了什么powtenth并powtenthnorm看到这个答案:使用const非整数指数优化pow()？.

#include <stdlib.h>
#include <math.h>


float powfive (float x);
float powtenth (float x);
float powtenthnorm (float x);

// Returns (x^0.2 + y^0.2 + z^0.2)^5
float pnormfifth (float x, float y, float z) {
   float u = powtenth(y/x);
   float v = powtenth(z/x);
   return x * powfive (1.0f + u*u + v*v);
}

// Returns (x^0.1 + y^0.1 + z^0.1)^10
float pnormtenth (float x, float y, float z) {
   float u = powtenth(y/x);
   float v = powtenth(z/x);
   float w = powfive (1.0f + u + v);
   return x * w * w;
}

// Returns x^5
float powfive (float x) {
   float xsq = x*x;
   return xsq*xsq*x;
}

// Returns x^0.1.
float powtenth (float x) {
   static const float pow2_tenth[10] = {
      1.0,
      pow(2.0, 0.1),
      pow(4.0, 0.1),
      pow(8.0, 0.1),
      pow(16.0, 0.1),
      pow(32.0, 0.1),
      pow(64.0, 0.1),
      pow(128.0, 0.1),
      pow(256.0, 0.1),
      pow(512.0, 0.1)
   };

   float s;
   int iexp;

   s = frexpf (x, &iexp);
   s *= 2.0;
   iexp -= 1;

   div_t qr = div (iexp, 10);
   if (qr.rem < 0) {
      qr.quot -= 1;
      qr.rem += 10;
   }

   return ldexpf (powtenthnorm(s)*pow2_tenth[qr.rem], qr.quot);
}

// Returns x^0.1 for x in [1,2), to within 1.2e-7 (relative error).
// Want more precision? Add more Chebychev polynomial coefs.
float powtenthnorm (float x) {
   static const int N = 8;

   // Chebychev polynomial terms.
   // Non-zero terms calculated via
   //   integrate (2/pi)*ChebyshevT[n,u]/sqrt(1-u^2)*((u+3)/2)^0.1
   //   from -1 to 1
   // Zeroth term is similar except it uses 1/pi rather than 2/pi.
   static const float Cn[N] = {
       1.0386703502389972,
       3.55833786872637476e-2,
      -2.7458105122604368629e-3,
       2.9828558990819401155e-4,
      -3.70977182883591745389e-5,
       4.96412322412154169407e-6,
      -6.9550743747592849e-7,
       1.00572368333558264698e-7};

   float Tn[N];

   float u = 2.0*x - 3.0;


   Tn[0] = 1.0;
   Tn[1] = u;
   for (int ii = 2; ii < N; ++ii) {
      Tn[ii] = 2*u*Tn[ii-1] - Tn[ii-2];
   }

   float y = 0.0;
   for (int ii = N-1; ii > 0; --ii) {
      y += Cn[ii]*Tn[ii];
   }

   return y + Cn[0];
}

Answer 2

一些 C 优化——不多，但有所作为。

[编辑] 抱歉 - 回到烦恼：如果 FP 值变化很大，那么经常发生re = a（或 b 或 c）。对大幅度差异进行测试将无需调用pow()某些 x、y 或 z。这有助于平均时间，但对最坏情况时间没有帮助。

替换1.0/a为a_inverse循环之前设置的which。

替换pow()为powf()，否则您将调用double的版本pow()。

float re[10000] = { 0.0f}次要：除了刷新内存缓存之外，不需要初始化。

次要：使用指针而不是索引数组可能会节省一点。

次要：对于某些平台，使用 typedouble 可能会运行得更快。- 推翻我的上述pow()/powf()评论。

次要：尝试 3 个独立的 x、y、z 数组。每个float *。

显然，分析很有帮助，但我们假设这是众所周知的。