diff options
| author | Ian Lance Taylor <iant@golang.org> | 2018-01-09 01:23:08 +0000 |
|---|---|---|
| committer | Ian Lance Taylor <ian@gcc.gnu.org> | 2018-01-09 01:23:08 +0000 |
| commit | 1a2f01efa63036a5104f203a4789e682c0e0915d (patch) | |
| tree | 373e15778dc8295354584e1f86915ae493b604ff /libgo/go/math | |
| parent | 8799df67f2dab88f9fda11739c501780a85575e2 (diff) | |
| download | gcc-1a2f01efa63036a5104f203a4789e682c0e0915d.zip gcc-1a2f01efa63036a5104f203a4789e682c0e0915d.tar.gz gcc-1a2f01efa63036a5104f203a4789e682c0e0915d.tar.bz2 | |
libgo: update to Go1.10beta1
Update the Go library to the 1.10beta1 release.
Requires a few changes to the compiler for modifications to the map
runtime code, and to handle some nowritebarrier cases in the runtime.
Reviewed-on: https://go-review.googlesource.com/86455
gotools/:
* Makefile.am (go_cmd_vet_files): New variable.
(go_cmd_buildid_files, go_cmd_test2json_files): New variables.
(s-zdefaultcc): Change from constants to functions.
(noinst_PROGRAMS): Add vet, buildid, and test2json.
(cgo$(EXEEXT)): Link against $(LIBGOTOOL).
(vet$(EXEEXT)): New target.
(buildid$(EXEEXT)): New target.
(test2json$(EXEEXT)): New target.
(install-exec-local): Install all $(noinst_PROGRAMS).
(uninstall-local): Uninstasll all $(noinst_PROGRAMS).
(check-go-tool): Depend on $(noinst_PROGRAMS). Copy down
objabi.go.
(check-runtime): Depend on $(noinst_PROGRAMS).
(check-cgo-test, check-carchive-test): Likewise.
(check-vet): New target.
(check): Depend on check-vet. Look at cmd_vet-testlog.
(.PHONY): Add check-vet.
* Makefile.in: Rebuild.
From-SVN: r256365
Diffstat (limited to 'libgo/go/math')
37 files changed, 1986 insertions, 255 deletions
diff --git a/libgo/go/math/abs.go b/libgo/go/math/abs.go index 1b7c7c1..d74a090 100644 --- a/libgo/go/math/abs.go +++ b/libgo/go/math/abs.go @@ -18,14 +18,5 @@ func Abs(x float64) float64 { } func abs(x float64) float64 { - // TODO: once golang.org/issue/13095 is fixed, change this to: - // return Float64frombits(Float64bits(x) &^ (1 << 63)) - // But for now, this generates better code and can also be inlined: - if x < 0 { - return -x - } - if x == 0 { - return 0 // return correctly abs(-0) - } - return x + return Float64frombits(Float64bits(x) &^ (1 << 63)) } diff --git a/libgo/go/math/all_test.go b/libgo/go/math/all_test.go index 39a3a49..0412c19 100644 --- a/libgo/go/math/all_test.go +++ b/libgo/go/math/all_test.go @@ -8,6 +8,7 @@ import ( "fmt" . "math" "testing" + "unsafe" ) var vf = []float64{ @@ -210,6 +211,18 @@ var erfc = []float64{ 7.9630075582117758758440411e-01, 1.7806938696800922672994468e+00, } +var erfinv = []float64{ + 4.746037673358033586786350696e-01, + 8.559054432692110956388764172e-01, + -2.45427830571707336251331946e-02, + -4.78116683518973366268905506e-01, + 1.479804430319470983648120853e+00, + 2.654485787128896161882650211e-01, + 5.027444534221520197823192493e-01, + 2.466703532707627818954585670e-01, + 1.632011465103005426240343116e-01, + -1.06672334642196900710000389e+00, +} var exp = []float64{ 1.4533071302642137507696589e+02, 2.2958822575694449002537581e+03, @@ -516,6 +529,18 @@ var remainder = []float64{ 8.734595415957246977711748e-01, 1.314075231424398637614104e+00, } +var round = []float64{ + 5, + 8, + Copysign(0, -1), + -5, + 10, + 3, + 5, + 3, + 2, + -9, +} var signbit = []bool{ false, false, @@ -941,6 +966,40 @@ var erfcSC = []float64{ NaN(), } +var vferfinvSC = []float64{ + 1, + -1, + 0, + Inf(-1), + Inf(1), + NaN(), +} +var erfinvSC = []float64{ + Inf(+1), + Inf(-1), + 0, + NaN(), + NaN(), + NaN(), +} + +var vferfcinvSC = []float64{ + 0, + 2, + 1, + Inf(1), + Inf(-1), + NaN(), +} +var erfcinvSC = []float64{ + Inf(+1), + Inf(-1), + 0, + NaN(), + NaN(), + NaN(), +} + var vfexpSC = []float64{ Inf(-1), -2000, @@ -952,6 +1011,7 @@ var vfexpSC = []float64{ // Issue 18912 1.48852223e+09, 1.4885222e+09, + 1, } var expSC = []float64{ 0, @@ -962,6 +1022,7 @@ var expSC = []float64{ Inf(1), Inf(1), Inf(1), + 2.718281828459045, } var vfexp2SC = []float64{ @@ -1368,6 +1429,8 @@ var vfldexpSC = []fi{ {Inf(-1), 0}, {Inf(-1), -1024}, {NaN(), -1024}, + {10, int(1) << (uint64(unsafe.Sizeof(0)-1) * 8)}, + {10, -(int(1) << (uint64(unsafe.Sizeof(0)-1) * 8))}, } var ldexpSC = []float64{ 0, @@ -1381,6 +1444,8 @@ var ldexpSC = []float64{ Inf(-1), Inf(-1), NaN(), + Inf(1), + 0, } var vflgammaSC = []float64{ @@ -1581,6 +1646,17 @@ var vfpowSC = [][2]float64{ {NaN(), 1}, {NaN(), Pi}, {NaN(), NaN()}, + + // Issue #7394 overflow checks + {2, float64(1 << 32)}, + {2, -float64(1 << 32)}, + {-2, float64(1<<32 + 1)}, + {1 / 2, float64(1 << 45)}, + {1 / 2, -float64(1 << 45)}, + {Nextafter(1, 2), float64(1 << 63)}, + {Nextafter(1, -2), float64(1 << 63)}, + {Nextafter(-1, 2), float64(1 << 63)}, + {Nextafter(-1, -2), float64(1 << 63)}, } var powSC = []float64{ 0, // pow(-Inf, -Pi) @@ -1642,6 +1718,17 @@ var powSC = []float64{ NaN(), // pow(NaN, 1) NaN(), // pow(NaN, +Pi) NaN(), // pow(NaN, NaN) + + // Issue #7394 overflow checks + Inf(1), // pow(2, float64(1 << 32)) + 0, // pow(2, -float64(1 << 32)) + Inf(-1), // pow(-2, float64(1<<32 + 1)) + 0, // pow(1/2, float64(1 << 45)) + Inf(1), // pow(1/2, -float64(1 << 45)) + Inf(1), // pow(Nextafter(1, 2), float64(1 << 63)) + 0, // pow(Nextafter(1, -2), float64(1 << 63)) + 0, // pow(Nextafter(-1, 2), float64(1 << 63)) + Inf(1), // pow(Nextafter(-1, -2), float64(1 << 63)) } var vfpow10SC = []int{ @@ -1680,6 +1767,37 @@ var pow10SC = []float64{ Inf(1), // pow10(MaxInt32) } +var vfroundSC = [][2]float64{ + {0, 0}, + {1.390671161567e-309, 0}, // denormal + {0.49999999999999994, 0}, // 0.5-epsilon + {0.5, 1}, + {0.5000000000000001, 1}, // 0.5+epsilon + {-1.5, -2}, + {-2.5, -3}, + {NaN(), NaN()}, + {Inf(1), Inf(1)}, + {2251799813685249.5, 2251799813685250}, // 1 bit fraction + {2251799813685250.5, 2251799813685251}, + {4503599627370495.5, 4503599627370496}, // 1 bit fraction, rounding to 0 bit fraction + {4503599627370497, 4503599627370497}, // large integer +} +var vfroundEvenSC = [][2]float64{ + {0, 0}, + {1.390671161567e-309, 0}, // denormal + {0.49999999999999994, 0}, // 0.5-epsilon + {0.5, 0}, + {0.5000000000000001, 1}, // 0.5+epsilon + {-1.5, -2}, + {-2.5, -2}, + {NaN(), NaN()}, + {Inf(1), Inf(1)}, + {2251799813685249.5, 2251799813685250}, // 1 bit fraction + {2251799813685250.5, 2251799813685250}, + {4503599627370495.5, 4503599627370496}, // 1 bit fraction, rounding to 0 bit fraction + {4503599627370497, 4503599627370497}, // large integer +} + var vfsignbitSC = []float64{ Inf(-1), Copysign(0, -1), @@ -2093,6 +2211,54 @@ func TestErfc(t *testing.T) { } } +func TestErfinv(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := vf[i] / 10 + if f := Erfinv(a); !veryclose(erfinv[i], f) { + t.Errorf("Erfinv(%g) = %g, want %g", a, f, erfinv[i]) + } + } + for i := 0; i < len(vferfinvSC); i++ { + if f := Erfinv(vferfinvSC[i]); !alike(erfinvSC[i], f) { + t.Errorf("Erfinv(%g) = %g, want %g", vferfinvSC[i], f, erfinvSC[i]) + } + } + for x := -0.9; x <= 0.90; x += 1e-2 { + if f := Erf(Erfinv(x)); !close(x, f) { + t.Errorf("Erf(Erfinv(%g)) = %g, want %g", x, f, x) + } + } + for x := -0.9; x <= 0.90; x += 1e-2 { + if f := Erfinv(Erf(x)); !close(x, f) { + t.Errorf("Erfinv(Erf(%g)) = %g, want %g", x, f, x) + } + } +} + +func TestErfcinv(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := 1.0 - (vf[i] / 10) + if f := Erfcinv(a); !veryclose(erfinv[i], f) { + t.Errorf("Erfcinv(%g) = %g, want %g", a, f, erfinv[i]) + } + } + for i := 0; i < len(vferfcinvSC); i++ { + if f := Erfcinv(vferfcinvSC[i]); !alike(erfcinvSC[i], f) { + t.Errorf("Erfcinv(%g) = %g, want %g", vferfcinvSC[i], f, erfcinvSC[i]) + } + } + for x := 0.1; x <= 1.9; x += 1e-2 { + if f := Erfc(Erfcinv(x)); !close(x, f) { + t.Errorf("Erfc(Erfcinv(%g)) = %g, want %g", x, f, x) + } + } + for x := 0.1; x <= 1.9; x += 1e-2 { + if f := Erfcinv(Erfc(x)); !close(x, f) { + t.Errorf("Erfcinv(Erfc(%g)) = %g, want %g", x, f, x) + } + } +} + func TestExp(t *testing.T) { testExp(t, Exp, "Exp") testExp(t, ExpGo, "ExpGo") @@ -2590,6 +2756,32 @@ func TestRemainder(t *testing.T) { } } +func TestRound(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Round(vf[i]); !alike(round[i], f) { + t.Errorf("Round(%g) = %g, want %g", vf[i], f, round[i]) + } + } + for i := 0; i < len(vfroundSC); i++ { + if f := Round(vfroundSC[i][0]); !alike(vfroundSC[i][1], f) { + t.Errorf("Round(%g) = %g, want %g", vfroundSC[i][0], f, vfroundSC[i][1]) + } + } +} + +func TestRoundToEven(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := RoundToEven(vf[i]); !alike(round[i], f) { + t.Errorf("RoundToEven(%g) = %g, want %g", vf[i], f, round[i]) + } + } + for i := 0; i < len(vfroundEvenSC); i++ { + if f := RoundToEven(vfroundEvenSC[i][0]); !alike(vfroundEvenSC[i][1], f) { + t.Errorf("RoundToEven(%g) = %g, want %g", vfroundEvenSC[i][0], f, vfroundEvenSC[i][1]) + } + } +} + func TestSignbit(t *testing.T) { for i := 0; i < len(vf); i++ { if f := Signbit(vf[i]); signbit[i] != f { @@ -2906,10 +3098,12 @@ func BenchmarkCeil(b *testing.B) { GlobalF = x } +var copysignNeg = -1.0 + func BenchmarkCopysign(b *testing.B) { x := 0.0 for i := 0; i < b.N; i++ { - x = Copysign(.5, -1) + x = Copysign(.5, copysignNeg) } GlobalF = x } @@ -2946,6 +3140,22 @@ func BenchmarkErfc(b *testing.B) { GlobalF = x } +func BenchmarkErfinv(b *testing.B) { + x := 0.0 + for i := 0; i < b.N; i++ { + x = Erfinv(.5) + } + GlobalF = x +} + +func BenchmarkErfcinv(b *testing.B) { + x := 0.0 + for i := 0; i < b.N; i++ { + x = Erfcinv(.5) + } + GlobalF = x +} + func BenchmarkExp(b *testing.B) { x := 0.0 for i := 0; i < b.N; i++ { @@ -2986,10 +3196,12 @@ func BenchmarkExp2Go(b *testing.B) { GlobalF = x } +var absPos = .5 + func BenchmarkAbs(b *testing.B) { x := 0.0 for i := 0; i < b.N; i++ { - x = Abs(.5) + x = Abs(absPos) } GlobalF = x @@ -2998,7 +3210,7 @@ func BenchmarkAbs(b *testing.B) { func BenchmarkDim(b *testing.B) { x := 0.0 for i := 0; i < b.N; i++ { - x = Dim(10, 3) + x = Dim(GlobalF, x) } GlobalF = x } @@ -3221,6 +3433,24 @@ func BenchmarkPow10Neg(b *testing.B) { GlobalF = x } +var roundNeg = float64(-2.5) + +func BenchmarkRound(b *testing.B) { + x := 0.0 + for i := 0; i < b.N; i++ { + x = Round(roundNeg) + } + GlobalF = x +} + +func BenchmarkRoundToEven(b *testing.B) { + x := 0.0 + for i := 0; i < b.N; i++ { + x = RoundToEven(roundNeg) + } + GlobalF = x +} + func BenchmarkRemainder(b *testing.B) { x := 0.0 for i := 0; i < b.N; i++ { @@ -3229,10 +3459,12 @@ func BenchmarkRemainder(b *testing.B) { GlobalF = x } +var signbitPos = 2.5 + func BenchmarkSignbit(b *testing.B) { x := false for i := 0; i < b.N; i++ { - x = Signbit(2.5) + x = Signbit(signbitPos) } GlobalB = x } diff --git a/libgo/go/math/big/calibrate_test.go b/libgo/go/math/big/calibrate_test.go index f69ffbf..2b96e74 100644 --- a/libgo/go/math/big/calibrate_test.go +++ b/libgo/go/math/big/calibrate_test.go @@ -2,13 +2,20 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. +// Calibration used to determine thresholds for using +// different algorithms. Ideally, this would be converted +// to go generate to create thresholds.go + // This file prints execution times for the Mul benchmark // given different Karatsuba thresholds. The result may be // used to manually fine-tune the threshold constant. The // results are somewhat fragile; use repeated runs to get // a clear picture. -// Usage: go test -run=TestCalibrate -calibrate +// Calculates lower and upper thresholds for when basicSqr +// is faster than standard multiplication. + +// Usage: go test -run=TestCalibrate -v -calibrate package big @@ -21,6 +28,27 @@ import ( var calibrate = flag.Bool("calibrate", false, "run calibration test") +func TestCalibrate(t *testing.T) { + if *calibrate { + computeKaratsubaThresholds() + + // compute basicSqrThreshold where overhead becomes negligible + minSqr := computeSqrThreshold(10, 30, 1, 3) + // compute karatsubaSqrThreshold where karatsuba is faster + maxSqr := computeSqrThreshold(300, 500, 10, 3) + if minSqr != 0 { + fmt.Printf("found basicSqrThreshold = %d\n", minSqr) + } else { + fmt.Println("no basicSqrThreshold found") + } + if maxSqr != 0 { + fmt.Printf("found karatsubaSqrThreshold = %d\n", maxSqr) + } else { + fmt.Println("no karatsubaSqrThreshold found") + } + } +} + func karatsubaLoad(b *testing.B) { BenchmarkMul(b) } @@ -34,7 +62,7 @@ func measureKaratsuba(th int) time.Duration { return time.Duration(res.NsPerOp()) } -func computeThresholds() { +func computeKaratsubaThresholds() { fmt.Printf("Multiplication times for varying Karatsuba thresholds\n") fmt.Printf("(run repeatedly for good results)\n") @@ -81,8 +109,56 @@ func computeThresholds() { } } -func TestCalibrate(t *testing.T) { - if *calibrate { - computeThresholds() +func measureBasicSqr(words, nruns int, enable bool) time.Duration { + // more runs for better statistics + initBasicSqr, initKaratsubaSqr := basicSqrThreshold, karatsubaSqrThreshold + + if enable { + // set thresholds to use basicSqr at this number of words + basicSqrThreshold, karatsubaSqrThreshold = words-1, words+1 + } else { + // set thresholds to disable basicSqr for any number of words + basicSqrThreshold, karatsubaSqrThreshold = -1, -1 + } + + var testval int64 + for i := 0; i < nruns; i++ { + res := testing.Benchmark(func(b *testing.B) { benchmarkNatSqr(b, words) }) + testval += res.NsPerOp() + } + testval /= int64(nruns) + + basicSqrThreshold, karatsubaSqrThreshold = initBasicSqr, initKaratsubaSqr + + return time.Duration(testval) +} + +func computeSqrThreshold(from, to, step, nruns int) int { + fmt.Println("Calibrating thresholds for basicSqr via benchmarks of z.mul(x,x)") + fmt.Printf("Looking for a timing difference for x between %d - %d words by %d step\n", from, to, step) + var initPos bool + var threshold int + for i := from; i <= to; i += step { + baseline := measureBasicSqr(i, nruns, false) + testval := measureBasicSqr(i, nruns, true) + pos := baseline > testval + delta := baseline - testval + percent := delta * 100 / baseline + fmt.Printf("words = %3d deltaT = %10s (%4d%%) is basicSqr better: %v", i, delta, percent, pos) + if i == from { + initPos = pos + } + if threshold == 0 && pos != initPos { + threshold = i + fmt.Printf(" threshold found") + } + fmt.Println() + + } + if threshold != 0 { + fmt.Printf("Found threshold = %d between %d - %d\n", threshold, from, to) + } else { + fmt.Printf("Found NO threshold between %d - %d\n", from, to) } + return threshold } diff --git a/libgo/go/math/big/decimal.go b/libgo/go/math/big/decimal.go index 2dfa032..ae9ffb5 100644 --- a/libgo/go/math/big/decimal.go +++ b/libgo/go/math/big/decimal.go @@ -20,7 +20,7 @@ package big // A decimal represents an unsigned floating-point number in decimal representation. -// The value of a non-zero decimal d is d.mant * 10**d.exp with 0.5 <= d.mant < 1, +// The value of a non-zero decimal d is d.mant * 10**d.exp with 0.1 <= d.mant < 1, // with the most-significant mantissa digit at index 0. For the zero decimal, the // mantissa length and exponent are 0. // The zero value for decimal represents a ready-to-use 0.0. diff --git a/libgo/go/math/big/float.go b/libgo/go/math/big/float.go index 7e11f1a..c042854 100644 --- a/libgo/go/math/big/float.go +++ b/libgo/go/math/big/float.go @@ -415,8 +415,9 @@ func (z *Float) round(sbit uint) { // bits > z.prec: mantissa too large => round r := uint(bits - z.prec - 1) // rounding bit position; r >= 0 rbit := z.mant.bit(r) & 1 // rounding bit; be safe and ensure it's a single bit - if sbit == 0 { - // TODO(gri) if rbit != 0 we don't need to compute sbit for some rounding modes (optimization) + // The sticky bit is only needed for rounding ToNearestEven + // or when the rounding bit is zero. Avoid computation otherwise. + if sbit == 0 && (rbit == 0 || z.mode == ToNearestEven) { sbit = z.mant.sticky(r) } sbit &= 1 // be safe and ensure it's a single bit @@ -1310,8 +1311,11 @@ func (z *Float) umul(x, y *Float) { // TODO(gri) Optimize this for the common case. e := int64(x.exp) + int64(y.exp) - z.mant = z.mant.mul(x.mant, y.mant) - + if x == y { + z.mant = z.mant.sqr(x.mant) + } else { + z.mant = z.mant.mul(x.mant, y.mant) + } z.setExpAndRound(e-fnorm(z.mant), 0) } diff --git a/libgo/go/math/big/int.go b/libgo/go/math/big/int.go index 62f7fc5..0eda9cd 100644 --- a/libgo/go/math/big/int.go +++ b/libgo/go/math/big/int.go @@ -153,6 +153,11 @@ func (z *Int) Mul(x, y *Int) *Int { // x * (-y) == -(x * y) // (-x) * y == -(x * y) // (-x) * (-y) == x * y + if x == y { + z.abs = z.abs.sqr(x.abs) + z.neg = false + return z + } z.abs = z.abs.mul(x.abs, y.abs) z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign return z @@ -324,6 +329,16 @@ func (x *Int) Cmp(y *Int) (r int) { return } +// CmpAbs compares the absolute values of x and y and returns: +// +// -1 if |x| < |y| +// 0 if |x| == |y| +// +1 if |x| > |y| +// +func (x *Int) CmpAbs(y *Int) int { + return x.abs.cmp(y.abs) +} + // low32 returns the least significant 32 bits of x. func low32(x nat) uint32 { if len(x) == 0 { @@ -384,16 +399,26 @@ func (x *Int) IsUint64() bool { // ``0x'' or ``0X'' selects base 16; the ``0'' prefix selects base 8, and a // ``0b'' or ``0B'' prefix selects base 2. Otherwise the selected base is 10. // +// For bases <= 36, lower and upper case letters are considered the same: +// The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35. +// For bases > 36, the upper case letters 'A' to 'Z' represent the digit +// values 36 to 61. +// func (z *Int) SetString(s string, base int) (*Int, bool) { - r := strings.NewReader(s) + return z.setFromScanner(strings.NewReader(s), base) +} + +// setFromScanner implements SetString given an io.BytesScanner. +// For documentation see comments of SetString. +func (z *Int) setFromScanner(r io.ByteScanner, base int) (*Int, bool) { if _, _, err := z.scan(r, base); err != nil { return nil, false } - // entire string must have been consumed + // entire content must have been consumed if _, err := r.ReadByte(); err != io.EOF { return nil, false } - return z, true // err == io.EOF => scan consumed all of s + return z, true // err == io.EOF => scan consumed all content of r } // SetBytes interprets buf as the bytes of a big-endian unsigned @@ -447,7 +472,7 @@ func (z *Int) Exp(x, y, m *Int) *Int { // GCD sets z to the greatest common divisor of a and b, which both must // be > 0, and returns z. -// If x and y are not nil, GCD sets x and y such that z = a*x + b*y. +// If x or y are not nil, GCD sets their value such that z = a*x + b*y. // If either a or b is <= 0, GCD sets z = x = y = 0. func (z *Int) GCD(x, y, a, b *Int) *Int { if a.Sign() <= 0 || b.Sign() <= 0 { @@ -461,17 +486,14 @@ func (z *Int) GCD(x, y, a, b *Int) *Int { return z } if x == nil && y == nil { - return z.binaryGCD(a, b) + return z.lehmerGCD(a, b) } A := new(Int).Set(a) B := new(Int).Set(b) X := new(Int) - Y := new(Int).SetInt64(1) - lastX := new(Int).SetInt64(1) - lastY := new(Int) q := new(Int) temp := new(Int) @@ -484,15 +506,8 @@ func (z *Int) GCD(x, y, a, b *Int) *Int { temp.Set(X) X.Mul(X, q) - X.neg = !X.neg - X.Add(X, lastX) + X.Sub(lastX, X) lastX.Set(temp) - - temp.Set(Y) - Y.Mul(Y, q) - Y.neg = !Y.neg - Y.Add(Y, lastY) - lastY.Set(temp) } if x != nil { @@ -500,74 +515,139 @@ func (z *Int) GCD(x, y, a, b *Int) *Int { } if y != nil { - *y = *lastY + // y = (z - a*x)/b + y.Mul(a, lastX) + y.Sub(A, y) + y.Div(y, b) } *z = *A return z } -// binaryGCD sets z to the greatest common divisor of a and b, which both must -// be > 0, and returns z. -// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm B. -func (z *Int) binaryGCD(a, b *Int) *Int { - u := z - v := new(Int) - - // use one Euclidean iteration to ensure that u and v are approx. the same size - switch { - case len(a.abs) > len(b.abs): - // must set v before u since u may be alias for a or b (was issue #11284) - v.Rem(a, b) - u.Set(b) - case len(a.abs) < len(b.abs): - v.Rem(b, a) - u.Set(a) - default: - v.Set(b) - u.Set(a) - } - // a, b must not be used anymore (may be aliases with u) - - // v might be 0 now - if len(v.abs) == 0 { - return u +// lehmerGCD sets z to the greatest common divisor of a and b, +// which both must be > 0, and returns z. +// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L. +// This implementation uses the improved condition by Collins requiring only one +// quotient and avoiding the possibility of single Word overflow. +// See Jebelean, "Improving the multiprecision Euclidean algorithm", +// Design and Implementation of Symbolic Computation Systems, pp 45-58. +func (z *Int) lehmerGCD(a, b *Int) *Int { + // ensure a >= b + if a.abs.cmp(b.abs) < 0 { + a, b = b, a } - // u > 0 && v > 0 - // determine largest k such that u = u' << k, v = v' << k - k := u.abs.trailingZeroBits() - if vk := v.abs.trailingZeroBits(); vk < k { - k = vk - } - u.Rsh(u, k) - v.Rsh(v, k) + // don't destroy incoming values of a and b + B := new(Int).Set(b) // must be set first in case b is an alias of z + A := z.Set(a) - // determine t (we know that u > 0) + // temp variables for multiprecision update t := new(Int) - if u.abs[0]&1 != 0 { - // u is odd - t.Neg(v) - } else { - t.Set(u) - } + r := new(Int) + s := new(Int) + w := new(Int) + + // loop invariant A >= B + for len(B.abs) > 1 { + // initialize the digits + var a1, a2, u0, u1, u2, v0, v1, v2 Word + + m := len(B.abs) // m >= 2 + n := len(A.abs) // n >= m >= 2 + + // extract the top Word of bits from A and B + h := nlz(A.abs[n-1]) + a1 = (A.abs[n-1] << h) | (A.abs[n-2] >> (_W - h)) + // B may have implicit zero words in the high bits if the lengths differ + switch { + case n == m: + a2 = (B.abs[n-1] << h) | (B.abs[n-2] >> (_W - h)) + case n == m+1: + a2 = (B.abs[n-2] >> (_W - h)) + default: + a2 = 0 + } + + // Since we are calculating with full words to avoid overflow, + // we use 'even' to track the sign of the cosequences. + // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 + // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 + // The first iteration starts with k=1 (odd). + even := false + // variables to track the cosequences + u0, u1, u2 = 0, 1, 0 + v0, v1, v2 = 0, 0, 1 + + // Calculate the quotient and cosequences using Collins' stopping condition. + // Note that overflow of a Word is not possible when computing the remainder + // sequence and cosequences since the cosequence size is bounded by the input size. + // See section 4.2 of Jebelean for details. + for a2 >= v2 && a1-a2 >= v1+v2 { + q := a1 / a2 + a1, a2 = a2, a1-q*a2 + u0, u1, u2 = u1, u2, u1+q*u2 + v0, v1, v2 = v1, v2, v1+q*v2 + even = !even + } + + // multiprecision step + if v0 != 0 { + // simulate the effect of the single precision steps using the cosequences + // A = u0*A + v0*B + // B = u1*A + v1*B + + t.abs = t.abs.setWord(u0) + s.abs = s.abs.setWord(v0) + t.neg = !even + s.neg = even + + t.Mul(A, t) + s.Mul(B, s) + + r.abs = r.abs.setWord(u1) + w.abs = w.abs.setWord(v1) + r.neg = even + w.neg = !even + + r.Mul(A, r) + w.Mul(B, w) + + A.Add(t, s) + B.Add(r, w) - for len(t.abs) > 0 { - // reduce t - t.Rsh(t, t.abs.trailingZeroBits()) - if t.neg { - v, t = t, v - v.neg = len(v.abs) > 0 && !v.neg // 0 has no sign } else { - u, t = t, u + // single-digit calculations failed to simluate any quotients + // do a standard Euclidean step + t.Rem(A, B) + A, B, t = B, t, A } - t.Sub(u, v) } - return z.Lsh(u, k) + if len(B.abs) > 0 { + // standard Euclidean algorithm base case for B a single Word + if len(A.abs) > 1 { + // A is longer than a single Word + t.Rem(A, B) + A, B, t = B, t, A + } + if len(B.abs) > 0 { + // A and B are both a single Word + a1, a2 := A.abs[0], B.abs[0] + for a2 != 0 { + a1, a2 = a2, a1%a2 + } + A.abs[0] = a1 + } + } + *z = *A + return z } // Rand sets z to a pseudo-random number in [0, n) and returns z. +// +// As this uses the math/rand package, it must not be used for +// security-sensitive work. Use crypto/rand.Int instead. func (z *Int) Rand(rnd *rand.Rand, n *Int) *Int { z.neg = false if n.neg || len(n.abs) == 0 { @@ -659,10 +739,9 @@ func Jacobi(x, y *Int) int { // to calculate the square root of any quadratic residue mod p quickly for 3 // mod 4 primes. func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { - z.Set(p) // z = p - z.Add(z, intOne) // z = p + 1 - z.Rsh(z, 2) // z = (p + 1) / 4 - z.Exp(x, z, p) // z = x^z mod p + e := new(Int).Add(p, intOne) // e = p + 1 + e.Rsh(e, 2) // e = (p + 1) / 4 + z.Exp(x, e, p) // z = x^e mod p return z } diff --git a/libgo/go/math/big/int_test.go b/libgo/go/math/big/int_test.go index 42e810b..270fec6 100644 --- a/libgo/go/math/big/int_test.go +++ b/libgo/go/math/big/int_test.go @@ -7,6 +7,7 @@ package big import ( "bytes" "encoding/hex" + "fmt" "math/rand" "strconv" "strings" @@ -674,6 +675,37 @@ func checkGcd(aBytes, bBytes []byte) bool { return x.Cmp(d) == 0 } +// euclidGCD is a reference implementation of Euclid's GCD +// algorithm for testing against optimized algorithms. +// Requirements: a, b > 0 +func euclidGCD(a, b *Int) *Int { + + A := new(Int).Set(a) + B := new(Int).Set(b) + t := new(Int) + + for len(B.abs) > 0 { + // A, B = B, A % B + t.Rem(A, B) + A, B, t = B, t, A + } + return A +} + +func checkLehmerGcd(aBytes, bBytes []byte) bool { + a := new(Int).SetBytes(aBytes) + b := new(Int).SetBytes(bBytes) + + if a.Sign() <= 0 || b.Sign() <= 0 { + return true // can only test positive arguments + } + + d := new(Int).lehmerGCD(a, b) + d0 := euclidGCD(a, b) + + return d.Cmp(d0) == 0 +} + var gcdTests = []struct { d, x, y, a, b string }{ @@ -707,38 +739,28 @@ func testGcd(t *testing.T, d, x, y, a, b *Int) { D := new(Int).GCD(X, Y, a, b) if D.Cmp(d) != 0 { - t.Errorf("GCD(%s, %s): got d = %s, want %s", a, b, D, d) + t.Errorf("GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, D, d) } if x != nil && X.Cmp(x) != 0 { - t.Errorf("GCD(%s, %s): got x = %s, want %s", a, b, X, x) + t.Errorf("GCD(%s, %s, %s, %s): got x = %s, want %s", x, y, a, b, X, x) } if y != nil && Y.Cmp(y) != 0 { - t.Errorf("GCD(%s, %s): got y = %s, want %s", a, b, Y, y) - } - - // binaryGCD requires a > 0 && b > 0 - if a.Sign() <= 0 || b.Sign() <= 0 { - return - } - - D.binaryGCD(a, b) - if D.Cmp(d) != 0 { - t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, D, d) + t.Errorf("GCD(%s, %s, %s, %s): got y = %s, want %s", x, y, a, b, Y, y) } // check results in presence of aliasing (issue #11284) a2 := new(Int).Set(a) b2 := new(Int).Set(b) - a2.binaryGCD(a2, b2) // result is same as 1st argument + a2.GCD(X, Y, a2, b2) // result is same as 1st argument if a2.Cmp(d) != 0 { - t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, a2, d) + t.Errorf("aliased z = a GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, a2, d) } a2 = new(Int).Set(a) b2 = new(Int).Set(b) - b2.binaryGCD(a2, b2) // result is same as 2nd argument + b2.GCD(X, Y, a2, b2) // result is same as 2nd argument if b2.Cmp(d) != 0 { - t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, b2, d) + t.Errorf("aliased z = b GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, b2, d) } } @@ -759,6 +781,10 @@ func TestGcd(t *testing.T) { if err := quick.Check(checkGcd, nil); err != nil { t.Error(err) } + + if err := quick.Check(checkLehmerGcd, nil); err != nil { + t.Error(err) + } } type intShiftTest struct { @@ -903,6 +929,63 @@ func TestLshRsh(t *testing.T) { } } +// Entries must be sorted by value in ascending order. +var cmpAbsTests = []string{ + "0", + "1", + "2", + "10", + "10000000", + "2783678367462374683678456387645876387564783686583485", + "2783678367462374683678456387645876387564783686583486", + "32957394867987420967976567076075976570670947609750670956097509670576075067076027578341538", +} + +func TestCmpAbs(t *testing.T) { + values := make([]*Int, len(cmpAbsTests)) + var prev *Int + for i, s := range cmpAbsTests { + x, ok := new(Int).SetString(s, 0) + if !ok { + t.Fatalf("SetString(%s, 0) failed", s) + } + if prev != nil && prev.Cmp(x) >= 0 { + t.Fatal("cmpAbsTests entries not sorted in ascending order") + } + values[i] = x + prev = x + } + + for i, x := range values { + for j, y := range values { + // try all combinations of signs for x, y + for k := 0; k < 4; k++ { + var a, b Int + a.Set(x) + b.Set(y) + if k&1 != 0 { + a.Neg(&a) + } + if k&2 != 0 { + b.Neg(&b) + } + + got := a.CmpAbs(&b) + want := 0 + switch { + case i > j: + want = 1 + case i < j: + want = -1 + } + if got != want { + t.Errorf("absCmp |%s|, |%s|: got %d; want %d", &a, &b, got, want) + } + } + } + } +} + var int64Tests = []string{ // int64 "0", @@ -1383,6 +1466,12 @@ func testModSqrt(t *testing.T, elt, mod, sq, sqrt *Int) bool { t.Errorf("ModSqrt returned inconsistent value %s", z) } + // test x aliasing z + z = sqrtChk.ModSqrt(sqrtChk.Set(sq), mod) + if z != &sqrtChk || z.Cmp(sqrt) != 0 { + t.Errorf("ModSqrt returned inconsistent value %s", z) + } + // make sure we actually got a square root if sqrt.Cmp(elt) == 0 { return true // we found the "desired" square root @@ -1536,6 +1625,26 @@ func TestSqrt(t *testing.T) { } } +// We can't test this together with the other Exp tests above because +// it requires a different receiver setup. +func TestIssue22830(t *testing.T) { + one := new(Int).SetInt64(1) + base, _ := new(Int).SetString("84555555300000000000", 10) + mod, _ := new(Int).SetString("66666670001111111111", 10) + want, _ := new(Int).SetString("17888885298888888889", 10) + + var tests = []int64{ + 0, 1, -1, + } + + for _, n := range tests { + m := NewInt(n) + if got := m.Exp(base, one, mod); got.Cmp(want) != 0 { + t.Errorf("(%v).Exp(%s, 1, %s) = %s, want %s", n, base, mod, got, want) + } + } +} + func BenchmarkSqrt(b *testing.B) { n, _ := new(Int).SetString("1"+strings.Repeat("0", 1001), 10) b.ResetTimer() @@ -1544,3 +1653,24 @@ func BenchmarkSqrt(b *testing.B) { t.Sqrt(n) } } + +func benchmarkIntSqr(b *testing.B, nwords int) { + x := new(Int) + x.abs = rndNat(nwords) + t := new(Int) + b.ResetTimer() + for i := 0; i < b.N; i++ { + t.Mul(x, x) + } +} + +func BenchmarkIntSqr(b *testing.B) { + for _, n := range sqrBenchSizes { + if isRaceBuilder && n > 1e3 { + continue + } + b.Run(fmt.Sprintf("%d", n), func(b *testing.B) { + benchmarkIntSqr(b, n) + }) + } +} diff --git a/libgo/go/math/big/intconv.go b/libgo/go/math/big/intconv.go index 91a62ce..6cca827 100644 --- a/libgo/go/math/big/intconv.go +++ b/libgo/go/math/big/intconv.go @@ -12,16 +12,11 @@ import ( "io" ) -// TODO(gri) Should rename itoa to utoa (there's no sign). That -// would permit the introduction of itoa which is like utoa but -// reserves a byte for a possible sign that's passed in. That -// would permit Int.Text to be implemented w/o the need for -// string copy if the number is negative. - // Text returns the string representation of x in the given base. -// Base must be between 2 and 36, inclusive. The result uses the -// lower-case letters 'a' to 'z' for digit values >= 10. No base -// prefix (such as "0x") is added to the string. +// Base must be between 2 and 62, inclusive. The result uses the +// lower-case letters 'a' to 'z' for digit values 10 to 35, and +// the upper-case letters 'A' to 'Z' for digit values 36 to 61. +// No prefix (such as "0x") is added to the string. func (x *Int) Text(base int) string { if x == nil { return "<nil>" diff --git a/libgo/go/math/big/intconv_test.go b/libgo/go/math/big/intconv_test.go index 5142081..2e01ee3 100644 --- a/libgo/go/math/big/intconv_test.go +++ b/libgo/go/math/big/intconv_test.go @@ -56,6 +56,10 @@ var stringTests = []struct { {"-0b111", "-7", 0, -7, true}, {"0b1001010111", "599", 0, 0x257, true}, {"1001010111", "1001010111", 2, 0x257, true}, + {"A", "a", 36, 10, true}, + {"A", "A", 37, 36, true}, + {"ABCXYZ", "abcxyz", 36, 623741435, true}, + {"ABCXYZ", "ABCXYZ", 62, 33536793425, true}, } func TestIntText(t *testing.T) { @@ -135,8 +139,16 @@ func TestGetString(t *testing.T) { } } - if got := fmt.Sprintf(format(test.base), z); got != test.out { - t.Errorf("#%db got %s; want %s", i, got, test.out) + f := format(test.base) + got := fmt.Sprintf(f, z) + if f == "%d" { + if got != fmt.Sprintf("%d", test.val) { + t.Errorf("#%db got %s; want %d", i, got, test.val) + } + } else { + if got != test.out { + t.Errorf("#%dc got %s; want %s", i, got, test.out) + } } } } diff --git a/libgo/go/math/big/intmarsh.go b/libgo/go/math/big/intmarsh.go index ee1e414..c1422e2 100644 --- a/libgo/go/math/big/intmarsh.go +++ b/libgo/go/math/big/intmarsh.go @@ -6,7 +6,10 @@ package big -import "fmt" +import ( + "bytes" + "fmt" +) // Gob codec version. Permits backward-compatible changes to the encoding. const intGobVersion byte = 1 @@ -52,8 +55,7 @@ func (x *Int) MarshalText() (text []byte, err error) { // UnmarshalText implements the encoding.TextUnmarshaler interface. func (z *Int) UnmarshalText(text []byte) error { - // TODO(gri): get rid of the []byte/string conversion - if _, ok := z.SetString(string(text), 0); !ok { + if _, ok := z.setFromScanner(bytes.NewReader(text), 0); !ok { return fmt.Errorf("math/big: cannot unmarshal %q into a *big.Int", text) } return nil diff --git a/libgo/go/math/big/nat.go b/libgo/go/math/big/nat.go index 889eacb..3bb818f 100644 --- a/libgo/go/math/big/nat.go +++ b/libgo/go/math/big/nat.go @@ -249,7 +249,7 @@ func karatsubaSub(z, x nat, n int) { // Operands that are shorter than karatsubaThreshold are multiplied using // "grade school" multiplication; for longer operands the Karatsuba algorithm // is used. -var karatsubaThreshold int = 40 // computed by calibrate.go +var karatsubaThreshold = 40 // computed by calibrate_test.go // karatsuba multiplies x and y and leaves the result in z. // Both x and y must have the same length n and n must be a @@ -473,6 +473,61 @@ func (z nat) mul(x, y nat) nat { return z.norm() } +// basicSqr sets z = x*x and is asymptotically faster than basicMul +// by about a factor of 2, but slower for small arguments due to overhead. +// Requirements: len(x) > 0, len(z) >= 2*len(x) +// The (non-normalized) result is placed in z[0 : 2 * len(x)]. +func basicSqr(z, x nat) { + n := len(x) + t := make(nat, 2*n) // temporary variable to hold the products + z[1], z[0] = mulWW(x[0], x[0]) // the initial square + for i := 1; i < n; i++ { + d := x[i] + // z collects the squares x[i] * x[i] + z[2*i+1], z[2*i] = mulWW(d, d) + // t collects the products x[i] * x[j] where j < i + t[2*i] = addMulVVW(t[i:2*i], x[0:i], d) + } + t[2*n-1] = shlVU(t[1:2*n-1], t[1:2*n-1], 1) // double the j < i products + addVV(z, z, t) // combine the result +} + +// Operands that are shorter than basicSqrThreshold are squared using +// "grade school" multiplication; for operands longer than karatsubaSqrThreshold +// the Karatsuba algorithm is used. +var basicSqrThreshold = 20 // computed by calibrate_test.go +var karatsubaSqrThreshold = 400 // computed by calibrate_test.go + +// z = x*x +func (z nat) sqr(x nat) nat { + n := len(x) + switch { + case n == 0: + return z[:0] + case n == 1: + d := x[0] + z = z.make(2) + z[1], z[0] = mulWW(d, d) + return z.norm() + } + + if alias(z, x) { + z = nil // z is an alias for x - cannot reuse + } + z = z.make(2 * n) + + if n < basicSqrThreshold { + basicMul(z, x, x) + return z.norm() + } + if n < karatsubaSqrThreshold { + basicSqr(z, x) + return z.norm() + } + + return z.mul(x, x) +} + // mulRange computes the product of all the unsigned integers in the // range [a, b] inclusively. If a > b (empty range), the result is 1. func (z nat) mulRange(a, b uint64) nat { @@ -566,8 +621,8 @@ func (z nat) divLarge(u, uIn, v nat) (q, r nat) { // determine if z can be reused // TODO(gri) should find a better solution - this if statement // is very costly (see e.g. time pidigits -s -n 10000) - if alias(z, uIn) || alias(z, v) { - z = nil // z is an alias for uIn or v - cannot reuse + if alias(z, u) || alias(z, uIn) || alias(z, v) { + z = nil // z is an alias for u or uIn or v - cannot reuse } q = z.make(m + 1) @@ -936,7 +991,7 @@ func (z nat) expNN(x, y, m nat) nat { // otherwise the arguments would alias. var zz, r nat for j := 0; j < w; j++ { - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz if v&mask != 0 { @@ -956,7 +1011,7 @@ func (z nat) expNN(x, y, m nat) nat { v = y[i] for j := 0; j < _W; j++ { - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz if v&mask != 0 { @@ -989,7 +1044,7 @@ func (z nat) expNNWindowed(x, y, m nat) nat { powers[1] = x for i := 2; i < 1<<n; i += 2 { p2, p, p1 := &powers[i/2], &powers[i], &powers[i+1] - *p = p.mul(*p2, *p2) + *p = p.sqr(*p2) zz, r = zz.div(r, *p, m) *p, r = r, *p *p1 = p1.mul(*p, x) @@ -1006,22 +1061,22 @@ func (z nat) expNNWindowed(x, y, m nat) nat { // Unrolled loop for significant performance // gain. Use go test -bench=".*" in crypto/rsa // to check performance before making changes. - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz zz, r = zz.div(r, z, m) z, r = r, z - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz zz, r = zz.div(r, z, m) z, r = r, z - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz zz, r = zz.div(r, z, m) z, r = r, z - zz = zz.mul(z, z) + zz = zz.sqr(z) zz, z = z, zz zz, r = zz.div(r, z, m) z, r = r, z diff --git a/libgo/go/math/big/nat_test.go b/libgo/go/math/big/nat_test.go index 200a247..c25cdf0 100644 --- a/libgo/go/math/big/nat_test.go +++ b/libgo/go/math/big/nat_test.go @@ -619,3 +619,49 @@ func TestSticky(t *testing.T) { } } } + +func testBasicSqr(t *testing.T, x nat) { + got := make(nat, 2*len(x)) + want := make(nat, 2*len(x)) + basicSqr(got, x) + basicMul(want, x, x) + if got.cmp(want) != 0 { + t.Errorf("basicSqr(%v), got %v, want %v", x, got, want) + } +} + +func TestBasicSqr(t *testing.T) { + for _, a := range prodNN { + if a.x != nil { + testBasicSqr(t, a.x) + } + if a.y != nil { + testBasicSqr(t, a.y) + } + if a.z != nil { + testBasicSqr(t, a.z) + } + } +} + +func benchmarkNatSqr(b *testing.B, nwords int) { + x := rndNat(nwords) + var z nat + b.ResetTimer() + for i := 0; i < b.N; i++ { + z.sqr(x) + } +} + +var sqrBenchSizes = []int{1, 2, 3, 5, 8, 10, 20, 30, 50, 80, 100, 200, 300, 500, 800, 1000} + +func BenchmarkNatSqr(b *testing.B) { + for _, n := range sqrBenchSizes { + if isRaceBuilder && n > 1e3 { + continue + } + b.Run(fmt.Sprintf("%d", n), func(b *testing.B) { + benchmarkNatSqr(b, n) + }) + } +} diff --git a/libgo/go/math/big/natconv.go b/libgo/go/math/big/natconv.go index 25a345e..21ccbd6 100644 --- a/libgo/go/math/big/natconv.go +++ b/libgo/go/math/big/natconv.go @@ -15,13 +15,14 @@ import ( "sync" ) -const digits = "0123456789abcdefghijklmnopqrstuvwxyz" +const digits = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ" -// Note: MaxBase = len(digits), but it must remain a rune constant +// Note: MaxBase = len(digits), but it must remain an untyped rune constant // for API compatibility. // MaxBase is the largest number base accepted for string conversions. -const MaxBase = 'z' - 'a' + 10 + 1 +const MaxBase = 10 + ('z' - 'a' + 1) + ('Z' - 'A' + 1) +const maxBaseSmall = 10 + ('z' - 'a' + 1) // maxPow returns (b**n, n) such that b**n is the largest power b**n <= _M. // For instance maxPow(10) == (1e19, 19) for 19 decimal digits in a 64bit Word. @@ -59,11 +60,11 @@ func pow(x Word, n int) (p Word) { // It returns the corresponding natural number res, the actual base b, // a digit count, and a read or syntax error err, if any. // -// number = [ prefix ] mantissa . -// prefix = "0" [ "x" | "X" | "b" | "B" ] . -// mantissa = digits | digits "." [ digits ] | "." digits . -// digits = digit { digit } . -// digit = "0" ... "9" | "a" ... "z" | "A" ... "Z" . +// number = [ prefix ] mantissa . +// prefix = "0" [ "x" | "X" | "b" | "B" ] . +// mantissa = digits | digits "." [ digits ] | "." digits . +// digits = digit { digit } . +// digit = "0" ... "9" | "a" ... "z" | "A" ... "Z" . // // Unless fracOk is set, the base argument must be 0 or a value between // 2 and MaxBase. If fracOk is set, the base argument must be one of @@ -80,6 +81,11 @@ func pow(x Word, n int) (p Word) { // is permitted. The result value is computed as if there were no period // present; and the count value is used to determine the fractional part. // +// For bases <= 36, lower and upper case letters are considered the same: +// The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35. +// For bases > 36, the upper case letters 'A' to 'Z' represent the digit +// values 36 to 61. +// // A result digit count > 0 corresponds to the number of (non-prefix) digits // parsed. A digit count <= 0 indicates the presence of a period (if fracOk // is set, only), and -count is the number of fractional digits found. @@ -173,7 +179,11 @@ func (z nat) scan(r io.ByteScanner, base int, fracOk bool) (res nat, b, count in case 'a' <= ch && ch <= 'z': d1 = Word(ch - 'a' + 10) case 'A' <= ch && ch <= 'Z': - d1 = Word(ch - 'A' + 10) + if b <= maxBaseSmall { + d1 = Word(ch - 'A' + 10) + } else { + d1 = Word(ch - 'A' + maxBaseSmall) + } default: d1 = MaxBase + 1 } @@ -469,7 +479,7 @@ func divisors(m int, b Word, ndigits int, bb Word) []divisor { table[0].bbb = nat(nil).expWW(bb, Word(leafSize)) table[0].ndigits = ndigits * leafSize } else { - table[i].bbb = nat(nil).mul(table[i-1].bbb, table[i-1].bbb) + table[i].bbb = nat(nil).sqr(table[i-1].bbb) table[i].ndigits = 2 * table[i-1].ndigits } diff --git a/libgo/go/math/big/natconv_test.go b/libgo/go/math/big/natconv_test.go index 898a39f..9f38bd9 100644 --- a/libgo/go/math/big/natconv_test.go +++ b/libgo/go/math/big/natconv_test.go @@ -13,6 +13,12 @@ import ( "testing" ) +func TestMaxBase(t *testing.T) { + if MaxBase != len(digits) { + t.Fatalf("%d != %d", MaxBase, len(digits)) + } +} + // log2 computes the integer binary logarithm of x. // The result is the integer n for which 2^n <= x < 2^(n+1). // If x == 0, the result is -1. @@ -61,6 +67,7 @@ var strTests = []struct { {nat{0xdeadbeef}, 16, "deadbeef"}, {nat{0x229be7}, 17, "1a2b3c"}, {nat{0x309663e6}, 32, "o9cov6"}, + {nat{0x309663e6}, 62, "TakXI"}, } func TestString(t *testing.T) { @@ -110,6 +117,7 @@ var natScanTests = []struct { {s: "?"}, {base: 10}, {base: 36}, + {base: 62}, {s: "?", base: 10}, {s: "0x"}, {s: "345", base: 2}, @@ -124,6 +132,7 @@ var natScanTests = []struct { {"0", 0, false, nil, 10, 1, true, 0}, {"0", 10, false, nil, 10, 1, true, 0}, {"0", 36, false, nil, 36, 1, true, 0}, + {"0", 62, false, nil, 62, 1, true, 0}, {"1", 0, false, nat{1}, 10, 1, true, 0}, {"1", 10, false, nat{1}, 10, 1, true, 0}, {"0 ", 0, false, nil, 10, 1, true, ' '}, @@ -135,8 +144,11 @@ var natScanTests = []struct { {"03271", 0, false, nat{03271}, 8, 4, true, 0}, {"10ab", 0, false, nat{10}, 10, 2, true, 'a'}, {"1234567890", 0, false, nat{1234567890}, 10, 10, true, 0}, + {"A", 36, false, nat{10}, 36, 1, true, 0}, + {"A", 37, false, nat{36}, 37, 1, true, 0}, {"xyz", 36, false, nat{(33*36+34)*36 + 35}, 36, 3, true, 0}, - {"xyz?", 36, false, nat{(33*36+34)*36 + 35}, 36, 3, true, '?'}, + {"XYZ?", 36, false, nat{(33*36+34)*36 + 35}, 36, 3, true, '?'}, + {"XYZ?", 62, false, nat{(59*62+60)*62 + 61}, 62, 3, true, '?'}, {"0x", 16, false, nil, 16, 1, true, 'x'}, {"0xdeadbeef", 0, false, nat{0xdeadbeef}, 16, 8, true, 0}, {"0XDEADBEEF", 0, false, nat{0xdeadbeef}, 16, 8, true, 0}, diff --git a/libgo/go/math/big/prime.go b/libgo/go/math/big/prime.go index 3e9690e..848affb 100644 --- a/libgo/go/math/big/prime.go +++ b/libgo/go/math/big/prime.go @@ -108,7 +108,7 @@ NextRandom: continue } for j := uint(1); j < k; j++ { - y = y.mul(y, y) + y = y.sqr(y) quotient, y = quotient.div(y, y, n) if y.cmp(nm1) == 0 { continue NextRandom @@ -194,7 +194,7 @@ func (n nat) probablyPrimeLucas() bool { // If n is a non-square we expect to find a d in just a few attempts on average. // After 40 attempts, take a moment to check if n is indeed a square. t1 = t1.sqrt(n) - t1 = t1.mul(t1, t1) + t1 = t1.sqr(t1) if t1.cmp(n) == 0 { return false } @@ -259,7 +259,7 @@ func (n nat) probablyPrimeLucas() bool { t1 = t1.sub(t1, natP) t2, vk = t2.div(vk, t1, n) // V(k'+1) = V(2k+2) = V(k+1)² - 2. - t1 = t1.mul(vk1, vk1) + t1 = t1.sqr(vk1) t1 = t1.add(t1, nm2) t2, vk1 = t2.div(vk1, t1, n) } else { @@ -270,7 +270,7 @@ func (n nat) probablyPrimeLucas() bool { t1 = t1.sub(t1, natP) t2, vk1 = t2.div(vk1, t1, n) // V(k') = V(2k) = V(k)² - 2 - t1 = t1.mul(vk, vk) + t1 = t1.sqr(vk) t1 = t1.add(t1, nm2) t2, vk = t2.div(vk, t1, n) } @@ -312,7 +312,7 @@ func (n nat) probablyPrimeLucas() bool { } // k' = 2k // V(k') = V(2k) = V(k)² - 2 - t1 = t1.mul(vk, vk) + t1 = t1.sqr(vk) t1 = t1.sub(t1, natTwo) t2, vk = t2.div(vk, t1, n) } diff --git a/libgo/go/math/big/rat.go b/libgo/go/math/big/rat.go index 56ce33d..b33fc69 100644 --- a/libgo/go/math/big/rat.go +++ b/libgo/go/math/big/rat.go @@ -422,7 +422,7 @@ func (z *Rat) norm() *Rat { neg := z.a.neg z.a.neg = false z.b.neg = false - if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 { + if f := NewInt(0).lehmerGCD(&z.a, &z.b); f.Cmp(intOne) != 0 { z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs) z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs) if z.b.abs.cmp(natOne) == 0 { @@ -490,6 +490,13 @@ func (z *Rat) Sub(x, y *Rat) *Rat { // Mul sets z to the product x*y and returns z. func (z *Rat) Mul(x, y *Rat) *Rat { + if x == y { + // a squared Rat is positive and can't be reduced + z.a.neg = false + z.a.abs = z.a.abs.sqr(x.a.abs) + z.b.abs = z.b.abs.sqr(x.b.abs) + return z + } z.a.Mul(&x.a, &y.a) z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs) return z.norm() diff --git a/libgo/go/math/big/ratconv.go b/libgo/go/math/big/ratconv.go index 4bc6ef7..157d8d0 100644 --- a/libgo/go/math/big/ratconv.go +++ b/libgo/go/math/big/ratconv.go @@ -202,6 +202,11 @@ func scanExponent(r io.ByteScanner, binExpOk bool) (exp int64, base int, err err // String returns a string representation of x in the form "a/b" (even if b == 1). func (x *Rat) String() string { + return string(x.marshal()) +} + +// marshal implements String returning a slice of bytes +func (x *Rat) marshal() []byte { var buf []byte buf = x.a.Append(buf, 10) buf = append(buf, '/') @@ -210,7 +215,7 @@ func (x *Rat) String() string { } else { buf = append(buf, '1') } - return string(buf) + return buf } // RatString returns a string representation of x in the form "a/b" if b != 1, diff --git a/libgo/go/math/big/ratmarsh.go b/libgo/go/math/big/ratmarsh.go index b82e8d4..fbc7b60 100644 --- a/libgo/go/math/big/ratmarsh.go +++ b/libgo/go/math/big/ratmarsh.go @@ -59,8 +59,10 @@ func (z *Rat) GobDecode(buf []byte) error { // MarshalText implements the encoding.TextMarshaler interface. func (x *Rat) MarshalText() (text []byte, err error) { - // TODO(gri): get rid of the []byte/string conversion - return []byte(x.RatString()), nil + if x.IsInt() { + return x.a.MarshalText() + } + return x.marshal(), nil } // UnmarshalText implements the encoding.TextUnmarshaler interface. diff --git a/libgo/go/math/big/sqrt.go b/libgo/go/math/big/sqrt.go new file mode 100644 index 0000000..00433cf --- /dev/null +++ b/libgo/go/math/big/sqrt.go @@ -0,0 +1,151 @@ +// Copyright 2017 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 big + +import "math" + +var ( + half = NewFloat(0.5) + two = NewFloat(2.0) + three = NewFloat(3.0) +) + +// Sqrt sets z to the rounded square root of x, and returns it. +// +// If z's precision is 0, it is changed to x's precision before the +// operation. Rounding is performed according to z's precision and +// rounding mode. +// +// The function panics if z < 0. The value of z is undefined in that +// case. +func (z *Float) Sqrt(x *Float) *Float { + if debugFloat { + x.validate() + } + + if z.prec == 0 { + z.prec = x.prec + } + + if x.Sign() == -1 { + // following IEEE754-2008 (section 7.2) + panic(ErrNaN{"square root of negative operand"}) + } + + // handle ±0 and +∞ + if x.form != finite { + z.acc = Exact + z.form = x.form + z.neg = x.neg // IEEE754-2008 requires √±0 = ±0 + return z + } + + // MantExp sets the argument's precision to the receiver's, and + // when z.prec > x.prec this will lower z.prec. Restore it after + // the MantExp call. + prec := z.prec + b := x.MantExp(z) + z.prec = prec + + // Compute √(z·2**b) as + // √( z)·2**(½b) if b is even + // √(2z)·2**(⌊½b⌋) if b > 0 is odd + // √(½z)·2**(⌈½b⌉) if b < 0 is odd + switch b % 2 { + case 0: + // nothing to do + case 1: + z.Mul(two, z) + case -1: + z.Mul(half, z) + } + // 0.25 <= z < 2.0 + + // Solving x² - z = 0 directly requires a Quo call, but it's + // faster for small precisions. + // + // Solving 1/x² - z = 0 avoids the Quo call and is much faster for + // high precisions. + // + // 128bit precision is an empirically chosen threshold. + if z.prec <= 128 { + z.sqrtDirect(z) + } else { + z.sqrtInverse(z) + } + + // re-attach halved exponent + return z.SetMantExp(z, b/2) +} + +// Compute √x (up to prec 128) by solving +// t² - x = 0 +// for t, starting with a 53 bits precision guess from math.Sqrt and +// then using at most two iterations of Newton's method. +func (z *Float) sqrtDirect(x *Float) { + // let + // f(t) = t² - x + // then + // g(t) = f(t)/f'(t) = ½(t² - x)/t + // and the next guess is given by + // t2 = t - g(t) = ½(t² + x)/t + u := new(Float) + ng := func(t *Float) *Float { + u.prec = t.prec + u.Mul(t, t) // u = t² + u.Add(u, x) // = t² + x + u.Mul(half, u) // = ½(t² + x) + return t.Quo(u, t) // = ½(t² + x)/t + } + + xf, _ := x.Float64() + sq := NewFloat(math.Sqrt(xf)) + + switch { + case z.prec > 128: + panic("sqrtDirect: only for z.prec <= 128") + case z.prec > 64: + sq.prec *= 2 + sq = ng(sq) + fallthrough + default: + sq.prec *= 2 + sq = ng(sq) + } + + z.Set(sq) +} + +// Compute √x (to z.prec precision) by solving +// 1/t² - x = 0 +// for t (using Newton's method), and then inverting. +func (z *Float) sqrtInverse(x *Float) { + // let + // f(t) = 1/t² - x + // then + // g(t) = f(t)/f'(t) = -½t(1 - xt²) + // and the next guess is given by + // t2 = t - g(t) = ½t(3 - xt²) + u := new(Float) + ng := func(t *Float) *Float { + u.prec = t.prec + u.Mul(t, t) // u = t² + u.Mul(x, u) // = xt² + u.Sub(three, u) // = 3 - xt² + u.Mul(t, u) // = t(3 - xt²) + return t.Mul(half, u) // = ½t(3 - xt²) + } + + xf, _ := x.Float64() + sqi := NewFloat(1 / math.Sqrt(xf)) + for prec := z.prec + 32; sqi.prec < prec; { + sqi.prec *= 2 + sqi = ng(sqi) + } + // sqi = 1/√x + + // x/√x = √x + z.Mul(x, sqi) +} diff --git a/libgo/go/math/big/sqrt_test.go b/libgo/go/math/big/sqrt_test.go new file mode 100644 index 0000000..86595a5 --- /dev/null +++ b/libgo/go/math/big/sqrt_test.go @@ -0,0 +1,131 @@ +// Copyright 2017 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 big + +import ( + "fmt" + "math" + "math/rand" + "runtime" + "testing" +) + +// TestFloatSqrt64 tests that Float.Sqrt of numbers with 53bit mantissa +// behaves like float math.Sqrt. +func TestFloatSqrt64(t *testing.T) { + // This test fails for gccgo on 386 with a one ULP difference, + // presumably due to the use of extended precision floating + // point. + if runtime.Compiler == "gccgo" && runtime.GOARCH == "386" { + t.Skip("skipping on gccgo for 386; gets a one ULP difference") + } + + for i := 0; i < 1e5; i++ { + r := rand.Float64() + + got := new(Float).SetPrec(53) + got.Sqrt(NewFloat(r)) + want := NewFloat(math.Sqrt(r)) + if got.Cmp(want) != 0 { + t.Fatalf("Sqrt(%g) =\n got %g;\nwant %g", r, got, want) + } + } +} + +func TestFloatSqrt(t *testing.T) { + for _, test := range []struct { + x string + want string + }{ + // Test values were generated on Wolfram Alpha using query + // 'sqrt(N) to 350 digits' + // 350 decimal digits give up to 1000 binary digits. + {"0.03125", "0.17677669529663688110021109052621225982120898442211850914708496724884155980776337985629844179095519659187673077886403712811560450698134215158051518713749197892665283324093819909447499381264409775757143376369499645074628431682460775184106467733011114982619404115381053858929018135497032545349940642599871090667456829147610370507757690729404938184321879"}, + {"0.125", "0.35355339059327376220042218105242451964241796884423701829416993449768311961552675971259688358191039318375346155772807425623120901396268430316103037427498395785330566648187639818894998762528819551514286752738999290149256863364921550368212935466022229965238808230762107717858036270994065090699881285199742181334913658295220741015515381458809876368643757"}, + {"0.5", "0.70710678118654752440084436210484903928483593768847403658833986899536623923105351942519376716382078636750692311545614851246241802792536860632206074854996791570661133296375279637789997525057639103028573505477998580298513726729843100736425870932044459930477616461524215435716072541988130181399762570399484362669827316590441482031030762917619752737287514"}, + {"2.0", "1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457503"}, + {"3.0", "1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756756261414154067030299699450949989524788116555120943736485280932319023055820679748201010846749232650153123432669033228866506722546689218379712270471316603678615880190499865373798593894676503475065760507566183481296061009476021871903250831458295239598"}, + {"4.0", "2.0"}, + + {"1p512", "1p256"}, + {"4p1024", "2p512"}, + {"9p2048", "3p1024"}, + + {"1p-1024", "1p-512"}, + {"4p-2048", "2p-1024"}, + {"9p-4096", "3p-2048"}, + } { + for _, prec := range []uint{24, 53, 64, 65, 100, 128, 129, 200, 256, 400, 600, 800, 1000} { + x := new(Float).SetPrec(prec) + x.Parse(test.x, 10) + + got := new(Float).SetPrec(prec).Sqrt(x) + want := new(Float).SetPrec(prec) + want.Parse(test.want, 10) + if got.Cmp(want) != 0 { + t.Errorf("prec = %d, Sqrt(%v) =\ngot %g;\nwant %g", + prec, test.x, got, want) + } + + // Square test. + // If got holds the square root of x to precision p, then + // got = √x + k + // for some k such that |k| < 2**(-p). Thus, + // got² = (√x + k)² = x + 2k√n + k² + // and the error must satisfy + // err = |got² - x| ≈ | 2k√n | < 2**(-p+1)*√n + // Ignoring the k² term for simplicity. + + // err = |got² - x| + // (but do intermediate steps with 32 guard digits to + // avoid introducing spurious rounding-related errors) + sq := new(Float).SetPrec(prec+32).Mul(got, got) + diff := new(Float).Sub(sq, x) + err := diff.Abs(diff).SetPrec(prec) + + // maxErr = 2**(-p+1)*√x + one := new(Float).SetPrec(prec).SetInt64(1) + maxErr := new(Float).Mul(new(Float).SetMantExp(one, -int(prec)+1), got) + + if err.Cmp(maxErr) >= 0 { + t.Errorf("prec = %d, Sqrt(%v) =\ngot err %g;\nwant maxErr %g", + prec, test.x, err, maxErr) + } + } + } +} + +func TestFloatSqrtSpecial(t *testing.T) { + for _, test := range []struct { + x *Float + want *Float + }{ + {NewFloat(+0), NewFloat(+0)}, + {NewFloat(-0), NewFloat(-0)}, + {NewFloat(math.Inf(+1)), NewFloat(math.Inf(+1))}, + } { + got := new(Float).Sqrt(test.x) + if got.neg != test.want.neg || got.form != test.want.form { + t.Errorf("Sqrt(%v) = %v (neg: %v); want %v (neg: %v)", + test.x, got, got.neg, test.want, test.want.neg) + } + } + +} + +// Benchmarks + +func BenchmarkFloatSqrt(b *testing.B) { + for _, prec := range []uint{64, 128, 256, 1e3, 1e4, 1e5, 1e6} { + x := NewFloat(2) + z := new(Float).SetPrec(prec) + b.Run(fmt.Sprintf("%v", prec), func(b *testing.B) { + b.ReportAllocs() + for n := 0; n < b.N; n++ { + z.Sqrt(x) + } + }) + } +} diff --git a/libgo/go/math/bits.go b/libgo/go/math/bits.go index d85ee9c..77bcdbe 100644 --- a/libgo/go/math/bits.go +++ b/libgo/go/math/bits.go @@ -8,9 +8,12 @@ const ( uvnan = 0x7FF8000000000001 uvinf = 0x7FF0000000000000 uvneginf = 0xFFF0000000000000 + uvone = 0x3FF0000000000000 mask = 0x7FF shift = 64 - 11 - 1 bias = 1023 + signMask = 1 << 63 + fracMask = 1<<shift - 1 ) // Inf returns positive infinity if sign >= 0, negative infinity if sign < 0. diff --git a/libgo/go/math/bits/example_test.go b/libgo/go/math/bits/example_test.go index 78750da..dfe7ece 100644 --- a/libgo/go/math/bits/example_test.go +++ b/libgo/go/math/bits/example_test.go @@ -2,6 +2,8 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. +// Code generated by go run make_examples.go. DO NOT EDIT. + // +build ignore package bits_test @@ -11,70 +13,194 @@ import ( "math/bits" ) +func ExampleLeadingZeros8() { + fmt.Printf("LeadingZeros8(%08b) = %d\n", 1, bits.LeadingZeros8(1)) + // Output: + // LeadingZeros8(00000001) = 7 +} + func ExampleLeadingZeros16() { - fmt.Println(bits.LeadingZeros16(0)) - fmt.Println(bits.LeadingZeros16(1)) - fmt.Println(bits.LeadingZeros16(256)) - fmt.Println(bits.LeadingZeros16(65535)) + fmt.Printf("LeadingZeros16(%016b) = %d\n", 1, bits.LeadingZeros16(1)) // Output: - // 16 - // 15 - // 7 - // 0 + // LeadingZeros16(0000000000000001) = 15 } func ExampleLeadingZeros32() { - fmt.Println(bits.LeadingZeros32(0)) - fmt.Println(bits.LeadingZeros32(1)) + fmt.Printf("LeadingZeros32(%032b) = %d\n", 1, bits.LeadingZeros32(1)) // Output: - // 32 - // 31 + // LeadingZeros32(00000000000000000000000000000001) = 31 } func ExampleLeadingZeros64() { - fmt.Println(bits.LeadingZeros64(0)) - fmt.Println(bits.LeadingZeros64(1)) + fmt.Printf("LeadingZeros64(%064b) = %d\n", 1, bits.LeadingZeros64(1)) // Output: - // 64 - // 63 + // LeadingZeros64(0000000000000000000000000000000000000000000000000000000000000001) = 63 } -func ExampleOnesCount() { - fmt.Printf("%b\n", 14) - fmt.Println(bits.OnesCount(14)) +func ExampleTrailingZeros8() { + fmt.Printf("TrailingZeros8(%08b) = %d\n", 14, bits.TrailingZeros8(14)) // Output: - // 1110 - // 3 + // TrailingZeros8(00001110) = 1 +} + +func ExampleTrailingZeros16() { + fmt.Printf("TrailingZeros16(%016b) = %d\n", 14, bits.TrailingZeros16(14)) + // Output: + // TrailingZeros16(0000000000001110) = 1 +} + +func ExampleTrailingZeros32() { + fmt.Printf("TrailingZeros32(%032b) = %d\n", 14, bits.TrailingZeros32(14)) + // Output: + // TrailingZeros32(00000000000000000000000000001110) = 1 +} + +func ExampleTrailingZeros64() { + fmt.Printf("TrailingZeros64(%064b) = %d\n", 14, bits.TrailingZeros64(14)) + // Output: + // TrailingZeros64(0000000000000000000000000000000000000000000000000000000000001110) = 1 } func ExampleOnesCount8() { - fmt.Printf("%b\n", 14) - fmt.Println(bits.OnesCount8(14)) + fmt.Printf("OnesCount8(%08b) = %d\n", 14, bits.OnesCount8(14)) // Output: - // 1110 - // 3 + // OnesCount8(00001110) = 3 } func ExampleOnesCount16() { - fmt.Printf("%b\n", 14) - fmt.Println(bits.OnesCount16(14)) + fmt.Printf("OnesCount16(%016b) = %d\n", 14, bits.OnesCount16(14)) // Output: - // 1110 - // 3 + // OnesCount16(0000000000001110) = 3 } func ExampleOnesCount32() { - fmt.Printf("%b\n", 14) - fmt.Println(bits.OnesCount32(14)) + fmt.Printf("OnesCount32(%032b) = %d\n", 14, bits.OnesCount32(14)) // Output: - // 1110 - // 3 + // OnesCount32(00000000000000000000000000001110) = 3 } func ExampleOnesCount64() { - fmt.Printf("%b\n", 14) - fmt.Println(bits.OnesCount64(14)) + fmt.Printf("OnesCount64(%064b) = %d\n", 14, bits.OnesCount64(14)) + // Output: + // OnesCount64(0000000000000000000000000000000000000000000000000000000000001110) = 3 +} + +func ExampleRotateLeft8() { + fmt.Printf("%08b\n", 15) + fmt.Printf("%08b\n", bits.RotateLeft8(15, 2)) + fmt.Printf("%08b\n", bits.RotateLeft8(15, -2)) + // Output: + // 00001111 + // 00111100 + // 11000011 +} + +func ExampleRotateLeft16() { + fmt.Printf("%016b\n", 15) + fmt.Printf("%016b\n", bits.RotateLeft16(15, 2)) + fmt.Printf("%016b\n", bits.RotateLeft16(15, -2)) + // Output: + // 0000000000001111 + // 0000000000111100 + // 1100000000000011 +} + +func ExampleRotateLeft32() { + fmt.Printf("%032b\n", 15) + fmt.Printf("%032b\n", bits.RotateLeft32(15, 2)) + fmt.Printf("%032b\n", bits.RotateLeft32(15, -2)) + // Output: + // 00000000000000000000000000001111 + // 00000000000000000000000000111100 + // 11000000000000000000000000000011 +} + +func ExampleRotateLeft64() { + fmt.Printf("%064b\n", 15) + fmt.Printf("%064b\n", bits.RotateLeft64(15, 2)) + fmt.Printf("%064b\n", bits.RotateLeft64(15, -2)) + // Output: + // 0000000000000000000000000000000000000000000000000000000000001111 + // 0000000000000000000000000000000000000000000000000000000000111100 + // 1100000000000000000000000000000000000000000000000000000000000011 +} + +func ExampleReverse8() { + fmt.Printf("%08b\n", 19) + fmt.Printf("%08b\n", bits.Reverse8(19)) + // Output: + // 00010011 + // 11001000 +} + +func ExampleReverse16() { + fmt.Printf("%016b\n", 19) + fmt.Printf("%016b\n", bits.Reverse16(19)) + // Output: + // 0000000000010011 + // 1100100000000000 +} + +func ExampleReverse32() { + fmt.Printf("%032b\n", 19) + fmt.Printf("%032b\n", bits.Reverse32(19)) + // Output: + // 00000000000000000000000000010011 + // 11001000000000000000000000000000 +} + +func ExampleReverse64() { + fmt.Printf("%064b\n", 19) + fmt.Printf("%064b\n", bits.Reverse64(19)) + // Output: + // 0000000000000000000000000000000000000000000000000000000000010011 + // 1100100000000000000000000000000000000000000000000000000000000000 +} + +func ExampleReverseBytes16() { + fmt.Printf("%016b\n", 15) + fmt.Printf("%016b\n", bits.ReverseBytes16(15)) + // Output: + // 0000000000001111 + // 0000111100000000 +} + +func ExampleReverseBytes32() { + fmt.Printf("%032b\n", 15) + fmt.Printf("%032b\n", bits.ReverseBytes32(15)) + // Output: + // 00000000000000000000000000001111 + // 00001111000000000000000000000000 +} + +func ExampleReverseBytes64() { + fmt.Printf("%064b\n", 15) + fmt.Printf("%064b\n", bits.ReverseBytes64(15)) + // Output: + // 0000000000000000000000000000000000000000000000000000000000001111 + // 0000111100000000000000000000000000000000000000000000000000000000 +} + +func ExampleLen8() { + fmt.Printf("Len8(%08b) = %d\n", 8, bits.Len8(8)) + // Output: + // Len8(00001000) = 4 +} + +func ExampleLen16() { + fmt.Printf("Len16(%016b) = %d\n", 8, bits.Len16(8)) + // Output: + // Len16(0000000000001000) = 4 +} + +func ExampleLen32() { + fmt.Printf("Len32(%032b) = %d\n", 8, bits.Len32(8)) + // Output: + // Len32(00000000000000000000000000001000) = 4 +} + +func ExampleLen64() { + fmt.Printf("Len64(%064b) = %d\n", 8, bits.Len64(8)) // Output: - // 1110 - // 3 + // Len64(0000000000000000000000000000000000000000000000000000000000001000) = 4 } diff --git a/libgo/go/math/bits/make_examples.go b/libgo/go/math/bits/make_examples.go new file mode 100644 index 0000000..cd81cd6 --- /dev/null +++ b/libgo/go/math/bits/make_examples.go @@ -0,0 +1,112 @@ +// Copyright 2017 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. + +// +build ignore + +// This program generates example_test.go. + +package main + +import ( + "bytes" + "fmt" + "io/ioutil" + "log" + "math/bits" +) + +const header = `// Copyright 2017 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. + +// Code generated by go run make_examples.go. DO NOT EDIT. + +package bits_test + +import ( + "fmt" + "math/bits" +) +` + +func main() { + w := bytes.NewBuffer([]byte(header)) + + for _, e := range []struct { + name string + in int + out [4]interface{} + out2 [4]interface{} + }{ + { + name: "LeadingZeros", + in: 1, + out: [4]interface{}{bits.LeadingZeros8(1), bits.LeadingZeros16(1), bits.LeadingZeros32(1), bits.LeadingZeros64(1)}, + }, + { + name: "TrailingZeros", + in: 14, + out: [4]interface{}{bits.TrailingZeros8(14), bits.TrailingZeros16(14), bits.TrailingZeros32(14), bits.TrailingZeros64(14)}, + }, + { + name: "OnesCount", + in: 14, + out: [4]interface{}{bits.OnesCount8(14), bits.OnesCount16(14), bits.OnesCount32(14), bits.OnesCount64(14)}, + }, + { + name: "RotateLeft", + in: 15, + out: [4]interface{}{bits.RotateLeft8(15, 2), bits.RotateLeft16(15, 2), bits.RotateLeft32(15, 2), bits.RotateLeft64(15, 2)}, + out2: [4]interface{}{bits.RotateLeft8(15, -2), bits.RotateLeft16(15, -2), bits.RotateLeft32(15, -2), bits.RotateLeft64(15, -2)}, + }, + { + name: "Reverse", + in: 19, + out: [4]interface{}{bits.Reverse8(19), bits.Reverse16(19), bits.Reverse32(19), bits.Reverse64(19)}, + }, + { + name: "ReverseBytes", + in: 15, + out: [4]interface{}{nil, bits.ReverseBytes16(15), bits.ReverseBytes32(15), bits.ReverseBytes64(15)}, + }, + { + name: "Len", + in: 8, + out: [4]interface{}{bits.Len8(8), bits.Len16(8), bits.Len32(8), bits.Len64(8)}, + }, + } { + for i, size := range []int{8, 16, 32, 64} { + if e.out[i] == nil { + continue // function doesn't exist + } + f := fmt.Sprintf("%s%d", e.name, size) + fmt.Fprintf(w, "\nfunc Example%s() {\n", f) + switch e.name { + case "RotateLeft", "Reverse", "ReverseBytes": + fmt.Fprintf(w, "\tfmt.Printf(\"%%0%db\\n\", %d)\n", size, e.in) + if e.name == "RotateLeft" { + fmt.Fprintf(w, "\tfmt.Printf(\"%%0%db\\n\", bits.%s(%d, 2))\n", size, f, e.in) + fmt.Fprintf(w, "\tfmt.Printf(\"%%0%db\\n\", bits.%s(%d, -2))\n", size, f, e.in) + } else { + fmt.Fprintf(w, "\tfmt.Printf(\"%%0%db\\n\", bits.%s(%d))\n", size, f, e.in) + } + fmt.Fprintf(w, "\t// Output:\n") + fmt.Fprintf(w, "\t// %0*b\n", size, e.in) + fmt.Fprintf(w, "\t// %0*b\n", size, e.out[i]) + if e.name == "RotateLeft" && e.out2[i] != nil { + fmt.Fprintf(w, "\t// %0*b\n", size, e.out2[i]) + } + default: + fmt.Fprintf(w, "\tfmt.Printf(\"%s(%%0%db) = %%d\\n\", %d, bits.%s(%d))\n", f, size, e.in, f, e.in) + fmt.Fprintf(w, "\t// Output:\n") + fmt.Fprintf(w, "\t// %s(%0*b) = %d\n", f, size, e.in, e.out[i]) + } + fmt.Fprintf(w, "}\n") + } + } + + if err := ioutil.WriteFile("example_test.go", w.Bytes(), 0666); err != nil { + log.Fatal(err) + } +} diff --git a/libgo/go/math/cmplx/asin.go b/libgo/go/math/cmplx/asin.go index 61880a2..062f324 100644 --- a/libgo/go/math/cmplx/asin.go +++ b/libgo/go/math/cmplx/asin.go @@ -49,11 +49,8 @@ import "math" // Asin returns the inverse sine of x. func Asin(x complex128) complex128 { - if imag(x) == 0 { - if math.Abs(real(x)) > 1 { - return complex(math.Pi/2, 0) // DOMAIN error - } - return complex(math.Asin(real(x)), 0) + if imag(x) == 0 && math.Abs(real(x)) <= 1 { + return complex(math.Asin(real(x)), imag(x)) } ct := complex(-imag(x), real(x)) // i * x xx := x * x @@ -65,12 +62,8 @@ func Asin(x complex128) complex128 { // Asinh returns the inverse hyperbolic sine of x. func Asinh(x complex128) complex128 { - // TODO check range - if imag(x) == 0 { - if math.Abs(real(x)) > 1 { - return complex(math.Pi/2, 0) // DOMAIN error - } - return complex(math.Asinh(real(x)), 0) + if imag(x) == 0 && math.Abs(real(x)) <= 1 { + return complex(math.Asinh(real(x)), imag(x)) } xx := x * x x1 := complex(1+real(xx), imag(xx)) // 1 + x*x @@ -140,10 +133,6 @@ func Acosh(x complex128) complex128 { // Atan returns the inverse tangent of x. func Atan(x complex128) complex128 { - if real(x) == 0 && imag(x) > 1 { - return NaN() - } - x2 := real(x) * real(x) a := 1 - x2 - imag(x)*imag(x) if a == 0 { diff --git a/libgo/go/math/cmplx/cmath_test.go b/libgo/go/math/cmplx/cmath_test.go index 7a5c485..8d70562 100644 --- a/libgo/go/math/cmplx/cmath_test.go +++ b/libgo/go/math/cmplx/cmath_test.go @@ -435,6 +435,24 @@ var tanhSC = []complex128{ NaN(), } +// branch cut continuity checks +// points on each axis at |z| > 1 are checked for one-sided continuity from both the positive and negative side +// all possible branch cuts for the elementary functions are at one of these points + +var zero = 0.0 +var eps = 1.0 / (1 << 53) + +var branchPoints = [][2]complex128{ + {complex(2.0, zero), complex(2.0, eps)}, + {complex(2.0, -zero), complex(2.0, -eps)}, + {complex(-2.0, zero), complex(-2.0, eps)}, + {complex(-2.0, -zero), complex(-2.0, -eps)}, + {complex(zero, 2.0), complex(eps, 2.0)}, + {complex(-zero, 2.0), complex(-eps, 2.0)}, + {complex(zero, -2.0), complex(eps, -2.0)}, + {complex(-zero, -2.0), complex(-eps, -2.0)}, +} + // functions borrowed from pkg/math/all_test.go func tolerance(a, b, e float64) bool { d := a - b @@ -508,6 +526,11 @@ func TestAcos(t *testing.T) { t.Errorf("Acos(%g) = %g, want %g", vcAcosSC[i], f, acosSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Acos(pt[0]), Acos(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Acos(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestAcosh(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -520,6 +543,11 @@ func TestAcosh(t *testing.T) { t.Errorf("Acosh(%g) = %g, want %g", vcAcoshSC[i], f, acoshSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Acosh(pt[0]), Acosh(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Acosh(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestAsin(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -532,6 +560,11 @@ func TestAsin(t *testing.T) { t.Errorf("Asin(%g) = %g, want %g", vcAsinSC[i], f, asinSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Asin(pt[0]), Asin(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Asin(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestAsinh(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -544,6 +577,11 @@ func TestAsinh(t *testing.T) { t.Errorf("Asinh(%g) = %g, want %g", vcAsinhSC[i], f, asinhSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Asinh(pt[0]), Asinh(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Asinh(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestAtan(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -556,6 +594,11 @@ func TestAtan(t *testing.T) { t.Errorf("Atan(%g) = %g, want %g", vcAtanSC[i], f, atanSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Atan(pt[0]), Atan(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Atan(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestAtanh(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -568,6 +611,11 @@ func TestAtanh(t *testing.T) { t.Errorf("Atanh(%g) = %g, want %g", vcAtanhSC[i], f, atanhSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Atanh(pt[0]), Atanh(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Atanh(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestConj(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -635,6 +683,11 @@ func TestLog(t *testing.T) { t.Errorf("Log(%g) = %g, want %g", vcLogSC[i], f, logSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Log(pt[0]), Log(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Log(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestLog10(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -685,6 +738,11 @@ func TestPow(t *testing.T) { t.Errorf("Pow(%g, %g) = %g, want %g", vcPowSC[i][0], vcPowSC[i][0], f, powSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Pow(pt[0], 0.1), Pow(pt[1], 0.1); !cVeryclose(f0, f1) { + t.Errorf("Pow(%g, 0.1) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestRect(t *testing.T) { for i := 0; i < len(vc); i++ { @@ -733,6 +791,11 @@ func TestSqrt(t *testing.T) { t.Errorf("Sqrt(%g) = %g, want %g", vcSqrtSC[i], f, sqrtSC[i]) } } + for _, pt := range branchPoints { + if f0, f1 := Sqrt(pt[0]), Sqrt(pt[1]); !cVeryclose(f0, f1) { + t.Errorf("Sqrt(%g) not continuous, got %g want %g", pt[0], f0, f1) + } + } } func TestTan(t *testing.T) { for i := 0; i < len(vc); i++ { diff --git a/libgo/go/math/cmplx/sqrt.go b/libgo/go/math/cmplx/sqrt.go index 72f81e9..0fbdcde 100644 --- a/libgo/go/math/cmplx/sqrt.go +++ b/libgo/go/math/cmplx/sqrt.go @@ -57,13 +57,14 @@ import "math" // The result r is chosen so that real(r) ≥ 0 and imag(r) has the same sign as imag(x). func Sqrt(x complex128) complex128 { if imag(x) == 0 { + // Ensure that imag(r) has the same sign as imag(x) for imag(x) == signed zero. if real(x) == 0 { - return complex(0, 0) + return complex(0, imag(x)) } if real(x) < 0 { - return complex(0, math.Sqrt(-real(x))) + return complex(0, math.Copysign(math.Sqrt(-real(x)), imag(x))) } - return complex(math.Sqrt(real(x)), 0) + return complex(math.Sqrt(real(x)), imag(x)) } if real(x) == 0 { if imag(x) < 0 { diff --git a/libgo/go/math/const.go b/libgo/go/math/const.go index 20b7065..0fc8715 100644 --- a/libgo/go/math/const.go +++ b/libgo/go/math/const.go @@ -9,18 +9,18 @@ package math // Mathematical constants. const ( - E = 2.71828182845904523536028747135266249775724709369995957496696763 // http://oeis.org/A001113 - Pi = 3.14159265358979323846264338327950288419716939937510582097494459 // http://oeis.org/A000796 - Phi = 1.61803398874989484820458683436563811772030917980576286213544862 // http://oeis.org/A001622 + E = 2.71828182845904523536028747135266249775724709369995957496696763 // https://oeis.org/A001113 + Pi = 3.14159265358979323846264338327950288419716939937510582097494459 // https://oeis.org/A000796 + Phi = 1.61803398874989484820458683436563811772030917980576286213544862 // https://oeis.org/A001622 - Sqrt2 = 1.41421356237309504880168872420969807856967187537694807317667974 // http://oeis.org/A002193 - SqrtE = 1.64872127070012814684865078781416357165377610071014801157507931 // http://oeis.org/A019774 - SqrtPi = 1.77245385090551602729816748334114518279754945612238712821380779 // http://oeis.org/A002161 - SqrtPhi = 1.27201964951406896425242246173749149171560804184009624861664038 // http://oeis.org/A139339 + Sqrt2 = 1.41421356237309504880168872420969807856967187537694807317667974 // https://oeis.org/A002193 + SqrtE = 1.64872127070012814684865078781416357165377610071014801157507931 // https://oeis.org/A019774 + SqrtPi = 1.77245385090551602729816748334114518279754945612238712821380779 // https://oeis.org/A002161 + SqrtPhi = 1.27201964951406896425242246173749149171560804184009624861664038 // https://oeis.org/A139339 - Ln2 = 0.693147180559945309417232121458176568075500134360255254120680009 // http://oeis.org/A002162 + Ln2 = 0.693147180559945309417232121458176568075500134360255254120680009 // https://oeis.org/A002162 Log2E = 1 / Ln2 - Ln10 = 2.30258509299404568401799145468436420760110148862877297603332790 // http://oeis.org/A002392 + Ln10 = 2.30258509299404568401799145468436420760110148862877297603332790 // https://oeis.org/A002392 Log10E = 1 / Ln10 ) diff --git a/libgo/go/math/dim.go b/libgo/go/math/dim.go index 37ab538..cf06f30 100644 --- a/libgo/go/math/dim.go +++ b/libgo/go/math/dim.go @@ -11,11 +11,18 @@ package math // Dim(-Inf, -Inf) = NaN // Dim(x, NaN) = Dim(NaN, x) = NaN func Dim(x, y float64) float64 { - return dim(x, y) -} - -func dim(x, y float64) float64 { - return max(x-y, 0) + // The special cases result in NaN after the subtraction: + // +Inf - +Inf = NaN + // -Inf - -Inf = NaN + // NaN - y = NaN + // x - NaN = NaN + v := x - y + if v <= 0 { + // v is negative or 0 + return 0 + } + // v is positive or NaN + return v } // Max returns the larger of x or y. diff --git a/libgo/go/math/erfinv.go b/libgo/go/math/erfinv.go new file mode 100644 index 0000000..21b5578 --- /dev/null +++ b/libgo/go/math/erfinv.go @@ -0,0 +1,127 @@ +// Copyright 2017 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 + +/* + Inverse of the floating-point error function. +*/ + +// This implementation is based on the rational approximation +// of percentage points of normal distribution available from +// http://www.jstor.org/stable/2347330. + +const ( + // Coefficients for approximation to erf in |x| <= 0.85 + a0 = 1.1975323115670912564578e0 + a1 = 4.7072688112383978012285e1 + a2 = 6.9706266534389598238465e2 + a3 = 4.8548868893843886794648e3 + a4 = 1.6235862515167575384252e4 + a5 = 2.3782041382114385731252e4 + a6 = 1.1819493347062294404278e4 + a7 = 8.8709406962545514830200e2 + b0 = 1.0000000000000000000e0 + b1 = 4.2313330701600911252e1 + b2 = 6.8718700749205790830e2 + b3 = 5.3941960214247511077e3 + b4 = 2.1213794301586595867e4 + b5 = 3.9307895800092710610e4 + b6 = 2.8729085735721942674e4 + b7 = 5.2264952788528545610e3 + // Coefficients for approximation to erf in 0.85 < |x| <= 1-2*exp(-25) + c0 = 1.42343711074968357734e0 + c1 = 4.63033784615654529590e0 + c2 = 5.76949722146069140550e0 + c3 = 3.64784832476320460504e0 + c4 = 1.27045825245236838258e0 + c5 = 2.41780725177450611770e-1 + c6 = 2.27238449892691845833e-2 + c7 = 7.74545014278341407640e-4 + d0 = 1.4142135623730950488016887e0 + d1 = 2.9036514445419946173133295e0 + d2 = 2.3707661626024532365971225e0 + d3 = 9.7547832001787427186894837e-1 + d4 = 2.0945065210512749128288442e-1 + d5 = 2.1494160384252876777097297e-2 + d6 = 7.7441459065157709165577218e-4 + d7 = 1.4859850019840355905497876e-9 + // Coefficients for approximation to erf in 1-2*exp(-25) < |x| < 1 + e0 = 6.65790464350110377720e0 + e1 = 5.46378491116411436990e0 + e2 = 1.78482653991729133580e0 + e3 = 2.96560571828504891230e-1 + e4 = 2.65321895265761230930e-2 + e5 = 1.24266094738807843860e-3 + e6 = 2.71155556874348757815e-5 + e7 = 2.01033439929228813265e-7 + f0 = 1.414213562373095048801689e0 + f1 = 8.482908416595164588112026e-1 + f2 = 1.936480946950659106176712e-1 + f3 = 2.103693768272068968719679e-2 + f4 = 1.112800997078859844711555e-3 + f5 = 2.611088405080593625138020e-5 + f6 = 2.010321207683943062279931e-7 + f7 = 2.891024605872965461538222e-15 +) + +// Erfinv returns the inverse error function of x. +// +// Special cases are: +// Erfinv(1) = +Inf +// Erfinv(-1) = -Inf +// Erfinv(x) = NaN if x < -1 or x > 1 +// Erfinv(NaN) = NaN +func Erfinv(x float64) float64 { + // special cases + if IsNaN(x) || x <= -1 || x >= 1 { + if x == -1 || x == 1 { + return Inf(int(x)) + } + return NaN() + } + + sign := false + if x < 0 { + x = -x + sign = true + } + + var ans float64 + if x <= 0.85 { // |x| <= 0.85 + r := 0.180625 - 0.25*x*x + z1 := ((((((a7*r+a6)*r+a5)*r+a4)*r+a3)*r+a2)*r+a1)*r + a0 + z2 := ((((((b7*r+b6)*r+b5)*r+b4)*r+b3)*r+b2)*r+b1)*r + b0 + ans = (x * z1) / z2 + } else { + var z1, z2 float64 + r := Sqrt(Ln2 - Log(1.0-x)) + if r <= 5.0 { + r -= 1.6 + z1 = ((((((c7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r + c0 + z2 = ((((((d7*r+d6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r + d0 + } else { + r -= 5.0 + z1 = ((((((e7*r+e6)*r+e5)*r+e4)*r+e3)*r+e2)*r+e1)*r + e0 + z2 = ((((((f7*r+f6)*r+f5)*r+f4)*r+f3)*r+f2)*r+f1)*r + f0 + } + ans = z1 / z2 + } + + if sign { + return -ans + } + return ans +} + +// Erfcinv returns the inverse of Erfc(x). +// +// Special cases are: +// Erfcinv(0) = +Inf +// Erfcinv(2) = -Inf +// Erfcinv(x) = NaN if x < 0 or x > 2 +// Erfcinv(NaN) = NaN +func Erfcinv(x float64) float64 { + return Erfinv(1 - x) +} diff --git a/libgo/go/math/example_test.go b/libgo/go/math/example_test.go index 12e9876..feaf9d8 100644 --- a/libgo/go/math/example_test.go +++ b/libgo/go/math/example_test.go @@ -9,6 +9,77 @@ import ( "math" ) +func ExampleAcos() { + fmt.Printf("%.2f", math.Acos(1)) + // Output: 0.00 +} + +func ExampleAcosh() { + fmt.Printf("%.2f", math.Acosh(1)) + // Output: 0.00 +} + +func ExampleAsin() { + fmt.Printf("%.2f", math.Asin(0)) + // Output: 0.00 +} + +func ExampleAsinh() { + fmt.Printf("%.2f", math.Asinh(0)) + // Output: 0.00 +} + +func ExampleAtan() { + fmt.Printf("%.2f", math.Atan(0)) + // Output: 0.00 +} + +func ExampleAtan2() { + fmt.Printf("%.2f", math.Atan2(0, 0)) + // Output: 0.00 +} + +func ExampleAtanh() { + fmt.Printf("%.2f", math.Atanh(0)) + // Output: 0.00 +} + +func ExampleCos() { + fmt.Printf("%.2f", math.Cos(math.Pi/2)) + // Output: 0.00 +} + +func ExampleCosh() { + fmt.Printf("%.2f", math.Cosh(0)) + // Output: 1.00 +} + +func ExampleSin() { + fmt.Printf("%.2f", math.Sin(math.Pi)) + // Output: 0.00 +} + +func ExampleSincos() { + sin, cos := math.Sincos(0) + fmt.Printf("%.2f, %.2f", sin, cos) + // Output: 0.00, 1.00 +} + +func ExampleSinh() { + fmt.Printf("%.2f", math.Sinh(0)) + // Output: 0.00 +} + +func ExampleTan() { + fmt.Printf("%.2f", math.Tan(0)) + // Output: 0.00 +} + +func ExampleTanh() { + fmt.Printf("%.2f", math.Tanh(0)) + // Output: 0.00 +} + func ExampleSqrt() { const ( a = 3 diff --git a/libgo/go/math/exp.go b/libgo/go/math/exp.go index cbc955d..f3d3cb6 100644 --- a/libgo/go/math/exp.go +++ b/libgo/go/math/exp.go @@ -50,7 +50,7 @@ func Exp(x float64) float64 { // the interval [0,0.34658]: // Write // R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... -// We use a special Remes algorithm on [0,0.34658] to generate +// We use a special Remez algorithm on [0,0.34658] to generate // a polynomial of degree 5 to approximate R. The maximum error // of this polynomial approximation is bounded by 2**-59. In // other words, @@ -183,7 +183,7 @@ func exp2(x float64) float64 { // exp1 returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2. func expmulti(hi, lo float64, k int) float64 { const ( - P1 = 1.66666666666666019037e-01 /* 0x3FC55555; 0x5555553E */ + P1 = 1.66666666666666657415e-01 /* 0x3FC55555; 0x55555555 */ P2 = -2.77777777770155933842e-03 /* 0xBF66C16C; 0x16BEBD93 */ P3 = 6.61375632143793436117e-05 /* 0x3F11566A; 0xAF25DE2C */ P4 = -1.65339022054652515390e-06 /* 0xBEBBBD41; 0xC5D26BF1 */ diff --git a/libgo/go/math/exp_asm.go b/libgo/go/math/exp_asm.go new file mode 100644 index 0000000..298420f --- /dev/null +++ b/libgo/go/math/exp_asm.go @@ -0,0 +1,12 @@ +// Copyright 2017 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. + +// +build amd64 amd64p32 +// +build ignore + +package math + +import "internal/cpu" + +var useFMA = cpu.X86.HasAVX && cpu.X86.HasFMA diff --git a/libgo/go/math/expm1.go b/libgo/go/math/expm1.go index 7494043..9c6e080 100644 --- a/libgo/go/math/expm1.go +++ b/libgo/go/math/expm1.go @@ -214,33 +214,33 @@ func expm1(x float64) float64 { r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))) t := 3 - r1*hfx e := hxs * ((r1 - t) / (6.0 - x*t)) - if k != 0 { - e = (x*(e-c) - c) - e -= hxs - switch { - case k == -1: - return 0.5*(x-e) - 0.5 - case k == 1: - if x < -0.25 { - return -2 * (e - (x + 0.5)) - } - return 1 + 2*(x-e) - case k <= -2 || k > 56: // suffice to return exp(x)-1 - y := 1 - (e - x) - y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent - return y - 1 - } - if k < 20 { - t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k - y := t - (e - x) - y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent - return y + if k == 0 { + return x - (x*e - hxs) // c is 0 + } + e = (x*(e-c) - c) + e -= hxs + switch { + case k == -1: + return 0.5*(x-e) - 0.5 + case k == 1: + if x < -0.25 { + return -2 * (e - (x + 0.5)) } - t := Float64frombits(uint64(0x3ff-k) << 52) // 2**-k - y := x - (e + t) - y++ + return 1 + 2*(x-e) + case k <= -2 || k > 56: // suffice to return exp(x)-1 + y := 1 - (e - x) + y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent + return y - 1 + } + if k < 20 { + t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k + y := t - (e - x) y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent return y } - return x - (x*e - hxs) // c is 0 + t = Float64frombits(uint64(0x3ff-k) << 52) // 2**-k + y := x - (e + t) + y++ + y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent + return y } diff --git a/libgo/go/math/floor.go b/libgo/go/math/floor.go index c40be41..9115968 100644 --- a/libgo/go/math/floor.go +++ b/libgo/go/math/floor.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. @@ -69,3 +69,78 @@ func trunc(x float64) float64 { d, _ := Modf(x) return d } + +// Round returns the nearest integer, rounding half away from zero. +// +// Special cases are: +// Round(±0) = ±0 +// Round(±Inf) = ±Inf +// Round(NaN) = NaN +func Round(x float64) float64 { + // Round is a faster implementation of: + // + // func Round(x float64) float64 { + // t := Trunc(x) + // if Abs(x-t) >= 0.5 { + // return t + Copysign(1, x) + // } + // return t + // } + bits := Float64bits(x) + e := uint(bits>>shift) & mask + if e < bias { + // Round abs(x) < 1 including denormals. + bits &= signMask // +-0 + if e == bias-1 { + bits |= uvone // +-1 + } + } else if e < bias+shift { + // Round any abs(x) >= 1 containing a fractional component [0,1). + // + // Numbers with larger exponents are returned unchanged since they + // must be either an integer, infinity, or NaN. + const half = 1 << (shift - 1) + e -= bias + bits += half >> e + bits &^= fracMask >> e + } + return Float64frombits(bits) +} + +// RoundToEven returns the nearest integer, rounding ties to even. +// +// Special cases are: +// RoundToEven(±0) = ±0 +// RoundToEven(±Inf) = ±Inf +// RoundToEven(NaN) = NaN +func RoundToEven(x float64) float64 { + // RoundToEven is a faster implementation of: + // + // func RoundToEven(x float64) float64 { + // t := math.Trunc(x) + // odd := math.Remainder(t, 2) != 0 + // if d := math.Abs(x - t); d > 0.5 || (d == 0.5 && odd) { + // return t + math.Copysign(1, x) + // } + // return t + // } + bits := Float64bits(x) + e := uint(bits>>shift) & mask + if e >= bias { + // Round abs(x) >= 1. + // - Large numbers without fractional components, infinity, and NaN are unchanged. + // - Add 0.499.. or 0.5 before truncating depending on whether the truncated + // number is even or odd (respectively). + const halfMinusULP = (1 << (shift - 1)) - 1 + e -= bias + bits += (halfMinusULP + (bits>>(shift-e))&1) >> e + bits &^= fracMask >> e + } else if e == bias-1 && bits&fracMask != 0 { + // Round 0.5 < abs(x) < 1. + bits = bits&signMask | uvone // +-1 + } else { + // Round abs(x) <= 0.5 including denormals. + bits &= signMask // +-0 + } + return Float64frombits(bits) +} diff --git a/libgo/go/math/pow.go b/libgo/go/math/pow.go index 7ba426f4..e255ee4 100644 --- a/libgo/go/math/pow.go +++ b/libgo/go/math/pow.go @@ -99,7 +99,16 @@ func pow(x, y float64) float64 { return NaN() } if yi >= 1<<63 { - return Exp(y * Log(x)) + // yi is a large even int that will lead to overflow (or underflow to 0) + // for all x except -1 (x == 1 was handled earlier) + switch { + case x == -1: + return 1 + case (Abs(x) < 1) == (y > 0): + return 0 + default: + return Inf(1) + } } // ans = a1 * 2**ae (= 1 for now). @@ -121,6 +130,15 @@ func pow(x, y float64) float64 { // accumulate powers of two into exp. x1, xe := Frexp(x) for i := int64(yi); i != 0; i >>= 1 { + if xe < -1<<12 || 1<<12 < xe { + // catch xe before it overflows the left shift below + // Since i !=0 it has at least one bit still set, so ae will accumulate xe + // on at least one more iteration, ae += xe is a lower bound on ae + // the lower bound on ae exceeds the size of a float64 exp + // so the final call to Ldexp will produce under/overflow (0/Inf) + ae += xe + break + } if i&1 == 1 { a1 *= x1 ae += xe diff --git a/libgo/go/math/rand/rand.go b/libgo/go/math/rand/rand.go index fe99c94..957bebd 100644 --- a/libgo/go/math/rand/rand.go +++ b/libgo/go/math/rand/rand.go @@ -135,6 +135,30 @@ func (r *Rand) Int31n(n int32) int32 { return v % n } +// int31n returns, as an int32, a non-negative pseudo-random number in [0,n). +// n must be > 0, but int31n does not check this; the caller must ensure it. +// int31n exists because Int31n is inefficient, but Go 1 compatibility +// requires that the stream of values produced by math/rand remain unchanged. +// int31n can thus only be used internally, by newly introduced APIs. +// +// For implementation details, see: +// https://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-modulo-reduction +// https://lemire.me/blog/2016/06/30/fast-random-shuffling +func (r *Rand) int31n(n int32) int32 { + v := r.Uint32() + prod := uint64(v) * uint64(n) + low := uint32(prod) + if low < uint32(n) { + thresh := uint32(-n) % uint32(n) + for low < thresh { + v = r.Uint32() + prod = uint64(v) * uint64(n) + low = uint32(prod) + } + } + return int32(prod >> 32) +} + // Intn returns, as an int, a non-negative pseudo-random number in [0,n). // It panics if n <= 0. func (r *Rand) Intn(n int) int { @@ -202,6 +226,31 @@ func (r *Rand) Perm(n int) []int { return m } +// Shuffle pseudo-randomizes the order of elements. +// n is the number of elements. Shuffle panics if n < 0. +// swap swaps the elements with indexes i and j. +func (r *Rand) Shuffle(n int, swap func(i, j int)) { + if n < 0 { + panic("invalid argument to Shuffle") + } + + // Fisher-Yates shuffle: https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle + // Shuffle really ought not be called with n that doesn't fit in 32 bits. + // Not only will it take a very long time, but with 2³¹! possible permutations, + // there's no way that any PRNG can have a big enough internal state to + // generate even a minuscule percentage of the possible permutations. + // Nevertheless, the right API signature accepts an int n, so handle it as best we can. + i := n - 1 + for ; i > 1<<31-1-1; i-- { + j := int(r.Int63n(int64(i + 1))) + swap(i, j) + } + for ; i > 0; i-- { + j := int(r.int31n(int32(i + 1))) + swap(i, j) + } +} + // Read generates len(p) random bytes and writes them into p. It // always returns len(p) and a nil error. // Read should not be called concurrently with any other Rand method. @@ -288,6 +337,11 @@ func Float32() float32 { return globalRand.Float32() } // from the default Source. func Perm(n int) []int { return globalRand.Perm(n) } +// Shuffle pseudo-randomizes the order of elements using the default Source. +// n is the number of elements. Shuffle panics if n < 0. +// swap swaps the elements with indexes i and j. +func Shuffle(n int, swap func(i, j int)) { globalRand.Shuffle(n, swap) } + // Read generates len(p) random bytes from the default Source and // writes them into p. It always returns len(p) and a nil error. // Read, unlike the Rand.Read method, is safe for concurrent use. diff --git a/libgo/go/math/rand/rand_test.go b/libgo/go/math/rand/rand_test.go index bf509e0..e663b84 100644 --- a/libgo/go/math/rand/rand_test.go +++ b/libgo/go/math/rand/rand_test.go @@ -53,7 +53,7 @@ func (this *statsResults) checkSimilarDistribution(expected *statsResults) error fmt.Println(s) return errors.New(s) } - if !nearEqual(this.stddev, expected.stddev, 0, expected.maxError) { + if !nearEqual(this.stddev, expected.stddev, expected.closeEnough, expected.maxError) { s := fmt.Sprintf("stddev %v != %v (allowed error %v, %v)", this.stddev, expected.stddev, expected.closeEnough, expected.maxError) fmt.Println(s) return errors.New(s) @@ -74,6 +74,7 @@ func getStatsResults(samples []float64) *statsResults { } func checkSampleDistribution(t *testing.T, samples []float64, expected *statsResults) { + t.Helper() actual := getStatsResults(samples) err := actual.checkSimilarDistribution(expected) if err != nil { @@ -82,6 +83,7 @@ func checkSampleDistribution(t *testing.T, samples []float64, expected *statsRes } func checkSampleSliceDistributions(t *testing.T, samples []float64, nslices int, expected *statsResults) { + t.Helper() chunk := len(samples) / nslices for i := 0; i < nslices; i++ { low := i * chunk @@ -374,7 +376,7 @@ func testReadUniformity(t *testing.T, n int, seed int64) { // Expect a uniform distribution of byte values, which lie in [0, 255]. var ( mean = 255.0 / 2 - stddev = math.Sqrt(255.0 * 255.0 / 12.0) + stddev = 256.0 / math.Sqrt(12.0) errorScale = stddev / math.Sqrt(float64(n)) ) @@ -448,6 +450,113 @@ func TestReadSeedReset(t *testing.T) { } } +func TestShuffleSmall(t *testing.T) { + // Check that Shuffle allows n=0 and n=1, but that swap is never called for them. + r := New(NewSource(1)) + for n := 0; n <= 1; n++ { + r.Shuffle(n, func(i, j int) { t.Fatalf("swap called, n=%d i=%d j=%d", n, i, j) }) + } +} + +// encodePerm converts from a permuted slice of length n, such as Perm generates, to an int in [0, n!). +// See https://en.wikipedia.org/wiki/Lehmer_code. +// encodePerm modifies the input slice. +func encodePerm(s []int) int { + // Convert to Lehmer code. + for i, x := range s { + r := s[i+1:] + for j, y := range r { + if y > x { + r[j]-- + } + } + } + // Convert to int in [0, n!). + m := 0 + fact := 1 + for i := len(s) - 1; i >= 0; i-- { + m += s[i] * fact + fact *= len(s) - i + } + return m +} + +// TestUniformFactorial tests several ways of generating a uniform value in [0, n!). +func TestUniformFactorial(t *testing.T) { + r := New(NewSource(testSeeds[0])) + top := 6 + if testing.Short() { + top = 4 + } + for n := 3; n <= top; n++ { + t.Run(fmt.Sprintf("n=%d", n), func(t *testing.T) { + // Calculate n!. + nfact := 1 + for i := 2; i <= n; i++ { + nfact *= i + } + + // Test a few different ways to generate a uniform distribution. + p := make([]int, n) // re-usable slice for Shuffle generator + tests := [...]struct { + name string + fn func() int + }{ + {name: "Int31n", fn: func() int { return int(r.Int31n(int32(nfact))) }}, + {name: "int31n", fn: func() int { return int(r.int31n(int32(nfact))) }}, + {name: "Perm", fn: func() int { return encodePerm(r.Perm(n)) }}, + {name: "Shuffle", fn: func() int { + // Generate permutation using Shuffle. + for i := range p { + p[i] = i + } + r.Shuffle(n, func(i, j int) { p[i], p[j] = p[j], p[i] }) + return encodePerm(p) + }}, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + // Gather chi-squared values and check that they follow + // the expected normal distribution given n!-1 degrees of freedom. + // See https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test and + // https://www.johndcook.com/Beautiful_Testing_ch10.pdf. + nsamples := 10 * nfact + if nsamples < 200 { + nsamples = 200 + } + samples := make([]float64, nsamples) + for i := range samples { + // Generate some uniformly distributed values and count their occurrences. + const iters = 1000 + counts := make([]int, nfact) + for i := 0; i < iters; i++ { + counts[test.fn()]++ + } + // Calculate chi-squared and add to samples. + want := iters / float64(nfact) + var χ2 float64 + for _, have := range counts { + err := float64(have) - want + χ2 += err * err + } + χ2 /= want + samples[i] = χ2 + } + + // Check that our samples approximate the appropriate normal distribution. + dof := float64(nfact - 1) + expected := &statsResults{mean: dof, stddev: math.Sqrt(2 * dof)} + errorScale := max(1.0, expected.stddev) + expected.closeEnough = 0.10 * errorScale + expected.maxError = 0.08 // TODO: What is the right value here? See issue 21211. + checkSampleDistribution(t, samples, expected) + }) + } + }) + } +} + // Benchmarks func BenchmarkInt63Threadsafe(b *testing.B) { @@ -512,6 +621,30 @@ func BenchmarkPerm30(b *testing.B) { } } +func BenchmarkPerm30ViaShuffle(b *testing.B) { + r := New(NewSource(1)) + for n := b.N; n > 0; n-- { + p := make([]int, 30) + for i := range p { + p[i] = i + } + r.Shuffle(30, func(i, j int) { p[i], p[j] = p[j], p[i] }) + } +} + +// BenchmarkShuffleOverhead uses a minimal swap function +// to measure just the shuffling overhead. +func BenchmarkShuffleOverhead(b *testing.B) { + r := New(NewSource(1)) + for n := b.N; n > 0; n-- { + r.Shuffle(52, func(i, j int) { + if i < 0 || i >= 52 || j < 0 || j >= 52 { + b.Fatalf("bad swap(%d, %d)", i, j) + } + }) + } +} + func BenchmarkRead3(b *testing.B) { r := New(NewSource(1)) buf := make([]byte, 3) |