java - Java の高速な超越/三角関数

Question

java.lang.Math の三角関数は非常に遅いため、迅速で適切な近似を行うライブラリはありますか? あまり精度を落とさずに数倍速く計算できるようです。(私のマシンでは、乗算には 1.5ns、java.lang.Math.sin では 46ns から 116ns かかります)。残念ながら、ハードウェア機能を使用する方法はまだありません。

更新: 関数は、たとえば GPS 計算には十分に正確である必要があります。つまり、少なくとも 7 桁の精度が必要であり、単純なルックアップテーブルは除外されます。基本的な x86 システムでは、java.lang.Math.sin よりもはるかに高速です。そうでなければ意味がありません。

pi/4 を超える値の場合、Java はハードウェア関数に加えてコストのかかる計算を行います。これには正当な理由がありますが、最終ビットの精度よりも速度を重視する場合があります。

score 15 · Accepted Answer

ハートによるコンピュータ近似。さまざまな精度で一連の関数のチェビシェフ節約近似式を表にします。

編集：私のコピーを棚から取り出したところ、非常に似たように聞こえる別の本であることが判明しました. テーブルを使用した 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;
}

score 13 · Accepted Answer

これは、三角関数をすばやく近似するための低レベルのトリックのコレクションです。C のサンプルコードは理解するのが難しいと思いますが、その手法は Java で簡単に実装できます。

これは、Java での invsqrt と atan2 の同等の実装です。

他の trig 関数についても同様のことを行うことができましたが、プロファイリングで 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;
    }
  }
}

score 5 · Accepted Answer

5

それはそれを作るかもしれません：http ：//sourceforge.net/projects/jafama/

于 2010-07-02T23:24:52.943 に答える

score 4 · Accepted Answer

x86 では、java.lang.Math の sin および cos 関数は、ハードウェア関数を直接呼び出すことはありません。これは、Intel がそれらを適切に実装することを常に行っているとは限らないためです。バグ #4857011 にわかりやすい説明があります。

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

不正確な結果についてよく考えてみてください。他のコードでこれを見つけるのにどれだけの時間を費やしているかは面白いです。

「でもコメントに罪って書いてある…」

score 1 · Accepted Answer

三角関数は、ルックアップテーブルの古典的な例です。優れたものを見る

ウィキペディアのルックアップテーブルの記事

J2ME のライブラリを検索している場合は、次を試すことができます。

固定小数点整数数学ライブラリ MathFP

score 0 · Accepted Answer

sin/cos テストでは、0 から 100 万の整数に対して実行していました。144 ns は十分に高速ではないと思います。

必要な速度に対する特定の要件はありますか?

1 回の操作にかかる時間に関して、要件を十分に満たすことができますか?

score 0 · Accepted Answer

java.lang.Math 関数は、ハードウェア関数を呼び出します。作成できる簡単な承認があるはずですが、それほど正確ではありません。

私のラボトップでは、sin と cos は約 144 ns かかります。

score -1 · Accepted Answer

これらのルーチンが遅すぎる場合、何をする必要があるか詳しく説明していただけますか? 何らかの方法で事前に座標変換を行うことができる場合があります。

java - Java の高速な超越/三角関数

12 に答える 12

Related

Reference