math: amd64 versions of exp and fabs

Benchmark of exp to 28 ns/op from 64 ns/op,
on 2.53GHz Intel Core 2 Duo.

R=rsc
CC=golang-dev
https://golang.org/cl/1594041

diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index 1447fc11d405d0d1e41638e8c31757628a1d5554..a2d11e43d732a6ccbf5e7d2fd47b54dd9d419757 100644 (file)
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -7,6 +7,8 @@ include ../../Make.$(GOARCH)
 TARG=math
 
 OFILES_amd64=\
+       exp_amd64.$O\
+       fabs_amd64.$O\
        fdim_amd64.$O\
        sqrt_amd64.$O\
 
@@ -48,6 +50,7 @@ ALLGOFILES=\
        copysign.go\
        erf.go\
        exp.go\
+       exp2.go\
        expm1.go\
        fabs.go\
        fdim.go\

diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 9a4801b006ca2a709e375935613aac3452d839d9..18a3f1b31300f2a489e0bc91b28e2316c985fcd9 100644 (file)
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -2343,9 +2343,9 @@ func BenchmarkExp2(b *testing.B) {
        }
 }
 
-func BenchmarkFloor(b *testing.B) {
+func BenchmarkFabs(b *testing.B) {
        for i := 0; i < b.N; i++ {
-               Floor(.5)
+               Fabs(.5)
        }
 }
 
@@ -2355,6 +2355,12 @@ func BenchmarkFdim(b *testing.B) {
        }
 }
 
+func BenchmarkFloor(b *testing.B) {
+       for i := 0; i < b.N; i++ {
+               Floor(.5)
+       }
+}
+
 func BenchmarkFmax(b *testing.B) {
        for i := 0; i < b.N; i++ {
                Fmax(10, 3)
@@ -2445,15 +2451,15 @@ func BenchmarkLogb(b *testing.B) {
        }
 }
 
-func BenchmarkLog10(b *testing.B) {
+func BenchmarkLog1p(b *testing.B) {
        for i := 0; i < b.N; i++ {
-               Log10(.5)
+               Log1p(.5)
        }
 }
 
-func BenchmarkLog1p(b *testing.B) {
+func BenchmarkLog10(b *testing.B) {
        for i := 0; i < b.N; i++ {
-               Log1p(.5)
+               Log10(.5)
        }
 }
 

diff --git a/src/pkg/math/exp.go b/src/pkg/math/exp.go
index b801d58acd20ad406fb73844e9d83ebc96741366..90409c341b6f39c084551113ad466467e6f78b1e 100644 (file)
--- a/src/pkg/math/exp.go
+++ b/src/pkg/math/exp.go
@@ -139,8 +139,3 @@ func Exp(x float64) float64 {
        // TODO(rsc): make sure Ldexp can handle boundary k
        return Ldexp(y, k)
 }
-
-// Exp2 returns 2**x, the base-2 exponential of x.
-//
-// Special cases are the same as Exp.
-func Exp2(x float64) float64 { return Exp(x * Ln2) }

diff --git a/src/pkg/math/exp2.go b/src/pkg/math/exp2.go
new file mode 100644 (file)
index 0000000..1e67f29
--- /dev/null
+++ b/src/pkg/math/exp2.go
@@ -0,0 +1,10 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+// Exp2 returns 2**x, the base-2 exponential of x.
+//
+// Special cases are the same as Exp.
+func Exp2(x float64) float64 { return Exp(x * Ln2) }

diff --git a/src/pkg/math/exp_amd64.s b/src/pkg/math/exp_amd64.s
new file mode 100644 (file)
index 0000000..844b5c9
--- /dev/null
+++ b/src/pkg/math/exp_amd64.s
@@ -0,0 +1,104 @@
+// Copyright 2010 The Go Authors.  All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// The method is based on a paper by Naoki Shibata: "Efficient evaluation
+// methods of elementary functions suitable for SIMD computation", Proc.
+// of International Supercomputing Conference 2010 (ISC'10), pp. 25 -- 32
+// (May 2010). The paper is available at
+// http://www.springerlink.com/content/340228x165742104/
+//
+// The original code and the constants below are from the author's
+// implementation available at http://freshmeat.net/projects/sleef.
+// The README file says, "The software is in public domain.
+// You can use the software without any obligation."
+//
+// This code is a simplified version of the original.
+
+#define LN2 0.6931471805599453094172321214581766 // log_e(2)
+#define LOG2E 1.4426950408889634073599246810018920 // 1/LN2
+#define LN2U 0.69314718055966295651160180568695068359375 // upper half LN2
+#define LN2L 0.28235290563031577122588448175013436025525412068e-12 // lower half LN2
+
+// func Exp(x float64) float64
+TEXT ·Exp(SB),7,$0
+// test bits for not-finite
+       MOVQ    x+0(FP), AX
+       MOVQ    $0x7ff0000000000000, BX
+       ANDQ    BX, AX
+       CMPQ    BX, AX
+       JEQ     not_finite
+       MOVSD   x+0(FP), X0
+       MOVSD   $LOG2E, X1
+       MULSD   X0, X1
+       CVTTSD2SQ X1, BX // BX = exponent
+       CVTSQ2SD BX, X1
+       MOVSD   $LN2U, X2
+       MULSD   X1, X2
+       SUBSD   X2, X0
+       MOVSD   $LN2L, X2
+       MULSD   X1, X2
+       SUBSD   X2, X0
+       // reduce argument
+       MOVSD   $0.0625, X1
+       MULSD   X1, X0
+       // Taylor series evaluation
+       MOVSD   $2.4801587301587301587e-5, X1
+       MULSD   X0, X1
+       MOVSD   $1.9841269841269841270e-4, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $1.3888888888888888889e-3, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $8.3333333333333333333e-3, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $4.1666666666666666667e-2, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $1.6666666666666666667e-1, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $0.5, X2
+       ADDSD   X2, X1
+       MULSD   X0, X1
+       MOVSD   $1.0, X2
+       ADDSD   X2, X1
+       MULSD   X1, X0
+       MOVSD   $2.0, X1
+       ADDSD   X0, X1
+       MULSD   X1, X0
+       MOVSD   $2.0, X1
+       ADDSD   X0, X1
+       MULSD   X1, X0
+       MOVSD   $2.0, X1
+       ADDSD   X0, X1
+       MULSD   X1, X0
+       MOVSD   $2.0, X1
+       ADDSD   X0, X1
+       MULSD   X1, X0
+       MOVSD   $1.0, X1
+       ADDSD   X1, X0
+       // return ldexp(fr, exp)
+       MOVQ    $0x3ff, AX // bias + 1
+       ADDQ    AX, BX
+       MOVQ    BX, X1
+       MOVQ    $52, AX // shift
+       MOVQ    AX, X2
+       PSLLQ   X2, X1
+       MULSD   X1, X0
+       MOVSD   X0, r+8(FP)
+       RET
+not_finite:
+// test bits for -Inf
+       MOVQ    x+0(FP), AX
+       MOVQ    $0xfff0000000000000, BX
+       CMPQ    BX, AX
+       JNE     not_neginf
+       XORQ    AX, AX
+       MOVQ    AX, r+8(FP)
+       RET
+not_neginf:
+       MOVQ    AX, r+8(FP)
+       RET

diff --git a/src/pkg/math/fabs_amd64.s b/src/pkg/math/fabs_amd64.s
new file mode 100644 (file)
index 0000000..8a9aedb
--- /dev/null
+++ b/src/pkg/math/fabs_amd64.s
@@ -0,0 +1,12 @@
+// Copyright 2010 The Go Authors.  All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// func Fabs(x float64) float64
+TEXT ·Fabs(SB),7,$0
+       MOVQ   $(1<<63), BX
+       MOVQ   BX, X0 // movsd $(-0.0), x0
+       MOVSD  x+0(FP), X1
+       ANDNPD X1, X0
+       MOVSD  X0, r+8(FP)
+       RET