math: special cases for Modf, Frexp and Ldexp; added Modf_386

Also moved Modf from bits.go into modf.go and added timing tests.

R=rsc
CC=golang-dev
https://golang.org/cl/202042
diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index f228084..0e89df7 100644
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -18,6 +18,7 @@
 	fmod_386.$O\
 	hypot_386.$O\
 	log_386.$O\
+	modf_386.$O\
 	sin_386.$O\
 	sqrt_386.$O\
 	tan_386.$O\
@@ -44,6 +45,7 @@
 	hypot.go\
 	log.go\
 	log1p.go\
+	modf.go\
 	pow.go\
 	pow10.go\
 	sin.go\
diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index fa9c4bd..c0ac152 100644
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -1,4 +1,4 @@
-// Copyright 2009-2010 The Go Authors. All rights reserved.
+// Copyright 2009 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.
 
@@ -220,6 +220,25 @@
 	8.734595415957246977711748e-01,
 	1.314075231424398637614104e+00,
 }
+
+type fi struct {
+	f float64
+	i int
+}
+
+var frexp = []fi{
+	fi{6.2237649061045918750e-01, 3},
+	fi{9.6735905932226306250e-01, 3},
+	fi{-5.5376011438400318000e-01, -1},
+	fi{-6.2632545228388436250e-01, 3},
+	fi{6.02268356699901081250e-01, 4},
+	fi{7.3159430981099115000e-01, 2},
+	fi{6.5363542893241332500e-01, 3},
+	fi{6.8198497760900255000e-01, 2},
+	fi{9.1265404584042750000e-01, 1},
+	fi{-5.4287029803597508250e-01, 4},
+}
+
 var log = []float64{
 	1.605231462693062999102599e+00,
 	2.0462560018708770653153909e+00,
@@ -256,6 +275,18 @@
 	1.8088493239630770262045333e-02,
 	-9.0865245631588989681559268e-02,
 }
+var modf = [][2]float64{
+	[2]float64{4.0000000000000000e+00, 9.7901192488367350108546816e-01},
+	[2]float64{7.0000000000000000e+00, 7.3887247457810456552351752e-01},
+	[2]float64{0.0000000000000000e+00, -2.7688005719200159404635997e-01},
+	[2]float64{-5.0000000000000000e+00, -1.060361827107492160848778e-02},
+	[2]float64{9.0000000000000000e+00, 6.3629370719841737980004837e-01},
+	[2]float64{2.0000000000000000e+00, 9.2637723924396464525443662e-01},
+	[2]float64{5.0000000000000000e+00, 2.2908343145930665230025625e-01},
+	[2]float64{2.0000000000000000e+00, 7.2793991043601025126008608e-01},
+	[2]float64{1.0000000000000000e+00, 8.2530809168085506044576505e-01},
+	[2]float64{-8.0000000000000000e+00, -6.8592476857560136238589621e-01},
+}
 var pow = []float64{
 	9.5282232631648411840742957e+04,
 	5.4811599352999901232411871e+07,
@@ -516,6 +547,19 @@
 	NaN(),
 }
 
+var vffrexpSC = []float64{
+	Inf(-1),
+	0,
+	Inf(1),
+	NaN(),
+}
+var frexpSC = []fi{
+	fi{Inf(-1), 0},
+	fi{0, 0},
+	fi{Inf(1), 0},
+	fi{NaN(), 0},
+}
+
 var vfhypotSC = [][2]float64{
 	[2]float64{Inf(-1), Inf(-1)},
 	[2]float64{Inf(-1), 0},
@@ -581,6 +625,17 @@
 	NaN(),
 }
 
+var vfmodfSC = []float64{
+	Inf(-1),
+	Inf(1),
+	NaN(),
+}
+var modfSC = [][2]float64{
+	[2]float64{Inf(-1), NaN()},
+	[2]float64{Inf(1), NaN()},
+	[2]float64{NaN(), NaN()},
+}
+
 var vfpowSC = [][2]float64{
 	[2]float64{-Pi, Pi},
 	[2]float64{-Pi, -Pi},
@@ -919,6 +974,19 @@
 	}
 }
 
+func TestFrexp(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f, j := Frexp(vf[i]); !veryclose(frexp[i].f, f) || frexp[i].i != j {
+			t.Errorf("Frexp(%g) = %g, %d, want %g, %d\n", vf[i], f, j, frexp[i].f, frexp[i].i)
+		}
+	}
+	for i := 0; i < len(vffrexpSC); i++ {
+		if f, j := Frexp(vffrexpSC[i]); !alike(frexpSC[i].f, f) || frexpSC[i].i != j {
+			t.Errorf("Frexp(%g) = %g, %d, want %g, %d\n", vffrexpSC[i], f, j, frexpSC[i].f, frexpSC[i].i)
+		}
+	}
+}
+
 func TestHypot(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
 		a := Fabs(1e200 * tanh[i] * Sqrt(2))
@@ -933,6 +1001,19 @@
 	}
 }
 
+func TestLdexp(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f := Ldexp(frexp[i].f, frexp[i].i); !veryclose(vf[i], f) {
+			t.Errorf("Ldexp(%g, %d) = %g, want %g\n", frexp[i].f, frexp[i].i, f, vf[i])
+		}
+	}
+	for i := 0; i < len(vffrexpSC); i++ {
+		if f := Ldexp(frexpSC[i].f, frexpSC[i].i); !alike(vffrexpSC[i], f) {
+			t.Errorf("Ldexp(%g, %d) = %g, want %g\n", frexpSC[i].f, frexpSC[i].i, f, vffrexpSC[i])
+		}
+	}
+}
+
 func TestLog(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
 		a := Fabs(vf[i])
@@ -985,6 +1066,19 @@
 	}
 }
 
+func TestModf(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f, g := Modf(vf[i]); !veryclose(modf[i][0], f) || !veryclose(modf[i][1], g) {
+			t.Errorf("Modf(%g) = %g, %g, want %g, %g\n", vf[i], f, g, modf[i][0], modf[i][1])
+		}
+	}
+	for i := 0; i < len(vfmodfSC); i++ {
+		if f, g := Modf(vfmodfSC[i]); !alike(modfSC[i][0], f) || !alike(modfSC[i][1], g) {
+			t.Errorf("Modf(%g) = %g, %g, want %g, %g\n", vfmodfSC[i], f, g, modfSC[i][0], modfSC[i][1])
+		}
+	}
+}
+
 func TestPow(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
 		if f := Pow(10, vf[i]); !close(pow[i], f) {
@@ -1220,12 +1314,24 @@
 	}
 }
 
+func BenchmarkFrexp(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Frexp(8)
+	}
+}
+
 func BenchmarkHypot(b *testing.B) {
 	for i := 0; i < b.N; i++ {
 		Hypot(3, 4)
 	}
 }
 
+func BenchmarkLdexp(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Ldexp(.5, 2)
+	}
+}
+
 func BenchmarkLog(b *testing.B) {
 	for i := 0; i < b.N; i++ {
 		Log(.5)
@@ -1244,6 +1350,12 @@
 	}
 }
 
+func BenchmarkModf(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Modf(1.5)
+	}
+}
+
 func BenchmarkPowInt(b *testing.B) {
 	for i := 0; i < b.N; i++ {
 		Pow(2, 2)
diff --git a/src/pkg/math/bits.go b/src/pkg/math/bits.go
index 9f93a19..ccbcf06 100644
--- a/src/pkg/math/bits.go
+++ b/src/pkg/math/bits.go
@@ -53,7 +53,14 @@
 // It returns frac and exp satisfying f == frac × 2<sup>exp</sup>,
 // with the absolute value of frac in the interval [½, 1).
 func Frexp(f float64) (frac float64, exp int) {
-	if f == 0 {
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	// special cases
+	switch {
+	case f == 0:
+		return
+	case f < -MaxFloat64 || f > MaxFloat64 || f != f: // IsInf(f, 0) || IsNaN(f):
+		frac = f
 		return
 	}
 	x := Float64bits(f)
@@ -67,6 +74,12 @@
 // Ldexp is the inverse of Frexp.
 // It returns frac × 2<sup>exp</sup>.
 func Ldexp(frac float64, exp int) float64 {
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	// special cases
+	if frac != frac { // IsNaN(frac)
+		return NaN()
+	}
 	x := Float64bits(frac)
 	exp += int(x>>shift) & mask
 	if exp <= 0 {
@@ -82,27 +95,3 @@
 	x |= uint64(exp) << shift
 	return Float64frombits(x)
 }
-
-// Modf returns integer and fractional floating-point numbers
-// that sum to f.
-// Integer and frac have the same sign as f.
-func Modf(f float64) (int float64, frac float64) {
-	if f < 1 {
-		if f < 0 {
-			int, frac = Modf(-f)
-			return -int, -frac
-		}
-		return 0, f
-	}
-
-	x := Float64bits(f)
-	e := uint(x>>shift)&mask - bias
-
-	// Keep the top 11+e bits, the integer part; clear the rest.
-	if e < 64-11 {
-		x &^= 1<<(64-11-e) - 1
-	}
-	int = Float64frombits(x)
-	frac = f - int
-	return
-}
diff --git a/src/pkg/math/modf.go b/src/pkg/math/modf.go
new file mode 100644
index 0000000..ae0c7c8
--- /dev/null
+++ b/src/pkg/math/modf.go
@@ -0,0 +1,33 @@
+// Copyright 2009 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
+
+// Modf returns integer and fractional floating-point numbers
+// that sum to f.  Both values have the same sign as f.
+//
+// Special cases are:
+//	Modf(+Inf) = +Inf, NaN
+//	Modf(-Inf) = -Inf, NaN
+//	Modf(NaN) = NaN, NaN
+func Modf(f float64) (int float64, frac float64) {
+	if f < 1 {
+		if f < 0 {
+			int, frac = Modf(-f)
+			return -int, -frac
+		}
+		return 0, f
+	}
+
+	x := Float64bits(f)
+	e := uint(x>>shift)&mask - bias
+
+	// Keep the top 11+e bits, the integer part; clear the rest.
+	if e < 64-11 {
+		x &^= 1<<(64-11-e) - 1
+	}
+	int = Float64frombits(x)
+	frac = f - int
+	return
+}
diff --git a/src/pkg/math/modf_386.s b/src/pkg/math/modf_386.s
new file mode 100644
index 0000000..5ccab98
--- /dev/null
+++ b/src/pkg/math/modf_386.s
@@ -0,0 +1,19 @@
+// 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 Modf(x float64) (int float64, frac float64)
+TEXT ·Modf(SB),7,$0
+	FMOVD   x+0(FP), F0  // F0=x
+	FMOVD   F0, F1       // F0=x, F1=x
+	FSTCW   -2(SP)       // save old Control Word
+	MOVW    -2(SP), AX
+	ORW     $0x0c00, AX  // Rounding Control set to truncate
+	MOVW    AX, -4(SP)   // store new Control Word
+	FLDCW   -4(SP)       // load new Control Word
+	FRNDINT              // F0=trunc(x), F1=x
+	FLDCW   -2(SP)       // load old Control Word
+	FSUBD   F0, F1       // F0=trunc(x), F1=x-trunc(x)
+	FMOVDP  F0, i+8(FP)  // F0=x-trunc(x)
+	FMOVDP  F0, f+16(FP)
+	RET
diff --git a/src/pkg/math/modf_decl.go b/src/pkg/math/modf_decl.go
new file mode 100644
index 0000000..7add2af
--- /dev/null
+++ b/src/pkg/math/modf_decl.go
@@ -0,0 +1,7 @@
+// 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
+
+func Modf(f float64) (int float64, frac float64)