【问题标题】:Calculation sine and cosine in one shot一次性计算正弦和余弦
【发布时间】:2015-08-04 16:04:19
【问题描述】:

我有一个科学代码,它使用同一参数的正弦和余弦(我基本上需要那个参数的复指数)。我想知道是否有可能比分别调用正弦和余弦函数更快。

另外,我只需要大约 0.1% 的精度。那么有什么方法可以找到默认的三角函数并截断幂级数以提高速度?

我想到的另一件事是,有没有办法执行余数运算以使结果始终为正?在我自己的算法中,我使用了x=fmod(x,2*pi);,但是如果 x 为负数,我需要添加 2pi(更小的域意味着我可以使用更短的幂级数)

编辑:LUT 原来是最好的方法,但我很高兴我了解了其他近似技术。我还将建议使用明确的中点近似。这就是我最终要做的:

const int N = 10000;//about 3e-4 error for 1000//3e-5 for 10 000//3e-6 for 100 000
double *cs = new double[N];
double *sn = new double[N];
for(int i  =0;i<N;i++){
    double A= (i+0.5)*2*pi/N;
    cs[i]=cos(A);
    sn[i]=sin(A);
}

以下部分近似(中点) sincos(2*pi*(wc2+t[j]*(cotp*t[j]-wc)))

double A=(wc2+t[j]*(cotp*t[j]-wc));
int B =(int)N*(A-floor(A));
re += cs[B]*f[j];
im += sn[B]*f[j];

另一种方法可能是使用切比雪夫分解。您可以使用正交性属性来查找系数。针对指数进行了优化,它看起来像这样:

double fastsin(double x){
    x=x-floor(x/2/pi)*2*pi-pi;//this line can be improved, both inside this 
                              //function and before you input it into the function

    double x2 = x*x;
    return (((0.00015025063885163012*x2- 
   0.008034350857376128)*x2+ 0.1659789684145034)*x2-0.9995812174943602)*x;} //7th order chebyshev approx

【问题讨论】:

标签: c++ algorithm trigonometry


【解决方案1】:

如果您寻求对幂级数具有良好(但不高)准确度的快速评估,您应该在切比雪夫多项式中使用展开式:将系数制成表格(对于 0.1% 的准确度,您需要很少的系数)并使用递归关系评估展开式对于这些多项式(真的很简单)。

参考资料:

  1. 列表系数:http://www.ams.org/mcom/1980-34-149/S0025-5718-1980-0551302-5/S0025-5718-1980-0551302-5.pdf
  2. 切比雪夫扩展评估:https://en.wikipedia.org/wiki/Chebyshev_polynomials

您需要 (a) 在 -pi/2..+pi/2 范围内获取“减少的”参数,然后 (b) 在参数实际应该在时处理结果中的符号整个基本区间的“另一半”-pi..+pi。这些方面不应该构成大问题:

  1. 确定(并“记住”为整数 1 或 -1)原始角度中的符号并继续使用绝对值。
  2. 使用模函数减少到区间 0..2PI
  3. 确定(并“记住”为整数 1 或 -1)它是否在“第二”半部分,如果是,则减去 pi*3/2,否则减去 pi/2。注意:这有效地互换了正弦和余弦(除了符号);在最终评估中考虑到这一点。

这完成了在 -pi/2..+pi/2 中获取角度的步骤 使用 Cheb 扩展评估正弦和余弦后,应用上述步骤 1 和 3 的“标志”以获得值中的正确符号。

【讨论】:

  • 我很高兴了解到这一点。我懒得阅读pdf,但我猜他们使用正交属性分解成ChebyshevU。这就是我所做的i.imgur.com/eSNoMYO.png
  • 另外,切比雪夫多项式可以很好地逼近什么样的函数?它们(对于 sincos)似乎比 pade 近似值要好得多。还是在pdf里有讲?
  • 切比雪夫多项式适用于不是“太狂野”的函数。当然,这是松散的语言。想一想:有渐近线、不连续性等,你就会明白。此外:近似值应适用于有限区间。
【解决方案2】:

只需创建一个查找表。下面将让您查找介于 -2PI 和 2PI 之间的任何弧度值的 sin 和 cos。

// LOOK UP TABLE
var LUT_SIN_COS = [];
var N = 14400;
var HALF_N = N >> 1;
var STEP = 4 * Math.PI / N;
var INV_STEP = 1 / STEP;
// BUILD LUT
for(var i=0, r = -2*Math.PI; i < N; i++, r += STEP) {
    LUT_SIN_COS[2*i] = Math.sin(r);
    LUT_SIN_COS[2*i + 1] = Math.cos(r);
}

您通过以下方式索引查找表:

var index = ((r * INV_STEP) + HALF_N) << 1;
var sin = LUT_SIN_COS[index];
var cos = LUT_SIN_COS[index + 1];

这是一个小提琴,它显示了您可以从不同大小的 LUTS http://jsfiddle.net/77h6tvhj/ 中获得的 % 错误

编辑 这是一个带有 ~benchmark~ 与浮点 sin 和 cos 的 ideone (c++)。 http://ideone.com/SGrFVG 无论 ideone.com 上的基准测试是否值得,LUT 都快 5 倍。

【讨论】:

  • IIRC,对于固定误差,线性插值允许使用 更小的 LUT。 sin 和 cos 是非常好的函数。特别是,如果您可以将较小的 LUT 放入 L1 缓存而不是 L2,它会快很多。当然,您需要两个 LUT 条目,但它们是相邻的,因此通常位于同一缓存行上。
  • @MSalters - 现代处理器上的 L1 和 L2 缓存命中之间的差异约为 2.5 倍(4 个周期对 10 个周期),我认为执行插值的逻辑会弥补差异。此外,如果您真的将 LUT 限制为 -pi/2..pi/2,那么您还将有额外的推断最终符号的开销。但是,如果您手头有代码,请随时将其粘贴到我链接的 ideone 中。
  • 做了一些微基准测试。 0.43 秒 sinf + cosf,0.24 秒查找表(正弦和余弦从 -PI/2 到 +PI/2 的 257 个值,线性插值),0.16 秒 XMScalarSinCos(时间为 10M 计算)。 XMScalarSinCos 使用高次多项式逼近。如您所见,在现代硬件上,浮点乘法和加法甚至比 L1 缓存查找还要快。
  • @Sonts - SIMD(SSE 指令集)令人印象深刻。虽然我不能否认在一般情况下 XM*(或任何其他合适的 SSE 数学库)函数比 LUT 更快,但瓶颈不是 L1 或 L2 缓存查找。在现代 CPU 上,L1 缓存读取约为 4 个周期(L2 约为 10 个周期),XM* 函数大约为 ~~30 个周期;我敢打赌 LUT 代码中的瓶颈是将 fp 值按摩到 int 索引的开销(此代码可能使用较慢的 x87 fp 指令,具体取决于编译器优化设置),以及构成 LUT 的内存可能未对齐。
  • @LouisRicci XMScalarSinCos 不是 SSE 内在函数,它是编译为 x87 指令的常规内联函数。例如,这里是正弦的主要公式: ( ( ( ( (-2.3889859e-08f * y2 + 2.7525562e-06f) * y2 - 0.00019840874f ) * y2 + 0.0083333310f ) * y2 - 0.16666667f ) * y2 + 1.0f ) * y.
【解决方案3】:

一种方法是学习如何实现CORDIC 算法。这在智力上并不难,而且很有趣。这给了你余弦和正弦。 Wikipedia 给出了一个MATLAB example,应该很容易在 C++ 中适应。

请注意,您可以通过降低参数n 来提高速度并降低精度。


关于您的第二个问题,已经有人问过here(在 C 中)。似乎没有简单的方法。

【讨论】:

  • 看起来很有希望。我仍然有兴趣学习(如果可能的话)进行仅给出正值的浮点余数运算。
  • 我不确定我明白你的意思。
  • 当我尝试实现幂级数时,我想使用 2pi 的正弦周期。但是当输入为 fmod(-3pi,2pi) 时,fmod/remainder 操作将给出 -pi。这意味着我尝试实现的幂级数应该在 (-2pi,2pi) 之间有效,或者我应该使用额外的 if 条件将 -pi 变为 +pi。我想知道是否可以让 fmod 始终输出正余数
  • 以上回答。将 MATLAB 代码转换为 C++ 会遇到问题吗?
  • 这是一个仅链接的答案,实际上,因为没有链接,答案毫无价值。
【解决方案4】:

在给定角度和余弦的情况下,您还可以使用平方根计算正弦。

以下示例假设角度范围为 0 到 2π:

 double c = cos(angle);
 double s = sqrt(1.0-c*c);
 if(angle>pi)s=-s;

【讨论】:

    【解决方案5】:

    对于单精度浮点数,Microsoft 对正弦使用 11 度多项式逼近,对余弦使用 10 度逼近:XMScalarSinCos。 他们还有更快的版本 XMScalarSinCosEst,它使用低次多项式。

    如果您不在 Windows 上,您会在 Boost 许可证下找到相同的代码 + 系数 on geometrictools.com

    【讨论】:

      猜你喜欢
      • 2016-06-16
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2016-06-19
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2017-10-06
      相关资源
      最近更新 更多