Java 的快速超越/三角函数

Question

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

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

对于超过 pi/4 的值，除了硬件之外，Java 还会进行一些昂贵的计算功能。 这样做有充分的理由，但有时您更关心的是速度而不是最后一位的准确性。

Answer 1

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

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

如果您必须推出自己的三角函数，请不要单独使用泰勒级数。 除非你的参数非常小，否则 CORDIC 算法要快得多。 您可以使用 CORDIC 开始，然后使用简短的泰勒级数完善结果。 请参阅此 StackOverflow 问题，了解如何实现三角函数。

Answer 2

在 x86 上，java.lang.Math sin 和 cos 函数不直接调用硬件函数，因为 Intel 在实现它们方面并不总是做得那么好。 bug #4857011 中有一个很好的解释。

https://bugs.java.com/bugdatabase/view_bug?bug_id=4857011

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

“但是评论说罪……”

Answer 3

如果您只需要一些近似值，您可以将 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):

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

Answer 4

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

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

Answer 5

三角函数是查找表的经典示例。 请参阅

• 维基百科上的查找表文章

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

• 定点整数数学库 MathFP

Answer 6

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

在我的实验室电脑上，sin 和 cos 大约需要 144 ns。

Answer 7

在正弦/余弦测试中，我对整数 0 到 100 万进行测试。 我认为 144 ns 对您来说还不够快。

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

您能否在每次操作的时间方面满足您的要求？

Answer 8

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

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

例如，以下是一些可能有用的泰勒级数展开式（取自 wikipedia）：

< img src="https://upload.wikimedia.org/math/d/0/d/d0d68a0e34cc259212ab6e4506e3d99c.png" alt="替代文本">

Answer 9

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

Answer 10

Hart 的计算机近似值。 将一系列不同精度函数的切比雪夫经济近似公式制成表格。

编辑：将我的副本下架，结果是一本不同的书，只是听起来很相似。 这是使用其表格的 sin 函数。 （在 C 中测试，因为这对我来说更方便。）我不知道这是否会比 Java 内置更快，但至少保证它不太准确。 :) 您可能需要首先缩小参数的范围； 请参阅John Cook 的建议。 书中还有反正弦和反正切。

#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;
}

Answer 11

这里是快速逼近三角函数的低级技巧的集合。 有 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;
    }
  }
}

Answer 12

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