java - 了解 strictMath java 库

Question

我感到无聊并决定在不引用任何 Math.java 函数的情况下重新制作平方根函数。我已经到了这一点:

package sqrt;
public class SquareRoot {

public static void main(String[] args) {
    System.out.println(sqrtOf(8));

}

public static double sqrtOf(double n){
    double x = log(n,2);
    return powerOf(2, x/2);
}

public static double log(double n, double base)
{
    return (Math.log(n)/Math.log(base));
}

public static double powerOf(double x, double y) {
    return powerOf(e(),y *  log(x, e()));
}   

public static int factorial(int n){
    if(n <= 1){
        return 1;
    }else{
        return n * factorial((n-1));
    }
}

public static double e(){
    return 1/factorial(1);
}
public static double e(int precision){
    return 1/factorial(precision);
}

}

您可能很清楚地看到，我在我的 powerOf() 函数中谈到了无限记忆自己的地步。我可以替换它并使用 Math.exp(y * log(x, e())，因此我深入研究了 Math 源代码以查看它如何处理我的问题，从而使我陷入困境。

public static double exp(double a) {
     return StrictMath.exp(a); // default impl. delegates to StrictMath
 }

这导致:

public static double exp(double x)
{
 if (x != x)
   return x;
 if (x > EXP_LIMIT_H)
   return Double.POSITIVE_INFINITY;
 if (x < EXP_LIMIT_L)
   return 0;

 // Argument reduction.
 double hi;
 double lo;
 int k;
 double t = abs(x);
 if (t > 0.5 * LN2)
   {
     if (t < 1.5 * LN2)
       {
         hi = t - LN2_H;
         lo = LN2_L;
         k = 1;
       }
     else
       {
         k = (int) (INV_LN2 * t + 0.5);
         hi = t - k * LN2_H;
         lo = k * LN2_L;
       }
     if (x < 0)
       {
         hi = -hi;
         lo = -lo;
         k = -k;
       }
     x = hi - lo;
   }
 else if (t < 1 / TWO_28)
   return 1;
 else
   lo = hi = k = 0;

// Now x is in primary range.
 t = x * x;
 double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
 if (k == 0)
   return 1 - (x * c / (c - 2) - x);
 double y = 1 - (lo - x * c / (2 - c) - hi);
 return scale(y, k);

引用的值:

 LN2 = 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL.
 LN2_H = 0.6931471803691238, // Long bits 0x3fe62e42fee00000L.
 LN2_L = 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L.
 INV_LN2 = 1.4426950408889634, // Long bits 0x3ff71547652b82feL.
 INV_LN2_H = 1.4426950216293335, // Long bits 0x3ff7154760000000L.
 INV_LN2_L = 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L.
 P1 = 0.16666666666666602, // Long bits 0x3fc555555555553eL.
 P2 = -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L.
 P3 = 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL.
 P4 = -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L.
 P5 = 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L.
 TWO_28 = 0x10000000, // Long bits 0x41b0000000000000L

这是我开始迷路的地方。但我可以做一些假设，到目前为止答案开始变得估计。然后我发现自己在这里:

private static double scale(double x, int n)
{
  if (Configuration.DEBUG && abs(n) >= 2048)
    throw new InternalError("Assertion failure");
  if (x == 0 || x == Double.NEGATIVE_INFINITY
      || ! (x < Double.POSITIVE_INFINITY) || n == 0)
    return x;
  long bits = Double.doubleToLongBits(x);
  int exp = (int) (bits >> 52) & 0x7ff;
  if (exp == 0) // Subnormal x.
    {
      x *= TWO_54;
      exp = ((int) (Double.doubleToLongBits(x) >> 52) & 0x7ff) - 54;
    }
  exp += n;
  if (exp > 0x7fe) // Overflow.
    return Double.POSITIVE_INFINITY * x;
  if (exp > 0) // Normal.
    return Double.longBitsToDouble((bits & 0x800fffffffffffffL)
                                   | ((long) exp << 52));
  if (exp <= -54)
    return 0 * x; // Underflow.
  exp += 54; // Subnormal result.
  x = Double.longBitsToDouble((bits & 0x800fffffffffffffL)
                              | ((long) exp << 52));
  return x * (1 / TWO_54);
}


 TWO_54 = 0x40000000000000L

我想说，虽然我非常了解数学和编程，但我发现自己处于两者的弗兰肯斯坦怪物组合的地步。我注意到了向位的内在转换(我对此几乎没有经验)，我希望有人可以向我解释可以这么说的“幕后”发生的过程。具体来说，我迷路的地方来自 wards 的 exp() 方法中的“现在 x 处于主要范围内”，以及所引用的值真正代表什么。我要求有人帮助我不仅理解方法本身，而且帮助我理解它们是如何得出答案的。根据需要随意深入。

编辑:如果有人可以制作这个标签:“strictMath”，那就太好了。我相信它的大小和派生自它的数学库证明了它的存在。

Answer 1

对指数函数:

会发生什么

exp(x) = 2^k * exp(x-k*log(2))

被利用为正 x .一些魔法被用来为大 x 获得更一致的结果。减少的地方x-k*log(2)将引入取消错误。

关于减少x在区间 0.5..1.5 上具有最小最大误差的有理逼近被使用，参见 Pade 近似和类似的。这是基于对称公式

exp(x) = exp(x/2)/exp(-x/2) = (c(x²)+x)/(c(x²)-x)

(注意代码中的 c 是 x+c(x)-2 )。使用泰勒级数时，c(x*x)=x*coth(x/2) 的近似值是基于

c(u)=2 + 1/6*u - 1/360*u^2 + 1/15120*u^3 - 1/604800*u^4 + 1/23950080*u^5 - 691/653837184000*u^6

scale(x,n)函数实现乘法 x*2^n通过直接操作 double 的位组件中的指数浮点格式。

---

计算平方根

要计算平方根，直接计算它们会更有利。首先通过

减少近似参数的区间

sqrt(x)=2^k*sqrt(x/4^k)

这可以通过直接操作 double 的位格式再次有效地完成。 .

在 x 之后减少到区间 0.5..2.0然后可以使用以下形式的公式

u = (x-1)/(x+1)

y = (c(u*u)+u) / (c(u*u)-u)

基于

sqrt(x)=sqrt(1+u)/sqrt(1-u)

和

c(v) = 1+sqrt(1-v) = 2 - 1/2*v - 1/8*v^2 - 1/16*v^3 - 5/128*v^4 - 7/256*v^5 - 21/1024*v^6 - 33/2048*v^7 - ...

在没有位操作的程序中，这看起来像

double my_sqrt(double x) {
    double c,u,v,y,scale=1;
    int k=0;
    if(x<0) return NaN;
    while(x>2  ) { x/=4; scale *=2; k++; }
    while(x<0.5) { x*=4; scale /=2; k--; }
    // rational approximation of sqrt
    u = (x-1)/(x+1); 
    v = u*u;
    c = 2 - v/2*(1 + v/4*(1 + v/2));
    y = 1 + 2*u/(c-u); // = (c+u)/(c-u);
    // one Halley iteration
    y = y*(1+8*x/(3*(3*y*y+x))) // = y*(y*y+3*x)/(3*y*y+x)
    // reconstruct original scale
    return y*scale;
}

可以用两个牛顿步代替哈雷步，或者在 c 中具有更好的均匀近似可以用牛顿步代替哈雷步，或者...