指數計算速度

Question

我目前正在測試Julia（我曾與Matlab合作）

在matlab中，N^3的計算速度比NxNxN慢。 N^2和NxN不會發生這種情況。他們使用不同的算法來計算高階指數，因爲他們更喜歡精確性而非速度。

我認爲朱莉婭做同樣的事情。

我想問一下，是否有辦法強制Julia使用乘法而不是默認算法來計算N的指數，至少對於立方體指數來說。

前段時間我做了一些關於matlab的測試。我將該代碼翻譯成了julia。

鏈接代碼： http://pastebin.com/bbeukhTc （我不能在這裏上傳的所有鏈接:(）Matlab的2014年劇本

結果：

Exponente1

經過時間是68.293793秒（最小的17.7倍）

指數2

已用時間爲24.236218秒。 （該smallests的6.3倍倍）

Exponente3

經過時間3.853348秒。

上朱莉婭0.46腳本的結果：

Exponente1

18.423204秒（8.22ķ分配：372.563 KB）（最小51.6x倍）

Exponente2

13.746904秒（9.02 k分配：407.332 KB）（最小的38.5倍）

Exponente3

0.356875秒（10.01ķ分配：450.441 KB）

在我的測試朱莉婭比Matlab的速度更快，但我使用的是相對舊版本。我無法測試其他版本。

Answer 1

檢查Julia的源代碼：

julia/base/math.jl：

^(x::Float64, y::Integer) = 
    box(Float64, powi_llvm(unbox(Float64,x), unbox(Int32,Int32(y)))) 
^(x::Float32, y::Integer) = 
    box(Float32, powi_llvm(unbox(Float32,x), unbox(Int32,Int32(y)))) 

julia/base/fastmath.jl：

pow_fast{T<:FloatTypes}(x::T, y::Integer) = pow_fast(x, Int32(y)) 
pow_fast{T<:FloatTypes}(x::T, y::Int32) = 
    box(T, Base.powi_llvm(unbox(T,x), unbox(Int32,y))) 

我們可以看到，朱莉婭使用powi_llvm

檢查llvm's source code：

define double @powi(double %F, i32 %power) { 
; CHECK: powi: 
; CHECK: bl __powidf2 
     %result = call double @llvm.powi.f64(double %F, i32 %power) 
    ret double %result 
} 

現在，__powidf2是有趣的函數這裏：

COMPILER_RT_ABI double 
__powidf2(double a, si_int b) 
{ 
    const int recip = b < 0; 
    double r = 1; 
    while (1) 
    { 
     if (b & 1) 
      r *= a; 
     b /= 2; 
     if (b == 0) 
      break; 
     a *= a; 
    } 
    return recip ? 1/r : r; 
} 

實施例1：給定一個= 2; b = 7：

-    r   = 1 
- iteration 1: r = 1 * 2 = 2; b = (int)(7/2) = 3; a = 2 * 2 = 4 
- iteration 2: r = 2 * 4 = 8; b = (int)(3/2) = 1; a = 4 * 4 = 16 
- iteration 3: r = 8 * 16 = 128; 

例2：給定a = 2; b = 8：

-    r   = 1 
- iteration 1: r   = 1; b = (int)(8/2) = 4; a = 2 * 2 = 4 
- iteration 2: r   = 1; b = (int)(4/2) = 2; a = 4 * 4 = 16 
- iteration 3: r   = 1; b = (int)(2/2) = 1; a = 16 * 16 = 256 
- iteration 4: r = 1 * 256 = 256; b = (int)(1/2) = 0; 

整數功率總是作爲序列乘法實現的。這就是爲什麼N^3比N^2慢。

jl_powi_llvm（在fastmath.jl中調用「jl_」由宏擴展連接）另一方面將指數轉換爲浮點並調用pow（）。 C源代碼：

JL_DLLEXPORT jl_value_t *jl_powi_llvm(jl_value_t *a, jl_value_t *b) 
{ 
    jl_value_t *ty = jl_typeof(a); 
    if (!jl_is_bitstype(ty)) 
     jl_error("powi_llvm: a is not a bitstype"); 
    if (!jl_is_bitstype(jl_typeof(b)) || jl_datatype_size(jl_typeof(b)) != 4) 
     jl_error("powi_llvm: b is not a 32-bit bitstype"); 
    jl_value_t *newv = newstruct((jl_datatype_t*)ty); 
    void *pa = jl_data_ptr(a), *pr = jl_data_ptr(newv); 
    int sz = jl_datatype_size(ty); 
    switch (sz) { 
    /* choose the right size c-type operation */ 
    case 4: 
     *(float*)pr = powf(*(float*)pa, (float)jl_unbox_int32(b)); 
     break; 
    case 8: 
     *(double*)pr = pow(*(double*)pa, (double)jl_unbox_int32(b)); 
     break; 
    default: 
     jl_error("powi_llvm: runtime floating point intrinsics are not implemented for bit sizes other than 32 and 64"); 
    } 
    return newv; 
} 

Answer 2

Lior的回答非常好。以下是您提出的問題的解決方案：是的，有一種強制使用乘法的方法，但要以準確性爲代價。這是@fastmath宏：

julia> @benchmark 1.1^3 
BenchmarkTools.Trial: 
    samples:   10000 
    evals/sample:  999 
    time tolerance: 5.00% 
    memory tolerance: 1.00% 
    memory estimate: 16.00 bytes 
    allocs estimate: 1 
    minimum time:  13.00 ns (0.00% GC) 
    median time:  14.00 ns (0.00% GC) 
    mean time:  15.74 ns (6.14% GC) 
    maximum time:  1.85 μs (98.16% GC) 

julia> @benchmark @fastmath 1.1^3 
BenchmarkTools.Trial: 
    samples:   10000 
    evals/sample:  1000 
    time tolerance: 5.00% 
    memory tolerance: 1.00% 
    memory estimate: 0.00 bytes 
    allocs estimate: 0 
    minimum time:  2.00 ns (0.00% GC) 
    median time:  3.00 ns (0.00% GC) 
    mean time:  2.59 ns (0.00% GC) 
    maximum time:  20.00 ns (0.00% GC) 

注意與@fastmath，性能要好得多。