]> Cypherpunks repositories - gostls13.git/commitdiff
math: added ilogb, logb, remainder, tests and special conditions
authorCharles L. Dorian <cldorian@gmail.com>
Thu, 4 Mar 2010 02:17:13 +0000 (18:17 -0800)
committerRuss Cox <rsc@golang.org>
Thu, 4 Mar 2010 02:17:13 +0000 (18:17 -0800)
Also added expm1_386 and remainder_386; shortened exp_386

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

src/pkg/math/Makefile
src/pkg/math/all_test.go
src/pkg/math/exp_386.s
src/pkg/math/expm1_386.s [new file with mode: 0644]
src/pkg/math/expm1_decl.go [new file with mode: 0644]
src/pkg/math/logb.go [new file with mode: 0644]
src/pkg/math/remainder.go [new file with mode: 0644]
src/pkg/math/remainder_386.s [new file with mode: 0644]
src/pkg/math/remainder_decl.go [new file with mode: 0644]

index e24c448f88d3d7f06606e118e3b52a643911a4ea..6650482a7e708567202bdc91f8dce4e7ca5daca7 100644 (file)
@@ -15,6 +15,7 @@ OFILES_386=\
        atan2_386.$O\
        exp_386.$O\
        exp2_386.$O\
+       expm1_386.$O\
        fabs_386.$O\
        floor_386.$O\
        frexp_386.$O\
@@ -24,6 +25,7 @@ OFILES_386=\
        log_386.$O\
        log1p_386.$O\
        modf_386.$O\
+       remainder_386.$O\
        sin_386.$O\
        sincos_386.$O\
        sqrt_386.$O\
@@ -52,6 +54,7 @@ ALLGOFILES=\
        fmod.go\
        frexp.go\
        hypot.go\
+       logb.go\
        lgamma.go\
        ldexp.go\
        log.go\
@@ -60,6 +63,7 @@ ALLGOFILES=\
        nextafter.go\
        pow.go\
        pow10.go\
+       remainder.go\
        sin.go\
        sincos.go\
        sinh.go\
index d80f4ee13378cbddac78e449e364d2cefdd09351..6279499713a1b79ed79b08831268c92db201fa9f 100644 (file)
@@ -310,6 +310,18 @@ var log = []float64{
        6.0174879014578057187016475e-01,
        2.161703872847352815363655e+00,
 }
+var logb = []float64{
+       3.0000000000000000e+00,
+       3.0000000000000000e+00,
+       -1.0000000000000000e+00,
+       3.0000000000000000e+00,
+       4.0000000000000000e+00,
+       2.0000000000000000e+00,
+       3.0000000000000000e+00,
+       2.0000000000000000e+00,
+       1.0000000000000000e+00,
+       4.0000000000000000e+00,
+}
 var log10 = []float64{
        6.9714316642508290997617083e-01,
        8.886776901739320576279124e-01,
@@ -382,6 +394,18 @@ var pow = []float64{
        6.688182138451414936380374e+01,
        2.0609869004248742886827439e-09,
 }
+var remainder = []float64{
+       4.197615023265299782906368e-02,
+       2.261127525421895434476482e+00,
+       3.231794108794261433104108e-02,
+       -2.120723654214984321697556e-02,
+       3.637062928015826201999516e-01,
+       1.220868282268106064236690e+00,
+       -4.581668629186133046005125e-01,
+       -9.117596417440410050403443e-01,
+       8.734595415957246977711748e-01,
+       1.314075231424398637614104e+00,
+}
 var sin = []float64{
        -9.6466616586009283766724726e-01,
        9.9338225271646545763467022e-01,
@@ -747,6 +771,19 @@ var hypotSC = []float64{
        NaN(),
 }
 
+var vfilogbSC = []float64{
+       Inf(-1),
+       0,
+       Inf(1),
+       NaN(),
+}
+var ilogbSC = []int{
+       MaxInt32,
+       MinInt32,
+       MaxInt32,
+       MaxInt32,
+}
+
 var vflgammaSC = []float64{
        Inf(-1),
        -3,
@@ -777,6 +814,19 @@ var logSC = []float64{
        NaN(),
 }
 
+var vflogbSC = []float64{
+       Inf(-1),
+       0,
+       Inf(1),
+       NaN(),
+}
+var logbSC = []float64{
+       Inf(1),
+       Inf(-1),
+       Inf(1),
+       NaN(),
+}
+
 var vflog1pSC = []float64{
        Inf(-1),
        -Pi,
@@ -1204,7 +1254,7 @@ func TestFmin(t *testing.T) {
 
 func TestFmod(t *testing.T) {
        for i := 0; i < len(vf); i++ {
-               if f := Fmod(10, vf[i]); fmod[i] != f { /*!close(fmod[i], f)*/
+               if f := Fmod(10, vf[i]); fmod[i] != f {
                        t.Errorf("Fmod(10, %g) = %g, want %g\n", vf[i], f, fmod[i])
                }
        }
@@ -1242,6 +1292,19 @@ func TestHypot(t *testing.T) {
        }
 }
 
+func TestIlogb(t *testing.T) {
+       for i := 0; i < len(vf); i++ {
+               if e := Ilogb(vf[i]); frexp[i].i != e {
+                       t.Errorf("Ilogb(%g) = %d, want %d\n", vf[i], e, frexp[i].i)
+               }
+       }
+       for i := 0; i < len(vflogbSC); i++ {
+               if e := Ilogb(vflogbSC[i]); ilogbSC[i] != e {
+                       t.Errorf("Ilogb(%g) = %d, want %d\n", vflogbSC[i], e, ilogbSC[i])
+               }
+       }
+}
+
 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) {
@@ -1285,6 +1348,19 @@ func TestLog(t *testing.T) {
        }
 }
 
+func TestLogb(t *testing.T) {
+       for i := 0; i < len(vf); i++ {
+               if f := Logb(vf[i]); logb[i] != f {
+                       t.Errorf("Logb(%g) = %g, want %g\n", vf[i], f, logb[i])
+               }
+       }
+       for i := 0; i < len(vflogbSC); i++ {
+               if f := Logb(vflogbSC[i]); !alike(logbSC[i], f) {
+                       t.Errorf("Logb(%g) = %g, want %g\n", vflogbSC[i], f, logbSC[i])
+               }
+       }
+}
+
 func TestLog10(t *testing.T) {
        for i := 0; i < len(vf); i++ {
                a := Fabs(vf[i])
@@ -1376,6 +1452,19 @@ func TestPow(t *testing.T) {
        }
 }
 
+func TestRemainder(t *testing.T) {
+       for i := 0; i < len(vf); i++ {
+               if f := Remainder(10, vf[i]); remainder[i] != f {
+                       t.Errorf("Remainder(10, %g) = %g, want %g\n", vf[i], f, remainder[i])
+               }
+       }
+       for i := 0; i < len(vffmodSC); i++ {
+               if f := Remainder(vffmodSC[i][0], vffmodSC[i][1]); !alike(fmodSC[i], f) {
+                       t.Errorf("Remainder(%g, %g) = %g, want %g\n", vffmodSC[i][0], vffmodSC[i][1], f, fmodSC[i])
+               }
+       }
+}
+
 func TestSin(t *testing.T) {
        for i := 0; i < len(vf); i++ {
                if f := Sin(vf[i]); !close(sin[i], f) {
@@ -1665,6 +1754,12 @@ func BenchmarkHypot(b *testing.B) {
        }
 }
 
+func BenchmarkIlogb(b *testing.B) {
+       for i := 0; i < b.N; i++ {
+               Ilogb(.5)
+       }
+}
+
 func BenchmarkLdexp(b *testing.B) {
        for i := 0; i < b.N; i++ {
                Ldexp(.5, 2)
@@ -1683,6 +1778,12 @@ func BenchmarkLog(b *testing.B) {
        }
 }
 
+func BenchmarkLogb(b *testing.B) {
+       for i := 0; i < b.N; i++ {
+               Logb(.5)
+       }
+}
+
 func BenchmarkLog10(b *testing.B) {
        for i := 0; i < b.N; i++ {
                Log10(.5)
@@ -1725,6 +1826,12 @@ func BenchmarkPowFrac(b *testing.B) {
        }
 }
 
+func BenchmarkRemainder(b *testing.B) {
+       for i := 0; i < b.N; i++ {
+               Remainder(10, 3)
+       }
+}
+
 func BenchmarkSin(b *testing.B) {
        for i := 0; i < b.N; i++ {
                Sin(.5)
index 5121f3ec1de4d55f956d854a53cfbe66c9179fd8..e0743e72a22d2e77b0116d316c6bd8073f4ca465 100644 (file)
@@ -10,8 +10,7 @@ TEXT ·Exp(SB),7,$0
        CMPL    AX, $0x7ff00000
        JEQ     not_finite
        FLDL2E                // F0=log2(e)
-       FMOVD   x+0(FP), F0   // F0=x, F1=log2(e)
-       FMULDP  F0, F1        // F0=x*log2(e)
+       FMULD   x+0(FP), F0   // F0=x*log2(e)
        FMOVD   F0, F1        // F0=x*log2(e), F1=x*log2(e)
        FRNDINT               // F0=int(x*log2(e)), F1=x*log2(e)
        FSUBD   F0, F1        // F0=int(x*log2(e)), F1=x*log2(e)-int(x*log2(e))
@@ -31,8 +30,8 @@ not_finite:
        JNE     not_neginf
        CMPL    CX, $0
        JNE     not_neginf
-       MOVL    $0, r+8(FP)
-       MOVL    $0, r+12(FP)
+       FLDZ                  // F0=0
+       FMOVDP  F0, r+8(FP)
        RET
 not_neginf:
        MOVL    CX, r+8(FP)
diff --git a/src/pkg/math/expm1_386.s b/src/pkg/math/expm1_386.s
new file mode 100644 (file)
index 0000000..8185f49
--- /dev/null
@@ -0,0 +1,55 @@
+// 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 Expm1(x float64) float64
+TEXT ·Expm1(SB),7,$0
+       FLDLN2               // F0=log(2) = 1/log2(e) ~ 0.693147
+       FMOVD   x+0(FP), F0  // F0=x, F1=1/log2(e)
+       FABS                 // F0=|x|, F1=1/log2(e) 
+       FUCOMPP F0, F1       // compare F0 to F1
+       FSTSW   AX
+       SAHF
+       JCC     use_exp      // jump if F0 >= F1
+       FLDL2E                // F0=log2(e)
+       FMULD   x+0(FP), F0   // F0=x*log2(e) (-1<F0<1)
+       F2XM1                 // F0=e**x-1 = 2**(x*log2(e))-1
+       FMOVDP  F0, r+8(FP)
+       RET
+use_exp:
+// test bits for not-finite
+       MOVL    x+4(FP), AX
+       ANDL    $0x7ff00000, AX
+       CMPL    AX, $0x7ff00000
+       JEQ     not_finite
+       FLDL2E                // F0=log2(e)
+       FMULD   x+0(FP), F0   // F0=x*log2(e)
+       FMOVD   F0, F1        // F0=x*log2(e), F1=x*log2(e)
+       FRNDINT               // F0=int(x*log2(e)), F1=x*log2(e)
+       FSUBD   F0, F1        // F0=int(x*log2(e)), F1=x*log2(e)-int(x*log2(e))
+       FXCHD   F0, F1        // F0=x*log2(e)-int(x*log2(e)), F1=int(x*log2(e))
+       F2XM1                 // F0=2**(x*log2(e)-int(x*log2(e)))-1, F1=int(x*log2(e))
+       FLD1                  // F0=1, F1=2**(x*log2(e)-int(x*log2(e)))-1, F2=int(x*log2(e))
+       FADDDP  F0, F1        // F0=2**(x*log2(e)-int(x*log2(e))), F1=int(x*log2(e))
+       FSCALE                // F0=e**x, F1=int(x*log2(e))
+       FMOVDP  F0, F1        // F0=e**x
+       FLD1                  // F0=1, F1=e**x
+       FSUBDP  F0, F1        // F0=e**x-1 
+       FMOVDP  F0, r+8(FP)
+       RET
+not_finite:
+// test bits for -Inf
+       MOVL    x+4(FP), BX
+       MOVL    x+0(FP), CX
+       CMPL    BX, $0xfff00000
+       JNE     not_neginf
+       CMPL    CX, $0
+       JNE     not_neginf
+       FLD1                 // F0=1
+       FCHS                 // F0=-1
+       FMOVDP  F0, r+8(FP)
+       RET
+not_neginf:
+       MOVL    CX, r+8(FP)
+       MOVL    BX, r+12(FP)
+       RET
diff --git a/src/pkg/math/expm1_decl.go b/src/pkg/math/expm1_decl.go
new file mode 100644 (file)
index 0000000..4dab70b
--- /dev/null
@@ -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 Expm1(x float64) float64
diff --git a/src/pkg/math/logb.go b/src/pkg/math/logb.go
new file mode 100644 (file)
index 0000000..acda15d
--- /dev/null
@@ -0,0 +1,47 @@
+// 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
+
+// Logb(x) returns the binary logarithm of non-zero x.
+//
+// Special cases are:
+//     Logb(±Inf) = +Inf
+//     Logb(0) = -Inf
+//     Logb(NaN) = NaN
+func Logb(x float64) float64 {
+       // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+       // when compiler does it for us
+       // special cases
+       switch {
+       case x == 0:
+               return Inf(-1)
+       case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
+               return Inf(1)
+       case x != x: // IsNaN(x):
+               return x
+       }
+       return float64(int((Float64bits(x)>>shift)&mask) - bias)
+}
+
+// Ilogb(x) returns the binary logarithm of non-zero x as an integer.
+//
+// Special cases are:
+//     Ilogb(±Inf) = MaxInt32
+//     Ilogb(0) = MinInt32
+//     Ilogb(NaN) = MaxInt32
+func Ilogb(x float64) int {
+       // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+       // when compiler does it for us
+       // special cases
+       switch {
+       case x == 0:
+               return MinInt32
+       case x != x: // IsNaN(x):
+               return MaxInt32
+       case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
+               return MaxInt32
+       }
+       return int((Float64bits(x)>>shift)&mask) - bias
+}
diff --git a/src/pkg/math/remainder.go b/src/pkg/math/remainder.go
new file mode 100644 (file)
index 0000000..be8724c
--- /dev/null
@@ -0,0 +1,85 @@
+// 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
+
+// The original C code and the the comment below are from
+// FreeBSD's /usr/src/lib/msun/src/e_remainder.c and came
+// with this notice.  The go code is a simplified version of
+// the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// __ieee754_remainder(x,y)
+// Return :
+//      returns  x REM y  =  x - [x/y]*y  as if in infinite
+//      precision arithmetic, where [x/y] is the (infinite bit)
+//      integer nearest x/y (in half way cases, choose the even one).
+// Method :
+//      Based on fmod() returning  x - [x/y]chopped * y  exactly.
+
+// Remainder returns the IEEE 754 floating-point remainder of x/y.
+//
+// Special cases are:
+//     Remainder(x, NaN) = NaN
+//     Remainder(NaN, y) = NaN
+//     Remainder(Inf, y) = NaN
+//     Remainder(x, 0) = NaN
+//     Remainder(x, Inf) = x
+func Remainder(x, y float64) float64 {
+       const (
+               Tiny    = 4.45014771701440276618e-308 // 0x0020000000000000
+               HalfMax = MaxFloat64 / 2
+       )
+       // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+       // when compiler does it for us
+       // special cases
+       switch {
+       case x != x || y != y || x < -MaxFloat64 || x > MaxFloat64 || y == 0: // IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0:
+               return NaN()
+       case y < -MaxFloat64 || y > MaxFloat64: // IsInf(y):
+               return x
+       }
+       sign := false
+       if x < 0 {
+               x = -x
+               sign = true
+       }
+       if y < 0 {
+               y = -y
+       }
+       if x == y {
+               return 0
+       }
+       if y <= HalfMax {
+               x = Fmod(x, y+y) // now x < 2y
+       }
+       if y < Tiny {
+               if x+x > y {
+                       x -= y
+                       if x+x >= y {
+                               x -= y
+                       }
+               }
+       } else {
+               yHalf := 0.5 * y
+               if x > yHalf {
+                       x -= y
+                       if x >= yHalf {
+                               x -= y
+                       }
+               }
+       }
+       if sign {
+               x = -x
+       }
+       return x
+}
diff --git a/src/pkg/math/remainder_386.s b/src/pkg/math/remainder_386.s
new file mode 100644 (file)
index 0000000..4cb9823
--- /dev/null
@@ -0,0 +1,15 @@
+// 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 Remainder(x, y float64) float64
+TEXT ·Remainder(SB),7,$0
+       FMOVD   y+8(FP), F0  // F0=y
+       FMOVD   x+0(FP), F0  // F0=x, F1=y
+       FPREM1               // F0=reduced_x, F1=y
+       FSTSW   AX           // AX=status word
+       ANDW    $0x0400, AX
+       JNE     -3(PC)       // jump if reduction incomplete
+       FMOVDP  F0, F1       // F0=x-q*y
+       FMOVDP  F0, r+16(FP)
+       RET
diff --git a/src/pkg/math/remainder_decl.go b/src/pkg/math/remainder_decl.go
new file mode 100644 (file)
index 0000000..1407d9a
--- /dev/null
@@ -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 Remainder(x, y float64) float64