超快的Java乘方算法在Ruby里面没用？？

有丘直方

你们知道我在搞RGSS的物理引擎
最近在研究JBox2D的源码
看到里面有这一段[pre lang="java"]package org.jbox2d.common;

/**
* Contains methods from MathUtils that rely on JVM features. These are separated out from
* MathUtils so that they can be overridden when compiling for GWT.
*/
class PlatformMathUtils {

  private static final float SHIFT23 = 1 << 23;
  private static final float INV_SHIFT23 = 1.0f / SHIFT23;

  public static final float fastPow(float a, float b) {
    float x = Float.floatToRawIntBits(a);
    x *= INV_SHIFT23;
    x -= 127;
    float y = x - (x >= 0 ? (int) x : (int) x - 1);
    b *= x + (y - y * y) * 0.346607f;
    y = b - (b >= 0 ? (int) b : (int) b - 1);
    y = (y - y * y) * 0.33971f;
    return Float.intBitsToFloat((int) ((b + 127 - y) * SHIFT23));
  }
}
[/pre]这是个用于快速计算乘方的算法，精度不高，根号2算出来大概是1.415
为了测试这个算法到底有多快，我写了一段代码进行测试[pre lang="java"]public class Test {

    private static final float SHIFT23 = 1 << 23;
    private static final float INV_SHIFT23 = 1.0f / SHIFT23;

    public static final float fastPow(float a, float b) {
      float x = Float.floatToRawIntBits(a);
      x *= INV_SHIFT23;
      x -= 127;
      float y = x - (x >= 0 ? (int) x : (int) x - 1);
      b *= x + (y - y * y) * 0.346607f;
      y = b - (b >= 0 ? (int) b : (int) b - 1);
      y = (y - y * y) * 0.33971f;
      return Float.intBitsToFloat((int) ((b + 127 - y) * SHIFT23));
    }
   
    public static long times(int timesCount, Runnable runnable) {
        long time = System.currentTimeMillis();
        for (int i = 0; i < timesCount; i++) {
            runnable.run();
        }
        return System.currentTimeMillis() - time;
    }
   
    public static void main(String[] args) {
        System.out.println(Test.times(100000000, () -> Test.fastPow(1.5f, 3.0f)));
        System.out.println(Test.times(100000000, () -> StrictMath.pow(1.5f, 3.0f)));
    }
}[/pre]原理是，执行一亿次新算法，输出消耗的时间，再执行一亿次StrictMath自带的算法，输出消耗时间
结果吓我一大跳：[pre lang=""]7
16918[/pre]新算法可以在7毫秒内执行一亿次乘方运算，比标准库自带算法快了2400多倍？
这么好的算法一定要移植到Ruby中用啊！[pre lang="ruby"]SHIFT23 = (1 << 23).to_f
INV_SHIFT23 = 1.0 / SHIFT23
def fast_pow(a, b)
  x = a.int_bits
  x *= INV_SHIFT23
  x -= 0x7f
  y = x - (x >= 0 ? x.to_i : x.to_i - 1)
  b *= x + (y - y * y) * 0.346607
  y = b - (b >= 0 ? b.to_i : b.to_i - 1)
  y = (y - y * y) * 0.33971
  Float.get_from_int_bits(((b + 0x7f - y) * SHIFT23).to_i)
end
class Float
  def int_bits
    [self].pack('g').unpack('B*').first.to_i(2)
  end
  class << self
    def get_from_int_bits(i)
      [sprintf("%032b", i)].pack('B*').unpack('g').first
    end
  end
end
def times(n, &block)
  time = Time.now
  n.times(&block)
  Time.now - time
end

puts(times(100000000) { fast_pow(1.5, 3.0) })
puts(times(100000000) { 1.5 ** 3.0 })[/pre]代码内还附带了根据IEEE 754浮点“单一格式”位布局返回指定浮点值的表示形式的方法，和返回对应于给定位表示形式的浮点数值的方法
运行结果再次让我大吃一惊：[pre lang=""]223.896929
8.523068[/pre]

guoxiaomi

虽然不知道为什么，但是你这个程序算出来的

1. fast_pow(1.5, 3.0) = 3.395968437194824
   
2. 1.5 ** 3.0 = 3.375

复制代码

测试了一下，其中fast_pow2是把第7和9行改成了 y = x % 1和 y = b % 1[pre lang="ruby"]Benchmark.bm(10) do |t|
  t.report('**_pow'){
    1_000_000.times{1.5 ** 3.0}
  }
  t.report('fast_pow'){
    1_000_000.times{fast_pow(1.5, 3.0)}
  }
  t.report('fast_pow2'){
    1_000_000.times{fast_pow2(1.5, 3.0)}
  }
  a = 1.5
  t.report('init_bits'){
    1_000_000.times{x = a.int_bits}
  }
  b = 3.0
  y = 0.06197100837140948
  t.report('get_from'){
    1_000_000.times{Float.get_from_int_bits(((b + 0x7f - y) * SHIFT23).to_i)}
  }
end[/pre]结果如下：[pre lang=""]                 user     system      total        real
**_pow       0.093000   0.000000   0.093000 (  0.092967)
fast_pow     2.485000   0.000000   2.485000 (  2.485887)
fast_pow2    2.406000   0.000000   2.406000 (  2.417023)
init_bits    0.625000   0.000000   0.625000 (  0.621773)
get_from     1.453000   0.000000   1.453000 (  1.463776)[/pre]