Java的快速超越/三角函数答案

【问题标题】：Fast transcendent / trigonometric functions for JavaJava的快速超越/三角函数
【发布时间】：2010-07-03 07:29:29
【问题描述】：

由于 java.lang.Math 中的三角函数非常慢：是否有一个库可以快速且良好地进行近似？在不损失太多精度的情况下，似乎可以将计算速度提高几倍。（在我的机器上，乘法需要 1.5ns，java.lang.Math.sin 需要 46ns 到 116ns）。不幸的是，目前还没有使用硬件功能的方法。

更新：函数应该足够准确，例如 GPS 计算。这意味着您需要至少 7 个十进制数字的准确性，这排除了简单的查找表。它应该比基本 x86 系统上的 java.lang.Math.sin 快得多。否则就没有意义了。

对于超过 pi/4 的值，除了硬件功能之外，Java 还执行 some expensive computations。这样做是有充分理由的，但有时您更关心速度而不是最后一点的准确性。

【问题讨论】：

有多快，多好？您总是可以只使用泰勒级数的前几个术语...这非常快速且随心所欲。
永远不要使用泰勒级数。看我的回答stackoverflow.com/questions/345085/…
完全错误，每个近似值都有时间和地点。当然，我通常不会使用超过两个或三个泰勒近似项，但对于正弦、余弦和指数，它们会很好地收敛。参加一两门数值分析课程，您可能会学到一些东西。
我支持我的评论。永远不要使用泰勒级数。它们针对非常接近单点的功能评估进行了优化。 x - x^3/6 在大约 pi/4 的值处开始丢失，即使这样，准确性也很粗糙。最小二乘拟合很容易做到。
泰勒级数可以惊人地准确：Taylor series approximations, illustrated.

标签： java optimization math trigonometry

【解决方案1】：

Computer Approximations 来自哈特。将Chebyshev-economized 一系列不同精度的函数的近似公式制成表格。

编辑： 将我的副本下架，结果是a different book，听起来非常相似。这是一个使用其表的 sin 函数。（在 C 中测试，因为这对我来说更方便。）我不知道这是否会比 Java 内置更快，但至少可以保证它不那么准确。 :) 您可能需要先对参数进行范围缩小；见John Cook's suggestions。本书还有arcsin和arctan。

#include <math.h>
#include <stdio.h>

// Return an approx to sin(pi/2 * x) where -1 <= x <= 1.
// In that range it has a max absolute error of 5e-9
// according to Hastings, Approximations For Digital Computers.
static double xsin (double x) {
  double x2 = x * x;
  return ((((.00015148419 * x2
             - .00467376557) * x2
            + .07968967928) * x2
           - .64596371106) * x2
          + 1.57079631847) * x;
}

int main () {
  double pi = 4 * atan (1);
  printf ("%.10f\n", xsin (0.77));
  printf ("%.10f\n", sin (0.77 * (pi/2)));
  return 0;
}

【讨论】：

看起来像是某种形式的泰勒级数近似。在此处查看通用公式：en.wikipedia.org/wiki/…

【解决方案2】：

Here 是一组用于快速逼近三角函数的低级技巧。有一些 C 语言的示例代码我觉得很难理解，但这些技术在 Java 中同样容易实现。

这是我在 Java 中对 invsqrt 和 atan2 的等效实现。

我可以为其他三角函数做类似的事情，但我发现没有必要，因为分析表明只有 sqrt 和 atan/atan2 是主要瓶颈。

public class FastTrig
{
  /** Fast approximation of 1.0 / sqrt(x).
   * See <a href="http://www.beyond3d.com/content/articles/8/">http://www.beyond3d.com/content/articles/8/</a>
   * @param x Positive value to estimate inverse of square root of
   * @return Approximately 1.0 / sqrt(x)
   **/
  public static double
  invSqrt(double x)
  {
    double xhalf = 0.5 * x; 
    long i = Double.doubleToRawLongBits(x);
    i = 0x5FE6EB50C7B537AAL - (i>>1); 
    x = Double.longBitsToDouble(i);
    x = x * (1.5 - xhalf*x*x); 
    return x; 
  }

  /** Approximation of arctangent.
   *  Slightly faster and substantially less accurate than
   *  {@link Math#atan2(double, double)}.
   **/
  public static double fast_atan2(double y, double x)
  {
    double d2 = x*x + y*y;

    // Bail out if d2 is NaN, zero or subnormal
    if (Double.isNaN(d2) ||
        (Double.doubleToRawLongBits(d2) < 0x10000000000000L))
    {
      return Double.NaN;
    }

    // Normalise such that 0.0 <= y <= x
    boolean negY = y < 0.0;
    if (negY) {y = -y;}
    boolean negX = x < 0.0;
    if (negX) {x = -x;}
    boolean steep = y > x;
    if (steep)
    {
      double t = x;
      x = y;
      y = t;
    }

    // Scale to unit circle (0.0 <= y <= x <= 1.0)
    double rinv = invSqrt(d2); // rinv ≅ 1.0 / hypot(x, y)
    x *= rinv; // x ≅ cos θ
    y *= rinv; // y ≅ sin θ, hence θ ≅ asin y

    // Hack: we want: ind = floor(y * 256)
    // We deliberately force truncation by adding floating-point numbers whose
    // exponents differ greatly.  The FPU will right-shift y to match exponents,
    // dropping all but the first 9 significant bits, which become the 9 LSBs
    // of the resulting mantissa.
    // Inspired by a similar piece of C code at
    // http://www.shellandslate.com/computermath101.html
    double yp = FRAC_BIAS + y;
    int ind = (int) Double.doubleToRawLongBits(yp);

    // Find φ (a first approximation of θ) from the LUT
    double φ = ASIN_TAB[ind];
    double cφ = COS_TAB[ind]; // cos(φ)

    // sin(φ) == ind / 256.0
    // Note that sφ is truncated, hence not identical to y.
    double sφ = yp - FRAC_BIAS;
    double sd = y * cφ - x * sφ; // sin(θ-φ) ≡ sinθ cosφ - cosθ sinφ

    // asin(sd) ≅ sd + ⅙sd³ (from first 2 terms of Maclaurin series)
    double d = (6.0 + sd * sd) * sd * ONE_SIXTH;
    double θ = φ + d;

    // Translate back to correct octant
    if (steep) { θ = Math.PI * 0.5 - θ; }
    if (negX) { θ = Math.PI - θ; }
    if (negY) { θ = -θ; }

    return θ;
  }

  private static final double ONE_SIXTH = 1.0 / 6.0;
  private static final int FRAC_EXP = 8; // LUT precision == 2 ** -8 == 1/256
  private static final int LUT_SIZE = (1 << FRAC_EXP) + 1;
  private static final double FRAC_BIAS =
    Double.longBitsToDouble((0x433L - FRAC_EXP) << 52);
  private static final double[] ASIN_TAB = new double[LUT_SIZE];
  private static final double[] COS_TAB = new double[LUT_SIZE];

  static
  {
    /* Populate trig tables */
    for (int ind = 0; ind < LUT_SIZE; ++ ind)
    {
      double v = ind / (double) (1 << FRAC_EXP);
      double asinv = Math.asin(v);
      COS_TAB[ind] = Math.cos(asinv);
      ASIN_TAB[ind] = asinv;
    }
  }
}

【讨论】：

我不会因为变量名称而投反对票，但是您想维护带有 משתנה 和 פונקציה 等标识符的代码吗？
@dotancohen：我的变量名是这样呈现给你的吗？我用 UTF-8 发布它们。听起来您的浏览器猜错了编码（CP1255？）
不，我的浏览器正确呈现希腊语。但是，如果希腊语标识符是公平的游戏，那么为什么不是希伯来语，甚至韩语呢？我试图说明，虽然我们的工具可能会被滥用，但我们真的不应该为了那些追随我们的人而滥用它们。即使在“仅供内部使用”的代码中，您的儿子或本科生也可能继承它！是的，我三岁的孩子今天说“public static void main”！
@dotancohen，因为我只使用那些希腊字母来表示已确定的含义（即，大多数讲英语的人都知道一些数学知识。）
喜欢 unicode 变量名！它们实际上使数学代码更清晰！

【解决方案3】：

这可能会成功：http://sourceforge.net/projects/jafama/

【讨论】：

【解决方案4】：

我很惊讶内置的 Java 函数会这么慢。当然，JVM 调用的是 CPU 上的本机三角函数，而不是在 Java 中实现算法。你确定你的瓶颈是调用触发函数而不是一些周围的代码吗？也许一些内存分配？

你能用 C++ 重写你的代码中进行数学运算的部分吗？仅仅调用 C++ 代码来计算三角函数可能不会加快速度，但是将一些上下文（例如外循环）移动到 C++ 可能会加快速度。

如果您必须滚动自己的三角函数，请不要单独使用泰勒级数。除非您的论点非常小，否则 CORDIC 算法要快得多。您可以使用 CORDIC 开始，然后使用简短的泰勒级数来完善结果。在 how to implement trig functions 上查看这个 StackOverflow 问题。

【讨论】：

【解决方案5】：

在 x86 上，java.lang.Math sin 和 cos 函数不直接调用硬件函数，因为 Intel 并不总是能很好地实现它们。在错误 #4857011 中有一个很好的解释。

http://bugs.sun.com/bugdatabase/view_bug.do?bug_id=4857011

您可能需要认真考虑一个不精确的结果。有趣的是，我经常花时间在其他代码中找到它。

“但是评论说罪...”

【讨论】：

我不明白您在链接后的评论。你在说什么样的错误？
特别喜欢“almabench 代码和这个bug 中提交的代码都没有真正检查结果以验证它们是否合理”。在结束消息中发表评论。
如果代码想要计算以弧度以外的方式测量的角度的正弦或余弦（可能 99% 的时间都涉及正弦和余弦），如果 @987654322 准确度将是最佳的@ 和 cos 函数假设 pi 的值是在角度到弧度转换中使用的值。该文件所指的“sin(x) 中的不准确性提高了Math.sin(x*(2.0*Math.Pi)) 计算 sin(2πx) 的准确性。同样，假设的改进实际上使 sin(2πx) 的评估通常不太准确。事实上，对于某些 x 值, ...
...使用 mod(Math.PI * 2.0) 减少将允许公式很好地计算 sin(2πx)，即使 x 很大，误差分布在零的任一侧。 “改进的”减少并没有消除乘法中的舍入误差，而是增加了额外的系统偏差，该偏差将角度低估了大约 3E-17 倍。

【解决方案6】：

如果您只需要一些近似值，您可以将 sin 和 cos 预先存储在一个数组中。例如，如果要存储 0° 到 360° 的值：

double sin[]=new double[360];
for(int i=0;i< sin.length;++i) sin[i]=Math.sin(i/180.0*Math.PI):

然后您可以使用度数/整数而不是弧度/双精度来使用此数组。

【讨论】：

是的，但这是非常不准确的。我在想一些更好的东西，比如多项式插值。
提醒人们像厄运一样的旧的预计算天......无论如何，您可以通过不只生成 360 个值来提高准确性，但例如0xffff 值。
@hstoerr：为什么不准确？它与数组的长度（即角度的粒度）一样精确。这是速度和内存之间的良好平衡，而性能在这里是最佳的。
如果您想要 7 位小数的精度，就像 GPS 计算所需要的那样，您需要 10000000 个值。您可能不想预先计算那么多，是吗？

【解决方案7】：

我还没有听说过任何库，可能是因为很少能看到触发繁重的 Java 应用程序。使用 JNI（精度相同，性能更好）、数值方法（可变精度/性能）或简单的近似表也很容易实现。

与任何优化一样，最好在重新发明轮子之前测试这些函数实际上是一个瓶颈。

【讨论】：

使用 JNI 进行单个 Math.sin 调用可能由于开销而无法工作。也许如果你将更多的程序放在 C 中，但你可以从 C 开始。
几年前面临类似的问题，调用空函数的 JNI 开销比调用 Math.sin() 慢。那是 1.3 或 1.4 的版本，所以它可能已经改变了，但是 afaik 现在并没有太大的不同。

【解决方案8】：

三角函数是查找表的经典示例。看到优秀的

Lookup table article at wikipedia

如果您正在搜索 J2ME 库，您可以尝试：

定点整数数学库MathFP

【讨论】：

【解决方案9】：

java.lang.Math 函数调用硬件函数。您应该可以做出一些简单的认可，但它们不会那么准确。

在我的 labtop 上，sin 和 cos 大约需要 144 ns。

【讨论】：

据我所知他们不使用硬件功能。 Math.sin 的 Javadoc 表示，结果必须精确到倒数第二位，而硬件实现可能无法满足这一要求。所以它在软件中。
我在我的系统上试过 - 2ns 用于乘法，46ns 用于 Math.sin。这不可能是硬件——罪并没有那么慢。
是的，可以。在 x87 FPU 上，乘法大约为 4 个周期，正弦在 100 范围内。因此，该结果与在硬件中评估它们的 2GHz 处理器完全一致。
好的，我必须用 C++ 或其他东西检查同样的事情。仍然：计算的时间取决于论点。如果你计算 0.1 的 sin 需要 46ns，如果你计算 6.28 的 sin 需要 115ns。那不是硬件，不是吗？
hstoerr：Bruce ONeel 引用的错误详细说明了为什么更大的参数会导致更长的计算时间。基本上，英特尔的 sin/cos 实现通过gardenhoses 为 [-pi/4,pi/4] 之外的参数吸高尔夫球，并且 JVM 必须手动将参数映射到这个范围内。

【解决方案10】：

在 sin/cos 测试中，我对 0 到 100 万的整数进行了测试。我假设 144 ns 对你来说不够快。

您对所需的速度有具体要求吗？

您能否以令人满意的每次操作时间来确定您的要求？

【讨论】：

【解决方案11】：

如果您想使用现有的东西，请查看Apache Commons Math package。

如果性能真的至关重要，那么您可以使用标准数学方法自行实现这些函数 - 特别是泰勒/麦克劳林级数。

例如，这里有几个可能有用的泰勒级数展开式（取自wikipedia）：

【讨论】：

commons math 是一个非常好的技巧，但我没有找到任何更快的替代 Math.sin 的方法，例如。有吗？
如果您想获得合理的准确性，泰勒级数可能不会更快。你必须做一些更聪明的事情，比如使用分段多项式。
在应用泰勒级数之前使用 CORDIC 减少参数。见stackoverflow.com/questions/345085/…
永远不要使用泰勒级数来近似函数。看我的评论stackoverflow.com/questions/345085/…

【解决方案12】：

如果这些例程太慢，您能否详细说明您需要做什么。您可能能够以某种方式提前进行一些坐标转换。

【讨论】：