diff --git a/src/math/acos_s390x.s b/src/math/acos_s390x.s
index 306f45a..d2288b8 100644
--- a/src/math/acos_s390x.s
+++ b/src/math/acos_s390x.s
@@ -42,7 +42,7 @@
 TEXT	·acosAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·acosrodataL13<>+0(SB), R9
-	WORD	$0xB3CD00C0	//lgdr %r12, %f0
+	LGDR	F0, R12
 	FMOVD	F0, F10
 	SRAD	$32, R12
 	WORD	$0xC0293FE6	//iilf	%r2,1072079005
diff --git a/src/math/acosh_s390x.s b/src/math/acosh_s390x.s
index 3575ed6..87a5d00 100644
--- a/src/math/acosh_s390x.s
+++ b/src/math/acosh_s390x.s
@@ -53,7 +53,7 @@
 TEXT	·acoshAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·acoshrodataL11<>+0(SB), R9
-	WORD	$0xB3CD0010	//lgdr %r1, %f0
+	LGDR	F0, R1
 	WORD	$0xC0295FEF	//iilf	%r2,1609564159
 	BYTE	$0xFF
 	BYTE	$0xFF
@@ -85,7 +85,7 @@
 	WORD	$0xC0398006	//iilf	%r3,2147909631
 	BYTE	$0x7F
 	BYTE	$0xFF
-	WORD	$0xB3CD0050	//lgdr %r5, %f0
+	LGDR	F0, R5
 	SRAD	$32, R5
 	MOVH	$0x0, R1
 	SUBW	R5, R3
@@ -105,7 +105,7 @@
 	SRAW	$8, R2, R2
 	ORW	$0x45000000, R2
 L5:
-	WORD	$0xB3C10001	//ldgr	%f0,%r1
+	LDGR	R1, F0
 	FMOVD	104(R9), F2
 	FMADD	F8, F0, F2
 	FMOVD	96(R9), F4
@@ -153,7 +153,7 @@
 	WORD	$0xC0398006	//iilf	%r3,2147909631
 	BYTE	$0x7F
 	BYTE	$0xFF
-	WORD	$0xB3CD0050	//lgdr %r5, %f0
+	LGDR	F0, R5
 	SRAD	$32, R5
 	MOVH	$0x0, R1
 	SUBW	R5, R3
diff --git a/src/math/all_test.go b/src/math/all_test.go
index bcc20a3..ed42941 100644
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -128,7 +128,7 @@
 var ceil = []float64{
 	5.0000000000000000e+00,
 	8.0000000000000000e+00,
-	0.0000000000000000e+00,
+	Copysign(0, -1),
 	-5.0000000000000000e+00,
 	1.0000000000000000e+01,
 	3.0000000000000000e+00,
@@ -175,6 +175,7 @@
 	-2.51772931436786954751e-01,
 	-7.3924135157173099849e-01,
 }
+
 var cosh = []float64{
 	7.2668796942212842775517446e+01,
 	1.1479413465659254502011135e+03,
@@ -644,7 +645,7 @@
 var trunc = []float64{
 	4.0000000000000000e+00,
 	7.0000000000000000e+00,
-	-0.0000000000000000e+00,
+	Copysign(0, -1),
 	-5.0000000000000000e+00,
 	9.0000000000000000e+00,
 	2.0000000000000000e+00,
@@ -1527,6 +1528,7 @@
 	0,
 	Inf(1),
 	NaN(),
+	4503599627370496.5, // Issue #29488
 }
 var log1pSC = []float64{
 	NaN(),
@@ -1536,6 +1538,7 @@
 	0,
 	Inf(1),
 	NaN(),
+	36.04365338911715, // Issue #29488
 }
 
 var vfmodfSC = []float64{
@@ -2158,7 +2161,7 @@
 
 func TestCeil(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if f := Ceil(vf[i]); ceil[i] != f {
+		if f := Ceil(vf[i]); !alike(ceil[i], f) {
 			t.Errorf("Ceil(%g) = %g, want %g", vf[i], f, ceil[i])
 		}
 	}
@@ -2385,7 +2388,7 @@
 
 func TestFloor(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if f := Floor(vf[i]); floor[i] != f {
+		if f := Floor(vf[i]); !alike(floor[i], f) {
 			t.Errorf("Floor(%g) = %g, want %g", vf[i], f, floor[i])
 		}
 	}
@@ -2916,7 +2919,7 @@
 
 func TestTrunc(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if f := Trunc(vf[i]); trunc[i] != f {
+		if f := Trunc(vf[i]); !alike(trunc[i], f) {
 			t.Errorf("Trunc(%g) = %g, want %g", vf[i], f, trunc[i])
 		}
 	}
@@ -3026,6 +3029,41 @@
 	}
 }
 
+// Check that trigReduce matches the standard reduction results for input values
+// below reduceThreshold.
+func TestTrigReduce(t *testing.T) {
+	inputs := make([]float64, len(vf))
+	// all of the standard inputs
+	copy(inputs, vf)
+	// all of the large inputs
+	large := float64(100000 * Pi)
+	for _, v := range vf {
+		inputs = append(inputs, v+large)
+	}
+	// Also test some special inputs, Pi and right below the reduceThreshold
+	inputs = append(inputs, Pi, Nextafter(ReduceThreshold, 0))
+	for _, x := range inputs {
+		// reduce the value to compare
+		j, z := TrigReduce(x)
+		xred := float64(j)*(Pi/4) + z
+
+		if f, fred := Sin(x), Sin(xred); !close(f, fred) {
+			t.Errorf("Sin(trigReduce(%g)) != Sin(%g), got %g, want %g", x, x, fred, f)
+		}
+		if f, fred := Cos(x), Cos(xred); !close(f, fred) {
+			t.Errorf("Cos(trigReduce(%g)) != Cos(%g), got %g, want %g", x, x, fred, f)
+		}
+		if f, fred := Tan(x), Tan(xred); !close(f, fred) {
+			t.Errorf(" Tan(trigReduce(%g)) != Tan(%g), got %g, want %g", x, x, fred, f)
+		}
+		f, g := Sincos(x)
+		fred, gred := Sincos(xred)
+		if !close(f, fred) || !close(g, gred) {
+			t.Errorf(" Sincos(trigReduce(%g)) != Sincos(%g), got %g, %g, want %g, %g", x, x, fred, gred, f, g)
+		}
+	}
+}
+
 // Check that math constants are accepted by compiler
 // and have right value (assumes strconv.ParseFloat works).
 // https://golang.org/issue/201
@@ -3635,3 +3673,41 @@
 	}
 	GlobalF = x
 }
+
+func BenchmarkFloat64bits(b *testing.B) {
+	y := uint64(0)
+	for i := 0; i < b.N; i++ {
+		y = Float64bits(roundNeg)
+	}
+	GlobalI = int(y)
+}
+
+var roundUint64 = uint64(5)
+
+func BenchmarkFloat64frombits(b *testing.B) {
+	x := 0.0
+	for i := 0; i < b.N; i++ {
+		x = Float64frombits(roundUint64)
+	}
+	GlobalF = x
+}
+
+var roundFloat32 = float32(-2.5)
+
+func BenchmarkFloat32bits(b *testing.B) {
+	y := uint32(0)
+	for i := 0; i < b.N; i++ {
+		y = Float32bits(roundFloat32)
+	}
+	GlobalI = int(y)
+}
+
+var roundUint32 = uint32(5)
+
+func BenchmarkFloat32frombits(b *testing.B) {
+	x := float32(0.0)
+	for i := 0; i < b.N; i++ {
+		x = Float32frombits(roundUint32)
+	}
+	GlobalF = float64(x)
+}
diff --git a/src/math/asin_s390x.s b/src/math/asin_s390x.s
index fd5ab04..dc54d05 100644
--- a/src/math/asin_s390x.s
+++ b/src/math/asin_s390x.s
@@ -46,7 +46,7 @@
 TEXT	·asinAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·asinrodataL15<>+0(SB), R9
-	WORD	$0xB3CD0070	//lgdr %r7, %f0
+	LGDR	F0, R7
 	FMOVD	F0, F8
 	SRAD	$32, R7
 	WORD	$0xC0193FE6 //iilf  %r1,1072079005
diff --git a/src/math/asinh_s390x.s b/src/math/asinh_s390x.s
index a9cee34..a3680c6 100644
--- a/src/math/asinh_s390x.s
+++ b/src/math/asinh_s390x.s
@@ -64,7 +64,7 @@
 TEXT	·asinhAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·asinhrodataL18<>+0(SB), R9
-	WORD	$0xB3CD00C0	//lgdr %r12, %f0
+	LGDR	F0, R12
 	WORD	$0xC0293FDF	//iilf	%r2,1071644671
 	BYTE	$0xFF
 	BYTE	$0xFF
@@ -93,7 +93,7 @@
 	WORD	$0xC0398006	//iilf	%r3,2147909631
 	BYTE	$0x7F
 	BYTE	$0xFF
-	WORD	$0xB3CD0050	//lgdr %r5, %f0
+	LGDR	F0, R5
 	SRAD	$32, R5
 	MOVH	$0x0, R2
 	SUBW	R5, R3
@@ -133,7 +133,7 @@
 	WORD	$0xC0398006	//iilf	%r3,2147909631
 	BYTE	$0x7F
 	BYTE	$0xFF
-	WORD	$0xB3CD0050	//lgdr %r5, %f0
+	LGDR	F0, R5
 	SRAD	$32, R5
 	MOVH	$0x0, R2
 	SUBW	R5, R3
@@ -146,7 +146,7 @@
 	BYTE	$0x59
 	ORW	$0x45000000, R1
 L6:
-	WORD	$0xB3C10022	//ldgr	%f2,%r2
+	LDGR	R2, F2
 	FMOVD	184(R9), F0
 	WFMADB	V8, V2, V0, V8
 	FMOVD	176(R9), F4
diff --git a/src/math/atan2_s390x.s b/src/math/atan2_s390x.s
index f37555b..c7a8a09 100644
--- a/src/math/atan2_s390x.s
+++ b/src/math/atan2_s390x.s
@@ -142,8 +142,8 @@
 	FMOVD	x+0(FP), F0
 	FMOVD	y+8(FP), F2
 	MOVD	$·atan2rodataL25<>+0(SB), R9
-	WORD	$0xB3CD0020	//lgdr	%r2,%f0
-	WORD	$0xB3CD0012	//lgdr	%r1,%f2
+	LGDR	F0, R2
+	LGDR	F2, R1
 	WORD	$0xEC2220BF	//risbgn	%r2,%r2,64-32,128+63,64+0+32
 	BYTE	$0x60
 	BYTE	$0x59
@@ -229,7 +229,7 @@
 	BYTE	$0x55
 	MOVD	$·atan2xpi2h<>+0(SB), R1
 	MOVD	·atan2xpim<>+0(SB), R3
-	WORD	$0xB3C10003	//ldgr	%f0,%r3
+	LDGR	R3, F0
 	WORD	$0xED021000	//madb	%f4,%f0,0(%r2,%r1)
 	BYTE	$0x40
 	BYTE	$0x1E
diff --git a/src/math/atan_s390x.s b/src/math/atan_s390x.s
index 9f4eaa2..713727d 100644
--- a/src/math/atan_s390x.s
+++ b/src/math/atan_s390x.s
@@ -54,7 +54,7 @@
 
 	MOVD	$·atanrodataL8<>+0(SB), R5
 	MOVH	$0x3FE0, R3
-	WORD	$0xB3CD0010	//lgdr	%r1,%f0
+	LGDR	F0, R1
 	WORD	$0xEC1120BF	//risbgn	%r1,%r1,64-32,128+63,64+0+32
 	BYTE	$0x60
 	BYTE	$0x59
diff --git a/src/math/atanh_s390x.s b/src/math/atanh_s390x.s
index 57b61a3..e7c6359 100644
--- a/src/math/atanh_s390x.s
+++ b/src/math/atanh_s390x.s
@@ -64,7 +64,7 @@
 TEXT    ·atanhAsm(SB), NOSPLIT, $0-16
     FMOVD   x+0(FP), F0
     MOVD    $·atanhrodataL10<>+0(SB), R5
-    WORD    $0xB3CD0010	//lgdr %r1, %f0
+    LGDR    F0, R1
     WORD    $0xC0393FEF //iilf  %r3,1072693247
     BYTE    $0xFF
     BYTE    $0xFF
@@ -128,7 +128,7 @@
     WORD    $0xED405088 //adb   %f4,.L12-.L10(%r5)
     BYTE    $0x00
     BYTE    $0x1A
-    WORD    $0xB3CD0044	//lgdr %r4, %f4
+    LGDR    F4, R4
     SRAD    $32, R4
     FMOVD   F4, F3
     WORD    $0xED305088 //sdb   %f3,.L12-.L10(%r5)
@@ -140,7 +140,7 @@
     BYTE    $0x00
     BYTE    $0x55
     SLD $32, R1, R1
-    WORD    $0xB3C10021 //ldgr  %f2,%r1
+    LDGR    R1, F2
     WFMADB  V4, V2, V16, V4
     SRAW    $8, R2, R1
     WFMADB  V4, V5, V6, V5
diff --git a/src/math/big/arith.go b/src/math/big/arith.go
index ad35240..f9db911 100644
--- a/src/math/big/arith.go
+++ b/src/math/big/arith.go
@@ -82,7 +82,7 @@
 	return uint(bits.LeadingZeros(uint(x)))
 }
 
-// q = (u1<<_W + u0 - r)/y
+// q = (u1<<_W + u0 - r)/v
 // Adapted from Warren, Hacker's Delight, p. 152.
 func divWW_g(u1, u0, v Word) (q, r Word) {
 	if u1 >= v {
diff --git a/src/math/big/arith_386.s b/src/math/big/arith_386.s
index 6c080f0..864fbc5 100644
--- a/src/math/big/arith_386.s
+++ b/src/math/big/arith_386.s
@@ -183,7 +183,7 @@
 	SHRL CX, DX:AX		// w>>s | w1<<ŝ
 	MOVL DX, (DI)(BX*4)	// z[i] = w>>s | w1<<ŝ
 	ADDL $1, BX		// i++
-	
+
 E9:	CMPL BX, BP
 	JL L9			// i < n-1
 
diff --git a/src/math/big/arith_amd64.s b/src/math/big/arith_amd64.s
index 1b950a4..e9c8887 100644
--- a/src/math/big/arith_amd64.s
+++ b/src/math/big/arith_amd64.s
@@ -324,10 +324,10 @@
 	MOVQ r+56(FP), CX	// c = r
 	MOVQ z_len+8(FP), R11
 	MOVQ $0, BX		// i = 0
-	
+
 	CMPQ R11, $4
 	JL E5
-	
+
 U5:	// i+4 <= n
 	// regular loop body unrolled 4x
 	MOVQ (0*8)(R8)(BX*8), AX
@@ -355,7 +355,7 @@
 	MOVQ AX, (3*8)(R10)(BX*8)
 	MOVQ DX, CX
 	ADDQ $4, BX		// i += 4
-	
+
 	LEAQ 4(BX), DX
 	CMPQ DX, R11
 	JLE U5
diff --git a/src/math/big/arith_arm.s b/src/math/big/arith_arm.s
index ba65fd2..33aa36f 100644
--- a/src/math/big/arith_arm.s
+++ b/src/math/big/arith_arm.s
@@ -123,7 +123,7 @@
 	MOVW	z_len+4(FP), R5
 	TEQ	$0, R5
 	BEQ	X7
-	
+
 	MOVW	z+0(FP), R1
 	MOVW	x+12(FP), R2
 	ADD	R5<<2, R2, R2
@@ -135,7 +135,7 @@
 	MOVW	$32, R4
 	SUB	R3, R4
 	MOVW	$0, R7
-	
+
 	MOVW.W	-4(R2), R6
 	MOVW	R6<<R3, R7
 	MOVW	R6>>R4, R6
diff --git a/src/math/big/arith_s390x.s b/src/math/big/arith_s390x.s
index 4520d16..9156d9d 100644
--- a/src/math/big/arith_s390x.s
+++ b/src/math/big/arith_s390x.s
@@ -54,7 +54,7 @@
 TEXT ·addVV(SB),NOSPLIT,$0
 	MOVD	addvectorfacility+0x00(SB),R1
 	BR	(R1)
-	
+
 TEXT ·addVV_check(SB),NOSPLIT, $0
 	MOVB	·hasVX(SB), R1
 	CMPBEQ	R1, $1, vectorimpl      // vectorfacility = 1, vector supported
@@ -89,7 +89,7 @@
 	BLT	v1
 	SUB     $12, R3                 // n -= 16
         BLT     A1                      // if n < 0 goto A1
-       
+
 	MOVD	R8, R5
 	MOVD	R9, R6
 	MOVD	R2, R7
@@ -291,7 +291,7 @@
 TEXT ·subVV(SB),NOSPLIT,$0
 	MOVD	subvectorfacility+0x00(SB),R1
 	BR	(R1)
-	
+
 TEXT ·subVV_check(SB),NOSPLIT,$0
 	MOVB	·hasVX(SB), R1
 	CMPBEQ	R1, $1, vectorimpl      // vectorfacility = 1, vector supported
@@ -321,7 +321,7 @@
 	MOVD	$0, R4		// c = 0
 	MOVD	$0, R0		// make sure it's zero
 	MOVD	$0, R10		// i = 0
-	
+
 	// s/JL/JMP/ below to disable the unrolled loop
 	SUB	$4, R3		// n -= 4
 	BLT	v1		// if n < 0 goto v1
@@ -413,7 +413,7 @@
 
 A1:	ADD	$12, R3		// n += 16
 	BLT	v1		// if n < 0 goto v1
-	
+
 U1:	// n >= 0
 	// regular loop body unrolled 4x
 	MOVD	0(R8)(R10*1), R5
@@ -532,7 +532,7 @@
 TEXT ·addVW(SB),NOSPLIT,$0
 	MOVD	addwvectorfacility+0x00(SB),R1
 	BR	(R1)
-	
+
 TEXT ·addVW_check(SB),NOSPLIT,$0
 	MOVB	·hasVX(SB), R1
 	CMPBEQ	R1, $1, vectorimpl      // vectorfacility = 1, vector supported
@@ -742,7 +742,7 @@
 TEXT ·subVW(SB),NOSPLIT,$0
 	MOVD	subwvectorfacility+0x00(SB),R1
 	BR	(R1)
-	
+
 TEXT ·subVW_check(SB),NOSPLIT,$0
 	MOVB	·hasVX(SB), R1
 	CMPBEQ	R1, $1, vectorimpl      // vectorfacility = 1, vector supported
diff --git a/src/math/big/float.go b/src/math/big/float.go
index 55b93c8..b3c3295 100644
--- a/src/math/big/float.go
+++ b/src/math/big/float.go
@@ -43,7 +43,7 @@
 // precision of the argument with the largest precision value before any
 // rounding takes place, and the rounding mode remains unchanged. Thus,
 // uninitialized Floats provided as result arguments will have their
-// precision set to a reasonable value determined by the operands and
+// precision set to a reasonable value determined by the operands, and
 // their mode is the zero value for RoundingMode (ToNearestEven).
 //
 // By setting the desired precision to 24 or 53 and using matching rounding
@@ -56,6 +56,12 @@
 // The zero (uninitialized) value for a Float is ready to use and represents
 // the number +0.0 exactly, with precision 0 and rounding mode ToNearestEven.
 //
+// Operations always take pointer arguments (*Float) rather
+// than Float values, and each unique Float value requires
+// its own unique *Float pointer. To "copy" a Float value,
+// an existing (or newly allocated) Float must be set to
+// a new value using the Float.Set method; shallow copies
+// of Floats are not supported and may lead to errors.
 type Float struct {
 	prec uint32
 	mode RoundingMode
@@ -293,7 +299,7 @@
 	z.round(sbit)
 }
 
-// SetMantExp sets z to mant × 2**exp and and returns z.
+// SetMantExp sets z to mant × 2**exp and returns z.
 // The result z has the same precision and rounding mode
 // as mant. SetMantExp is an inverse of MantExp but does
 // not require 0.5 <= |mant| < 1.0. Specifically:
@@ -321,7 +327,7 @@
 	return z
 }
 
-// Signbit returns true if x is negative or negative zero.
+// Signbit reports whether x is negative or negative zero.
 func (x *Float) Signbit() bool {
 	return x.neg
 }
diff --git a/src/math/big/int.go b/src/math/big/int.go
index 47a288a..dab9a5c 100644
--- a/src/math/big/int.go
+++ b/src/math/big/int.go
@@ -15,6 +15,13 @@
 
 // An Int represents a signed multi-precision integer.
 // The zero value for an Int represents the value 0.
+//
+// Operations always take pointer arguments (*Int) rather
+// than Int values, and each unique Int value requires
+// its own unique *Int pointer. To "copy" an Int value,
+// an existing (or newly allocated) Int must be set to
+// a new value using the Int.Set method; shallow copies
+// of Ints are not supported and may lead to errors.
 type Int struct {
 	neg bool // sign
 	abs nat  // absolute value of the integer
diff --git a/src/math/big/int_test.go b/src/math/big/int_test.go
index 9930ed0..7ef2b39 100644
--- a/src/math/big/int_test.go
+++ b/src/math/big/int_test.go
@@ -1727,3 +1727,29 @@
 		})
 	}
 }
+
+func benchmarkDiv(b *testing.B, aSize, bSize int) {
+	var r = rand.New(rand.NewSource(1234))
+	aa := randInt(r, uint(aSize))
+	bb := randInt(r, uint(bSize))
+	if aa.Cmp(bb) < 0 {
+		aa, bb = bb, aa
+	}
+	x := new(Int)
+	y := new(Int)
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		x.DivMod(aa, bb, y)
+	}
+}
+
+func BenchmarkDiv(b *testing.B) {
+	min, max, step := 10, 100000, 10
+	for i := min; i <= max; i *= step {
+		j := 2 * i
+		b.Run(fmt.Sprintf("%d/%d", j, i), func(b *testing.B) {
+			benchmarkDiv(b, j, i)
+		})
+	}
+}
diff --git a/src/math/big/nat.go b/src/math/big/nat.go
index a6f79ed..1e4a3b0 100644
--- a/src/math/big/nat.go
+++ b/src/math/big/nat.go
@@ -58,6 +58,10 @@
 	if n <= cap(z) {
 		return z[:n] // reuse z
 	}
+	if n == 1 {
+		// Most nats start small and stay that way; don't over-allocate.
+		return make(nat, 1)
+	}
 	// Choosing a good value for e has significant performance impact
 	// because it increases the chance that a value can be reused.
 	const e = 4 // extra capacity
@@ -680,43 +684,36 @@
 
 var natPool sync.Pool
 
-// q = (uIn-r)/v, with 0 <= r < y
+// q = (uIn-r)/vIn, with 0 <= r < y
 // Uses z as storage for q, and u as storage for r if possible.
 // See Knuth, Volume 2, section 4.3.1, Algorithm D.
 // Preconditions:
-//    len(v) >= 2
-//    len(uIn) >= len(v)
-func (z nat) divLarge(u, uIn, v nat) (q, r nat) {
-	n := len(v)
+//    len(vIn) >= 2
+//    len(uIn) >= len(vIn)
+//    u must not alias z
+func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) {
+	n := len(vIn)
 	m := len(uIn) - n
 
-	// 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, u) || alias(z, uIn) || alias(z, v) {
-		z = nil // z is an alias for u or uIn or v - cannot reuse
+	// D1.
+	shift := nlz(vIn[n-1])
+	// do not modify vIn, it may be used by another goroutine simultaneously
+	vp := getNat(n)
+	v := *vp
+	shlVU(v, vIn, shift)
+
+	// u may safely alias uIn or vIn, the value of uIn is used to set u and vIn was already used
+	u = u.make(len(uIn) + 1)
+	u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift)
+
+	// z may safely alias uIn or vIn, both values were used already
+	if alias(z, u) {
+		z = nil // z is an alias for u - cannot reuse
 	}
 	q = z.make(m + 1)
 
 	qhatvp := getNat(n + 1)
 	qhatv := *qhatvp
-	if alias(u, uIn) || alias(u, v) {
-		u = nil // u is an alias for uIn or v - cannot reuse
-	}
-	u = u.make(len(uIn) + 1)
-	u.clear() // TODO(gri) no need to clear if we allocated a new u
-
-	// D1.
-	var v1p *nat
-	shift := nlz(v[n-1])
-	if shift > 0 {
-		// do not modify v, it may be used by another goroutine simultaneously
-		v1p = getNat(n)
-		v1 := *v1p
-		shlVU(v1, v, shift)
-		v = v1
-	}
-	u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift)
 
 	// D2.
 	vn1 := v[n-1]
@@ -756,9 +753,8 @@
 
 		q[j] = qhat
 	}
-	if v1p != nil {
-		putNat(v1p)
-	}
+
+	putNat(vp)
 	putNat(qhatvp)
 
 	q = q.norm()
diff --git a/src/math/big/prime.go b/src/math/big/prime.go
index 4c2c152..d9a5f1e 100644
--- a/src/math/big/prime.go
+++ b/src/math/big/prime.go
@@ -51,7 +51,7 @@
 	}
 
 	if w&1 == 0 {
-		return false // n is even
+		return false // x is even
 	}
 
 	const primesA = 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 37
diff --git a/src/math/big/rat.go b/src/math/big/rat.go
index 46d58fc..5d0800c 100644
--- a/src/math/big/rat.go
+++ b/src/math/big/rat.go
@@ -13,6 +13,13 @@
 
 // A Rat represents a quotient a/b of arbitrary precision.
 // The zero value for a Rat represents the value 0.
+//
+// Operations always take pointer arguments (*Rat) rather
+// than Rat values, and each unique Rat value requires
+// its own unique *Rat pointer. To "copy" a Rat value,
+// an existing (or newly allocated) Rat must be set to
+// a new value using the Rat.Set method; shallow copies
+// of Rats are not supported and may lead to errors.
 type Rat struct {
 	// To make zero values for Rat work w/o initialization,
 	// a zero value of b (len(b) == 0) acts like b == 1.
diff --git a/src/math/big/ratconv.go b/src/math/big/ratconv.go
index 157d8d0..5656280 100644
--- a/src/math/big/ratconv.go
+++ b/src/math/big/ratconv.go
@@ -38,8 +38,8 @@
 }
 
 // SetString sets z to the value of s and returns z and a boolean indicating
-// success. s can be given as a fraction "a/b" or as a floating-point number
-// optionally followed by an exponent. The entire string (not just a prefix)
+// success. s can be given as a fraction "a/b" or as a decimal floating-point
+// number optionally followed by an exponent. The entire string (not just a prefix)
 // must be valid for success. If the operation failed, the value of z is
 // undefined but the returned value is nil.
 func (z *Rat) SetString(s string) (*Rat, bool) {
@@ -78,6 +78,7 @@
 	}
 
 	// mantissa
+	// TODO(gri) allow other bases besides 10 for mantissa and exponent? (issue #29799)
 	var ecorr int
 	z.a.abs, _, ecorr, err = z.a.abs.scan(r, 10, true)
 	if err != nil {
diff --git a/src/math/big/sqrt.go b/src/math/big/sqrt.go
index b989649..53403aa 100644
--- a/src/math/big/sqrt.go
+++ b/src/math/big/sqrt.go
@@ -7,8 +7,6 @@
 import "math"
 
 var (
-	half  = NewFloat(0.5)
-	two   = NewFloat(2.0)
 	three = NewFloat(3.0)
 )
 
@@ -57,9 +55,9 @@
 	case 0:
 		// nothing to do
 	case 1:
-		z.Mul(two, z)
+		z.exp++
 	case -1:
-		z.Mul(half, z)
+		z.exp--
 	}
 	// 0.25 <= z < 2.0
 
@@ -96,7 +94,7 @@
 		u.prec = t.prec
 		u.Mul(t, t)        // u = t²
 		u.Add(u, x)        //   = t² + x
-		u.Mul(half, u)     //   = ½(t² + x)
+		u.exp--            //   = ½(t² + x)
 		return t.Quo(u, t) //   = ½(t² + x)/t
 	}
 
@@ -133,11 +131,13 @@
 	ng := func(t *Float) *Float {
 		u.prec = t.prec
 		v.prec = t.prec
-		u.Mul(t, t)           // u = t²
-		u.Mul(x, u)           //   = xt²
-		v.Sub(three, u)       // v = 3 - xt²
-		u.Mul(t, v)           // u = t(3 - xt²)
-		return t.Mul(half, u) //   = ½t(3 - xt²)
+		u.Mul(t, t)     // u = t²
+		u.Mul(x, u)     //   = xt²
+		v.Sub(three, u) // v = 3 - xt²
+		u.Mul(t, v)     // u = t(3 - xt²)
+		u.exp--         //   = ½t(3 - xt²)
+		return t.Set(u)
+
 	}
 
 	xf, _ := x.Float64()
diff --git a/src/math/bits/bits.go b/src/math/bits/bits.go
index 989baac..b06c363 100644
--- a/src/math/bits/bits.go
+++ b/src/math/bits/bits.go
@@ -8,6 +8,8 @@
 // functions for the predeclared unsigned integer types.
 package bits
 
+import _ "unsafe" // for go:linkname
+
 const uintSize = 32 << (^uint(0) >> 32 & 1) // 32 or 64
 
 // UintSize is the size of a uint in bits.
@@ -63,7 +65,7 @@
 }
 
 // TrailingZeros16 returns the number of trailing zero bits in x; the result is 16 for x == 0.
-func TrailingZeros16(x uint16) (n int) {
+func TrailingZeros16(x uint16) int {
 	if x == 0 {
 		return 16
 	}
@@ -328,3 +330,206 @@
 	}
 	return n + int(len8tab[x])
 }
+
+// --- Add with carry ---
+
+// Add returns the sum with carry of x, y and carry: sum = x + y + carry.
+// The carry input must be 0 or 1; otherwise the behavior is undefined.
+// The carryOut output is guaranteed to be 0 or 1.
+func Add(x, y, carry uint) (sum, carryOut uint) {
+	yc := y + carry
+	sum = x + yc
+	if sum < x || yc < y {
+		carryOut = 1
+	}
+	return
+}
+
+// Add32 returns the sum with carry of x, y and carry: sum = x + y + carry.
+// The carry input must be 0 or 1; otherwise the behavior is undefined.
+// The carryOut output is guaranteed to be 0 or 1.
+func Add32(x, y, carry uint32) (sum, carryOut uint32) {
+	yc := y + carry
+	sum = x + yc
+	if sum < x || yc < y {
+		carryOut = 1
+	}
+	return
+}
+
+// Add64 returns the sum with carry of x, y and carry: sum = x + y + carry.
+// The carry input must be 0 or 1; otherwise the behavior is undefined.
+// The carryOut output is guaranteed to be 0 or 1.
+func Add64(x, y, carry uint64) (sum, carryOut uint64) {
+	yc := y + carry
+	sum = x + yc
+	if sum < x || yc < y {
+		carryOut = 1
+	}
+	return
+}
+
+// --- Subtract with borrow ---
+
+// Sub returns the difference of x, y and borrow: diff = x - y - borrow.
+// The borrow input must be 0 or 1; otherwise the behavior is undefined.
+// The borrowOut output is guaranteed to be 0 or 1.
+func Sub(x, y, borrow uint) (diff, borrowOut uint) {
+	yb := y + borrow
+	diff = x - yb
+	if diff > x || yb < y {
+		borrowOut = 1
+	}
+	return
+}
+
+// Sub32 returns the difference of x, y and borrow, diff = x - y - borrow.
+// The borrow input must be 0 or 1; otherwise the behavior is undefined.
+// The borrowOut output is guaranteed to be 0 or 1.
+func Sub32(x, y, borrow uint32) (diff, borrowOut uint32) {
+	yb := y + borrow
+	diff = x - yb
+	if diff > x || yb < y {
+		borrowOut = 1
+	}
+	return
+}
+
+// Sub64 returns the difference of x, y and borrow: diff = x - y - borrow.
+// The borrow input must be 0 or 1; otherwise the behavior is undefined.
+// The borrowOut output is guaranteed to be 0 or 1.
+func Sub64(x, y, borrow uint64) (diff, borrowOut uint64) {
+	yb := y + borrow
+	diff = x - yb
+	if diff > x || yb < y {
+		borrowOut = 1
+	}
+	return
+}
+
+// --- Full-width multiply ---
+
+// Mul returns the full-width product of x and y: (hi, lo) = x * y
+// with the product bits' upper half returned in hi and the lower
+// half returned in lo.
+func Mul(x, y uint) (hi, lo uint) {
+	if UintSize == 32 {
+		h, l := Mul32(uint32(x), uint32(y))
+		return uint(h), uint(l)
+	}
+	h, l := Mul64(uint64(x), uint64(y))
+	return uint(h), uint(l)
+}
+
+// Mul32 returns the 64-bit product of x and y: (hi, lo) = x * y
+// with the product bits' upper half returned in hi and the lower
+// half returned in lo.
+func Mul32(x, y uint32) (hi, lo uint32) {
+	tmp := uint64(x) * uint64(y)
+	hi, lo = uint32(tmp>>32), uint32(tmp)
+	return
+}
+
+// Mul64 returns the 128-bit product of x and y: (hi, lo) = x * y
+// with the product bits' upper half returned in hi and the lower
+// half returned in lo.
+func Mul64(x, y uint64) (hi, lo uint64) {
+	const mask32 = 1<<32 - 1
+	x0 := x & mask32
+	x1 := x >> 32
+	y0 := y & mask32
+	y1 := y >> 32
+	w0 := x0 * y0
+	t := x1*y0 + w0>>32
+	w1 := t & mask32
+	w2 := t >> 32
+	w1 += x0 * y1
+	hi = x1*y1 + w2 + w1>>32
+	lo = x * y
+	return
+}
+
+// --- Full-width divide ---
+
+// Div returns the quotient and remainder of (hi, lo) divided by y:
+// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper
+// half in parameter hi and the lower half in parameter lo.
+// Div panics for y == 0 (division by zero) or y <= hi (quotient overflow).
+func Div(hi, lo, y uint) (quo, rem uint) {
+	if UintSize == 32 {
+		q, r := Div32(uint32(hi), uint32(lo), uint32(y))
+		return uint(q), uint(r)
+	}
+	q, r := Div64(uint64(hi), uint64(lo), uint64(y))
+	return uint(q), uint(r)
+}
+
+// Div32 returns the quotient and remainder of (hi, lo) divided by y:
+// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper
+// half in parameter hi and the lower half in parameter lo.
+// Div32 panics for y == 0 (division by zero) or y <= hi (quotient overflow).
+func Div32(hi, lo, y uint32) (quo, rem uint32) {
+	if y != 0 && y <= hi {
+		panic(overflowError)
+	}
+	z := uint64(hi)<<32 | uint64(lo)
+	quo, rem = uint32(z/uint64(y)), uint32(z%uint64(y))
+	return
+}
+
+// Div64 returns the quotient and remainder of (hi, lo) divided by y:
+// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper
+// half in parameter hi and the lower half in parameter lo.
+// Div64 panics for y == 0 (division by zero) or y <= hi (quotient overflow).
+func Div64(hi, lo, y uint64) (quo, rem uint64) {
+	const (
+		two32  = 1 << 32
+		mask32 = two32 - 1
+	)
+	if y == 0 {
+		panic(divideError)
+	}
+	if y <= hi {
+		panic(overflowError)
+	}
+
+	s := uint(LeadingZeros64(y))
+	y <<= s
+
+	yn1 := y >> 32
+	yn0 := y & mask32
+	un32 := hi<<s | lo>>(64-s)
+	un10 := lo << s
+	un1 := un10 >> 32
+	un0 := un10 & mask32
+	q1 := un32 / yn1
+	rhat := un32 - q1*yn1
+
+	for q1 >= two32 || q1*yn0 > two32*rhat+un1 {
+		q1--
+		rhat += yn1
+		if rhat >= two32 {
+			break
+		}
+	}
+
+	un21 := un32*two32 + un1 - q1*y
+	q0 := un21 / yn1
+	rhat = un21 - q0*yn1
+
+	for q0 >= two32 || q0*yn0 > two32*rhat+un0 {
+		q0--
+		rhat += yn1
+		if rhat >= two32 {
+			break
+		}
+	}
+
+	return q1*two32 + q0, (un21*two32 + un0 - q0*y) >> s
+}
+
+//go:linkname overflowError runtime.overflowError
+var overflowError error
+
+//go:linkname divideError runtime.divideError
+var divideError error
diff --git a/src/math/bits/bits_test.go b/src/math/bits/bits_test.go
index 5c34f6d..1ec5107 100644
--- a/src/math/bits/bits_test.go
+++ b/src/math/bits/bits_test.go
@@ -6,6 +6,7 @@
 
 import (
 	. "math/bits"
+	"runtime"
 	"testing"
 	"unsafe"
 )
@@ -705,6 +706,385 @@
 	}
 }
 
+const (
+	_M   = 1<<UintSize - 1
+	_M32 = 1<<32 - 1
+	_M64 = 1<<64 - 1
+)
+
+func TestAddSubUint(t *testing.T) {
+	test := func(msg string, f func(x, y, c uint) (z, cout uint), x, y, c, z, cout uint) {
+		z1, cout1 := f(x, y, c)
+		if z1 != z || cout1 != cout {
+			t.Errorf("%s: got z:cout = %#x:%#x; want %#x:%#x", msg, z1, cout1, z, cout)
+		}
+	}
+	for _, a := range []struct{ x, y, c, z, cout uint }{
+		{0, 0, 0, 0, 0},
+		{0, 1, 0, 1, 0},
+		{0, 0, 1, 1, 0},
+		{0, 1, 1, 2, 0},
+		{12345, 67890, 0, 80235, 0},
+		{12345, 67890, 1, 80236, 0},
+		{_M, 1, 0, 0, 1},
+		{_M, 0, 1, 0, 1},
+		{_M, 1, 1, 1, 1},
+		{_M, _M, 0, _M - 1, 1},
+		{_M, _M, 1, _M, 1},
+	} {
+		test("Add", Add, a.x, a.y, a.c, a.z, a.cout)
+		test("Add symmetric", Add, a.y, a.x, a.c, a.z, a.cout)
+		test("Sub", Sub, a.z, a.x, a.c, a.y, a.cout)
+		test("Sub symmetric", Sub, a.z, a.y, a.c, a.x, a.cout)
+	}
+}
+
+func TestAddSubUint32(t *testing.T) {
+	test := func(msg string, f func(x, y, c uint32) (z, cout uint32), x, y, c, z, cout uint32) {
+		z1, cout1 := f(x, y, c)
+		if z1 != z || cout1 != cout {
+			t.Errorf("%s: got z:cout = %#x:%#x; want %#x:%#x", msg, z1, cout1, z, cout)
+		}
+	}
+	for _, a := range []struct{ x, y, c, z, cout uint32 }{
+		{0, 0, 0, 0, 0},
+		{0, 1, 0, 1, 0},
+		{0, 0, 1, 1, 0},
+		{0, 1, 1, 2, 0},
+		{12345, 67890, 0, 80235, 0},
+		{12345, 67890, 1, 80236, 0},
+		{_M32, 1, 0, 0, 1},
+		{_M32, 0, 1, 0, 1},
+		{_M32, 1, 1, 1, 1},
+		{_M32, _M32, 0, _M32 - 1, 1},
+		{_M32, _M32, 1, _M32, 1},
+	} {
+		test("Add32", Add32, a.x, a.y, a.c, a.z, a.cout)
+		test("Add32 symmetric", Add32, a.y, a.x, a.c, a.z, a.cout)
+		test("Sub32", Sub32, a.z, a.x, a.c, a.y, a.cout)
+		test("Sub32 symmetric", Sub32, a.z, a.y, a.c, a.x, a.cout)
+	}
+}
+
+func TestAddSubUint64(t *testing.T) {
+	test := func(msg string, f func(x, y, c uint64) (z, cout uint64), x, y, c, z, cout uint64) {
+		z1, cout1 := f(x, y, c)
+		if z1 != z || cout1 != cout {
+			t.Errorf("%s: got z:cout = %#x:%#x; want %#x:%#x", msg, z1, cout1, z, cout)
+		}
+	}
+	for _, a := range []struct{ x, y, c, z, cout uint64 }{
+		{0, 0, 0, 0, 0},
+		{0, 1, 0, 1, 0},
+		{0, 0, 1, 1, 0},
+		{0, 1, 1, 2, 0},
+		{12345, 67890, 0, 80235, 0},
+		{12345, 67890, 1, 80236, 0},
+		{_M64, 1, 0, 0, 1},
+		{_M64, 0, 1, 0, 1},
+		{_M64, 1, 1, 1, 1},
+		{_M64, _M64, 0, _M64 - 1, 1},
+		{_M64, _M64, 1, _M64, 1},
+	} {
+		test("Add64", Add64, a.x, a.y, a.c, a.z, a.cout)
+		test("Add64 symmetric", Add64, a.y, a.x, a.c, a.z, a.cout)
+		test("Sub64", Sub64, a.z, a.x, a.c, a.y, a.cout)
+		test("Sub64 symmetric", Sub64, a.z, a.y, a.c, a.x, a.cout)
+	}
+}
+
+func TestMulDiv(t *testing.T) {
+	testMul := func(msg string, f func(x, y uint) (hi, lo uint), x, y, hi, lo uint) {
+		hi1, lo1 := f(x, y)
+		if hi1 != hi || lo1 != lo {
+			t.Errorf("%s: got hi:lo = %#x:%#x; want %#x:%#x", msg, hi1, lo1, hi, lo)
+		}
+	}
+	testDiv := func(msg string, f func(hi, lo, y uint) (q, r uint), hi, lo, y, q, r uint) {
+		q1, r1 := f(hi, lo, y)
+		if q1 != q || r1 != r {
+			t.Errorf("%s: got q:r = %#x:%#x; want %#x:%#x", msg, q1, r1, q, r)
+		}
+	}
+	for _, a := range []struct {
+		x, y      uint
+		hi, lo, r uint
+	}{
+		{1 << (UintSize - 1), 2, 1, 0, 1},
+		{_M, _M, _M - 1, 1, 42},
+	} {
+		testMul("Mul", Mul, a.x, a.y, a.hi, a.lo)
+		testMul("Mul symmetric", Mul, a.y, a.x, a.hi, a.lo)
+		testDiv("Div", Div, a.hi, a.lo+a.r, a.y, a.x, a.r)
+		testDiv("Div symmetric", Div, a.hi, a.lo+a.r, a.x, a.y, a.r)
+	}
+}
+
+func TestMulDiv32(t *testing.T) {
+	testMul := func(msg string, f func(x, y uint32) (hi, lo uint32), x, y, hi, lo uint32) {
+		hi1, lo1 := f(x, y)
+		if hi1 != hi || lo1 != lo {
+			t.Errorf("%s: got hi:lo = %#x:%#x; want %#x:%#x", msg, hi1, lo1, hi, lo)
+		}
+	}
+	testDiv := func(msg string, f func(hi, lo, y uint32) (q, r uint32), hi, lo, y, q, r uint32) {
+		q1, r1 := f(hi, lo, y)
+		if q1 != q || r1 != r {
+			t.Errorf("%s: got q:r = %#x:%#x; want %#x:%#x", msg, q1, r1, q, r)
+		}
+	}
+	for _, a := range []struct {
+		x, y      uint32
+		hi, lo, r uint32
+	}{
+		{1 << 31, 2, 1, 0, 1},
+		{0xc47dfa8c, 50911, 0x98a4, 0x998587f4, 13},
+		{_M32, _M32, _M32 - 1, 1, 42},
+	} {
+		testMul("Mul32", Mul32, a.x, a.y, a.hi, a.lo)
+		testMul("Mul32 symmetric", Mul32, a.y, a.x, a.hi, a.lo)
+		testDiv("Div32", Div32, a.hi, a.lo+a.r, a.y, a.x, a.r)
+		testDiv("Div32 symmetric", Div32, a.hi, a.lo+a.r, a.x, a.y, a.r)
+	}
+}
+
+func TestMulDiv64(t *testing.T) {
+	testMul := func(msg string, f func(x, y uint64) (hi, lo uint64), x, y, hi, lo uint64) {
+		hi1, lo1 := f(x, y)
+		if hi1 != hi || lo1 != lo {
+			t.Errorf("%s: got hi:lo = %#x:%#x; want %#x:%#x", msg, hi1, lo1, hi, lo)
+		}
+	}
+	testDiv := func(msg string, f func(hi, lo, y uint64) (q, r uint64), hi, lo, y, q, r uint64) {
+		q1, r1 := f(hi, lo, y)
+		if q1 != q || r1 != r {
+			t.Errorf("%s: got q:r = %#x:%#x; want %#x:%#x", msg, q1, r1, q, r)
+		}
+	}
+	for _, a := range []struct {
+		x, y      uint64
+		hi, lo, r uint64
+	}{
+		{1 << 63, 2, 1, 0, 1},
+		{0x3626229738a3b9, 0xd8988a9f1cc4a61, 0x2dd0712657fe8, 0x9dd6a3364c358319, 13},
+		{_M64, _M64, _M64 - 1, 1, 42},
+	} {
+		testMul("Mul64", Mul64, a.x, a.y, a.hi, a.lo)
+		testMul("Mul64 symmetric", Mul64, a.y, a.x, a.hi, a.lo)
+		testDiv("Div64", Div64, a.hi, a.lo+a.r, a.y, a.x, a.r)
+		testDiv("Div64 symmetric", Div64, a.hi, a.lo+a.r, a.x, a.y, a.r)
+	}
+}
+
+const (
+	divZeroError  = "runtime error: integer divide by zero"
+	overflowError = "runtime error: integer overflow"
+)
+
+func TestDivPanicOverflow(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div should have panicked when y<=hi")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != overflowError {
+			t.Errorf("Div expected panic: %q, got: %q ", overflowError, e.Error())
+		}
+	}()
+	q, r := Div(1, 0, 1)
+	t.Errorf("undefined q, r = %v, %v calculated when Div should have panicked", q, r)
+}
+
+func TestDiv32PanicOverflow(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div32 should have panicked when y<=hi")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != overflowError {
+			t.Errorf("Div32 expected panic: %q, got: %q ", overflowError, e.Error())
+		}
+	}()
+	q, r := Div32(1, 0, 1)
+	t.Errorf("undefined q, r = %v, %v calculated when Div32 should have panicked", q, r)
+}
+
+func TestDiv64PanicOverflow(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div64 should have panicked when y<=hi")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != overflowError {
+			t.Errorf("Div64 expected panic: %q, got: %q ", overflowError, e.Error())
+		}
+	}()
+	q, r := Div64(1, 0, 1)
+	t.Errorf("undefined q, r = %v, %v calculated when Div64 should have panicked", q, r)
+}
+
+func TestDivPanicZero(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div should have panicked when y==0")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != divZeroError {
+			t.Errorf("Div expected panic: %q, got: %q ", divZeroError, e.Error())
+		}
+	}()
+	q, r := Div(1, 1, 0)
+	t.Errorf("undefined q, r = %v, %v calculated when Div should have panicked", q, r)
+}
+
+func TestDiv32PanicZero(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div32 should have panicked when y==0")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != divZeroError {
+			t.Errorf("Div32 expected panic: %q, got: %q ", divZeroError, e.Error())
+		}
+	}()
+	q, r := Div32(1, 1, 0)
+	t.Errorf("undefined q, r = %v, %v calculated when Div32 should have panicked", q, r)
+}
+
+func TestDiv64PanicZero(t *testing.T) {
+	// Expect a panic
+	defer func() {
+		if err := recover(); err == nil {
+			t.Error("Div64 should have panicked when y==0")
+		} else if e, ok := err.(runtime.Error); !ok || e.Error() != divZeroError {
+			t.Errorf("Div64 expected panic: %q, got: %q ", divZeroError, e.Error())
+		}
+	}()
+	q, r := Div64(1, 1, 0)
+	t.Errorf("undefined q, r = %v, %v calculated when Div64 should have panicked", q, r)
+}
+
+func BenchmarkAdd(b *testing.B) {
+	var z, c uint
+	for i := 0; i < b.N; i++ {
+		z, c = Add(uint(Input), uint(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkAdd32(b *testing.B) {
+	var z, c uint32
+	for i := 0; i < b.N; i++ {
+		z, c = Add32(uint32(Input), uint32(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkAdd64(b *testing.B) {
+	var z, c uint64
+	for i := 0; i < b.N; i++ {
+		z, c = Add64(uint64(Input), uint64(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkAdd64multiple(b *testing.B) {
+	var z0 = uint64(Input)
+	var z1 = uint64(Input)
+	var z2 = uint64(Input)
+	var z3 = uint64(Input)
+	for i := 0; i < b.N; i++ {
+		var c uint64
+		z0, c = Add64(z0, uint64(i), c)
+		z1, c = Add64(z1, uint64(i), c)
+		z2, c = Add64(z2, uint64(i), c)
+		z3, _ = Add64(z3, uint64(i), c)
+	}
+	Output = int(z0 + z1 + z2 + z3)
+}
+
+func BenchmarkSub(b *testing.B) {
+	var z, c uint
+	for i := 0; i < b.N; i++ {
+		z, c = Sub(uint(Input), uint(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkSub32(b *testing.B) {
+	var z, c uint32
+	for i := 0; i < b.N; i++ {
+		z, c = Sub32(uint32(Input), uint32(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkSub64(b *testing.B) {
+	var z, c uint64
+	for i := 0; i < b.N; i++ {
+		z, c = Sub64(uint64(Input), uint64(i), c)
+	}
+	Output = int(z + c)
+}
+
+func BenchmarkSub64multiple(b *testing.B) {
+	var z0 = uint64(Input)
+	var z1 = uint64(Input)
+	var z2 = uint64(Input)
+	var z3 = uint64(Input)
+	for i := 0; i < b.N; i++ {
+		var c uint64
+		z0, c = Sub64(z0, uint64(i), c)
+		z1, c = Sub64(z1, uint64(i), c)
+		z2, c = Sub64(z2, uint64(i), c)
+		z3, _ = Sub64(z3, uint64(i), c)
+	}
+	Output = int(z0 + z1 + z2 + z3)
+}
+
+func BenchmarkMul(b *testing.B) {
+	var hi, lo uint
+	for i := 0; i < b.N; i++ {
+		hi, lo = Mul(uint(Input), uint(i))
+	}
+	Output = int(hi + lo)
+}
+
+func BenchmarkMul32(b *testing.B) {
+	var hi, lo uint32
+	for i := 0; i < b.N; i++ {
+		hi, lo = Mul32(uint32(Input), uint32(i))
+	}
+	Output = int(hi + lo)
+}
+
+func BenchmarkMul64(b *testing.B) {
+	var hi, lo uint64
+	for i := 0; i < b.N; i++ {
+		hi, lo = Mul64(uint64(Input), uint64(i))
+	}
+	Output = int(hi + lo)
+}
+
+func BenchmarkDiv(b *testing.B) {
+	var q, r uint
+	for i := 0; i < b.N; i++ {
+		q, r = Div(1, uint(i), uint(Input))
+	}
+	Output = int(q + r)
+}
+
+func BenchmarkDiv32(b *testing.B) {
+	var q, r uint32
+	for i := 0; i < b.N; i++ {
+		q, r = Div32(1, uint32(i), uint32(Input))
+	}
+	Output = int(q + r)
+}
+
+func BenchmarkDiv64(b *testing.B) {
+	var q, r uint64
+	for i := 0; i < b.N; i++ {
+		q, r = Div64(1, uint64(i), uint64(Input))
+	}
+	Output = int(q + r)
+}
+
 // ----------------------------------------------------------------------------
 // Testing support
 
diff --git a/src/math/cbrt_s390x.s b/src/math/cbrt_s390x.s
index 85a2fcb..d79b48f 100644
--- a/src/math/cbrt_s390x.s
+++ b/src/math/cbrt_s390x.s
@@ -77,7 +77,7 @@
 TEXT	·cbrtAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·cbrtrodataL9<>+0(SB), R9
-	WORD	$0xB3CD0020	//lgdr %r2, %f0
+	LGDR	F0, R2
 	WORD	$0xC039000F	//iilf	%r3,1048575
 	BYTE	$0xFF
 	BYTE	$0xFF
@@ -103,7 +103,7 @@
 	BYTE	$0x00
 	BYTE	$0x1C
 	MOVH	$0x200, R4
-	WORD	$0xB3CD0022	//lgdr %r2, %f2
+	LGDR	F2, R2
 	SRAD	$32, R2
 L4:
 	WORD	$0xEC3239BE	//risbg	%r3,%r2,57,128+62,64-25
@@ -134,7 +134,7 @@
 	ADDW	R4, R1
 	SLW	$16, R1, R1
 	SLD	$32, R1, R1
-	WORD	$0xB3C10021	//ldgr	%f2,%r1
+	LDGR	R1, F2
 	WFMDB	V2, V2, V4
 	WFMDB	V4, V0, V6
 	WFMSDB	V4, V6, V2, V4
diff --git a/src/math/cmplx/isinf.go b/src/math/cmplx/isinf.go
index d5a65b4..6273cd3 100644
--- a/src/math/cmplx/isinf.go
+++ b/src/math/cmplx/isinf.go
@@ -6,7 +6,7 @@
 
 import "math"
 
-// IsInf returns true if either real(x) or imag(x) is an infinity.
+// IsInf reports whether either real(x) or imag(x) is an infinity.
 func IsInf(x complex128) bool {
 	if math.IsInf(real(x), 0) || math.IsInf(imag(x), 0) {
 		return true
diff --git a/src/math/cmplx/isnan.go b/src/math/cmplx/isnan.go
index 05d0cce..d3382c0 100644
--- a/src/math/cmplx/isnan.go
+++ b/src/math/cmplx/isnan.go
@@ -6,7 +6,7 @@
 
 import "math"
 
-// IsNaN returns true if either real(x) or imag(x) is NaN
+// IsNaN reports whether either real(x) or imag(x) is NaN
 // and neither is an infinity.
 func IsNaN(x complex128) bool {
 	switch {
diff --git a/src/math/erf_s390x.s b/src/math/erf_s390x.s
index 5b62bda..5be5d4d 100644
--- a/src/math/erf_s390x.s
+++ b/src/math/erf_s390x.s
@@ -100,7 +100,7 @@
 TEXT	·erfAsm(SB), NOSPLIT, $0-16
 	FMOVD	x+0(FP), F0
 	MOVD	$·erfrodataL13<>+0(SB), R5
-	WORD	$0xB3CD0010	//lgdr %r1, %f0
+	LGDR	F0, R1
 	FMOVD	F0, F6
 	SRAD	$48, R1
 	MOVH	$16383, R3
@@ -205,7 +205,7 @@
 	FMOVD	256(R5), F4
 	WFMADB	V1, V4, V3, V4
 	FDIV	F6, F2
-	WORD	$0xB3CD0014	//lgdr %r1, %f4
+	LGDR	F4, R1
 	FSUB	F3, F4
 	FMOVD	248(R5), F6
 	WFMSDB	V4, V6, V1, V4
@@ -230,7 +230,7 @@
 	BYTE	$0x59
 	MOVD	$·erftab2066<>+0(SB), R1
 	FMOVD	192(R5), F1
-	WORD	$0xB3C10033	//ldgr %f3,%r3
+	LDGR	R3, F3
 	WORD	$0xED221000	//madb %f2,%f2,0(%r2,%r1)
 	BYTE	$0x20
 	BYTE	$0x1E
diff --git a/src/math/erfc_s390x.s b/src/math/erfc_s390x.s
index 57710b2..0cb606d 100644
--- a/src/math/erfc_s390x.s
+++ b/src/math/erfc_s390x.s
@@ -219,7 +219,7 @@
 	WFMADB	V0, V5, V3, V5
 	WFMADB	V6, V7, V2, V7
 L11:
-	WORD	$0xB3CD0065	//lgdr %r6, %f5
+	LGDR	F5, R6
 	WFSDB	V0, V0, V2
 	WORD	$0xED509298	//sdb	%f5,.L55-.L38(%r9)
 	BYTE	$0x00
@@ -253,7 +253,7 @@
 	BYTE	$0x30
 	BYTE	$0x59
 	WFMADB	V4, V0, V2, V4
-	WORD	$0xB3C10024	//ldgr	%f2,%r4
+	LDGR	R4, F2
 	FMADD	F4, F2, F2
 	MOVW	R2, R6
 	CMPBLE	R6, $0, L20
@@ -504,7 +504,7 @@
 	CMPBGT	R6, R7, L24
 
 	WORD	$0xA5400010	//iihh	%r4,16
-	WORD	$0xB3C10024	//ldgr	%f2,%r4
+	LDGR	R4, F2
 	FMUL	F2, F2
 	BR	L1
 L23:
@@ -521,7 +521,7 @@
 	CMPBGT	R6, R7, L25
 	WORD	$0xA5408010	//iihh	%r4,32784
 	FMOVD	568(R9), F2
-	WORD	$0xB3C10004	//ldgr	%f0,%r4
+	LDGR	R4, F0
 	FMADD	F2, F0, F2
 	BR	L1
 L25:
diff --git a/src/math/example_test.go b/src/math/example_test.go
index a1f764b..25d6975 100644
--- a/src/math/example_test.go
+++ b/src/math/example_test.go
@@ -113,3 +113,25 @@
 	fmt.Printf("%.1f", c)
 	// Output: 100.0
 }
+
+func ExampleRound() {
+	p := math.Round(10.5)
+	fmt.Printf("%.1f\n", p)
+
+	n := math.Round(-10.5)
+	fmt.Printf("%.1f\n", n)
+	// Output:
+	// 11.0
+	// -11.0
+}
+
+func ExampleRoundToEven() {
+	u := math.RoundToEven(11.5)
+	fmt.Printf("%.1f\n", u)
+
+	d := math.RoundToEven(12.5)
+	fmt.Printf("%.1f\n", d)
+	// Output:
+	// 12.0
+	// 12.0
+}
diff --git a/src/math/exp_s390x.s b/src/math/exp_s390x.s
index 613ec24..cef1ce7 100644
--- a/src/math/exp_s390x.s
+++ b/src/math/exp_s390x.s
@@ -84,7 +84,7 @@
 	FMOVD	32(R5), F4
 	FMUL	F0, F0
 	WFMADB	V2, V4, V1, V4
-	WORD	$0xB3CD0016	//lgdr	%r1,%f6
+	LGDR	F6, R1
 	FMOVD	24(R5), F1
 	WFMADB	V2, V3, V1, V3
 	FMOVD	16(R5), F1
@@ -100,7 +100,7 @@
 	FMADD	F4, F2, F2
 	SLD	$48, R1, R2
 	WFMADB	V2, V0, V4, V2
-	WORD	$0xB3C10002	//ldgr	%f0,%r2
+	LDGR	R2, F0
 	FMADD	F0, F2, F0
 	FMOVD	F0, ret+8(FP)
 	RET
@@ -135,7 +135,7 @@
 	FMUL	F6, F6
 	WFMADB	V4, V1, V5, V1
 	FMOVD	48(R5), F7
-	WORD	$0xB3CD0013	//lgdr	%r1,%f3
+	LGDR	F3, R1
 	FMOVD	24(R5), F5
 	WFMADB	V4, V7, V5, V7
 	FMOVD	16(R5), F5
@@ -157,7 +157,7 @@
 	WORD	$0xEC21000F	//risbgn	%r2,%r1,64-64+0,64-64+0+16-1,64-0-16
 	BYTE	$0x30
 	BYTE	$0x59
-	WORD	$0xB3C10002	//ldgr	%f0,%r2
+	LDGR	R2, F0
 	FMADD	F0, F4, F0
 	MOVD	$·expx4ff<>+0(SB), R3
 	FMOVD	0(R3), F2
@@ -173,7 +173,7 @@
 	WORD	$0xEC21000F	//risbgn	%r2,%r1,64-64+0,64-64+0+16-1,64-0-16
 	BYTE	$0x30
 	BYTE	$0x59
-	WORD	$0xB3C10002	//ldgr	%f0,%r2
+	LDGR	R2, F0
 	FMADD	F0, F4, F0
 	MOVD	$·expx2ff<>+0(SB), R3
 	FMOVD	0(R3), F2
diff --git a/src/math/expm1_386.s b/src/math/expm1_386.s
index c1392cd..d020296 100644
--- a/src/math/expm1_386.s
+++ b/src/math/expm1_386.s
@@ -8,7 +8,7 @@
 TEXT ·Expm1(SB),NOSPLIT,$0
 	FLDLN2               // F0=log(2) = 1/log2(e) ~ 0.693147
 	FMOVD   x+0(FP), F0  // F0=x, F1=1/log2(e)
-	FABS                 // F0=|x|, F1=1/log2(e) 
+	FABS                 // F0=|x|, F1=1/log2(e)
 	FUCOMPP F0, F1       // compare F0 to F1
 	FSTSW   AX
 	SAHF
@@ -36,7 +36,7 @@
 	FSCALE                // F0=e**x, F1=int(x*log2(e))
 	FMOVDP  F0, F1        // F0=e**x
 	FLD1                  // F0=1, F1=e**x
-	FSUBDP  F0, F1        // F0=e**x-1 
+	FSUBDP  F0, F1        // F0=e**x-1
 	FMOVDP  F0, ret+8(FP)
 	RET
 not_finite:
diff --git a/src/math/expm1_s390x.s b/src/math/expm1_s390x.s
index 22e5eb1..c7c793b 100644
--- a/src/math/expm1_s390x.s
+++ b/src/math/expm1_s390x.s
@@ -89,7 +89,7 @@
 	FMADD	F2, F0, F6
 	WFMADB	V0, V5, V3, V5
 	WFMDB	V0, V0, V2
-	WORD	$0xB3CD0011	//lgdr	%r1,%f1
+	LGDR	F1, R1
 	WFMADB	V6, V2, V5, V6
 	FMOVD	40(R5), F3
 	FMOVD	32(R5), F5
@@ -108,7 +108,7 @@
 	FMADD	F4, F0, F0
 	SLD	$48, R1, R2
 	WFMSDB	V2, V0, V4, V0
-	WORD	$0xB3C10042	//ldgr	%f4,%r2
+	LDGR	R2, F4
 	WORD	$0xB3130000	//lcdbr	%f0,%f0
 	FSUB	F4, F6
 	WFMSDB	V0, V4, V6, V0
@@ -155,7 +155,7 @@
 	WFMADB	V1, V16, V3, V1
 	FMOVD	16(R5), F6
 	FMADD	F4, F1, F6
-	WORD	$0xB3CD0015	//lgdr	%r1,%f5
+	LGDR	F5, R1
 	WORD	$0xB3130066	//lcdbr	%f6,%f6
 	WORD	$0xEC3139BC	//risbg	%r3,%r1,57,128+60,3
 	BYTE	$0x03
@@ -171,7 +171,7 @@
 	WORD	$0xEC21000F	//risbgn	%r2,%r1,64-64+0,64-64+0+16-1,64-0-16
 	BYTE	$0x30
 	BYTE	$0x59
-	WORD	$0xB3C10002	//ldgr	%f0,%r2
+	LDGR	R2, F0
 	FMADD	F0, F4, F0
 	MOVD	$·expm1x4ff<>+0(SB), R3
 	FMOVD	0(R5), F4
@@ -189,7 +189,7 @@
 	WORD	$0xEC21000F	//risbgn	%r2,%r1,64-64+0,64-64+0+16-1,64-0-16
 	BYTE	$0x30
 	BYTE	$0x59
-	WORD	$0xB3C10002	//ldgr	%f0,%r2
+	LDGR	R2, F0
 	FMADD	F0, F4, F0
 	MOVD	$·expm1x2ff<>+0(SB), R3
 	FMOVD	0(R5), F4
diff --git a/src/math/export_test.go b/src/math/export_test.go
index 368308e..53d9205 100644
--- a/src/math/export_test.go
+++ b/src/math/export_test.go
@@ -9,3 +9,6 @@
 var Exp2Go = exp2
 var HypotGo = hypot
 var SqrtGo = sqrt
+var TrigReduce = trigReduce
+
+const ReduceThreshold = reduceThreshold
diff --git a/src/math/huge_test.go b/src/math/huge_test.go
new file mode 100644
index 0000000..0b45dbf
--- /dev/null
+++ b/src/math/huge_test.go
@@ -0,0 +1,99 @@
+// Copyright 2018 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.
+
+// Disabled for s390x because it uses assembly routines that are not
+// accurate for huge arguments.
+
+// +build !s390x
+
+package math_test
+
+import (
+	. "math"
+	"testing"
+)
+
+// Inputs to test trig_reduce
+var trigHuge = []float64{
+	1 << 120,
+	1 << 240,
+	1 << 480,
+	1234567891234567 << 180,
+	1234567891234567 << 300,
+	MaxFloat64,
+}
+
+// Results for trigHuge[i] calculated with https://github.com/robpike/ivy
+// using 4096 bits of working precision.   Values requiring less than
+// 102 decimal digits (1 << 120, 1 << 240, 1 << 480, 1234567891234567 << 180)
+// were confirmed via https://keisan.casio.com/
+var cosHuge = []float64{
+	-0.92587902285483787,
+	0.93601042593353793,
+	-0.28282777640193788,
+	-0.14616431394103619,
+	-0.79456058210671406,
+	-0.99998768942655994,
+}
+
+var sinHuge = []float64{
+	0.37782010936075202,
+	-0.35197227524865778,
+	0.95917070894368716,
+	0.98926032637023618,
+	-0.60718488235646949,
+	0.00496195478918406,
+}
+
+var tanHuge = []float64{
+	-0.40806638884180424,
+	-0.37603456702698076,
+	-3.39135965054779932,
+	-6.76813854009065030,
+	0.76417695016604922,
+	-0.00496201587444489,
+}
+
+// Check that trig values of huge angles return accurate results.
+// This confirms that argument reduction works for very large values
+// up to MaxFloat64.
+func TestHugeCos(t *testing.T) {
+	for i := 0; i < len(trigHuge); i++ {
+		f1 := cosHuge[i]
+		f2 := Cos(trigHuge[i])
+		if !close(f1, f2) {
+			t.Errorf("Cos(%g) = %g, want %g", trigHuge[i], f2, f1)
+		}
+	}
+}
+
+func TestHugeSin(t *testing.T) {
+	for i := 0; i < len(trigHuge); i++ {
+		f1 := sinHuge[i]
+		f2 := Sin(trigHuge[i])
+		if !close(f1, f2) {
+			t.Errorf("Sin(%g) = %g, want %g", trigHuge[i], f2, f1)
+		}
+	}
+}
+
+func TestHugeSinCos(t *testing.T) {
+	for i := 0; i < len(trigHuge); i++ {
+		f1, g1 := sinHuge[i], cosHuge[i]
+		f2, g2 := Sincos(trigHuge[i])
+		if !close(f1, f2) || !close(g1, g2) {
+			t.Errorf("Sincos(%g) = %g, %g, want %g, %g", trigHuge[i], f2, g2, f1, g1)
+		}
+	}
+}
+
+func TestHugeTan(t *testing.T) {
+	for i := 0; i < len(trigHuge); i++ {
+		f1 := tanHuge[i]
+		f2 := Tan(trigHuge[i])
+		if !close(f1, f2) {
+			t.Errorf("Tan(%g) = %g, want %g", trigHuge[i], f2, f1)
+		}
+	}
+}
diff --git a/src/math/log1p.go b/src/math/log1p.go
index b128a16..c4ec61b 100644
--- a/src/math/log1p.go
+++ b/src/math/log1p.go
@@ -151,12 +151,13 @@
 			u = 1.0 + x
 			iu = Float64bits(u)
 			k = int((iu >> 52) - 1023)
+			// correction term
 			if k > 0 {
 				c = 1.0 - (u - x)
 			} else {
-				c = x - (u - 1.0) // correction term
-				c /= u
+				c = x - (u - 1.0)
 			}
+			c /= u
 		} else {
 			u = x
 			iu = Float64bits(u)
diff --git a/src/math/log1p_s390x.s b/src/math/log1p_s390x.s
index c7e9860..ba4933d 100644
--- a/src/math/log1p_s390x.s
+++ b/src/math/log1p_s390x.s
@@ -96,7 +96,7 @@
 	MOVD	$·log1pc5<>+0(SB), R1
 	VLEG	$0, 0(R1), V16
 	MOVD	R2, R5
-	WORD	$0xB3CD0034	//lgdr	%r3,%f4
+	LGDR	F4, R3
 	WORD	$0xC0190006	//iilf	%r1,425983
 	BYTE	$0x7F
 	BYTE	$0xFF
@@ -118,7 +118,7 @@
 	MOVD	$·log1pxzero<>+0(SB), R1
 	FMOVD	0(R1), F2
 	BVS	LEXITTAGlog1p
-	WORD	$0xB3130044
+	WORD	$0xB3130044	// lcdbr %f4,%f4
 	WFCEDBS	V2, V4, V6
 	BEQ	L9
 	WFCHDBS	V4, V2, V2
@@ -129,11 +129,11 @@
 	RET
 
 L8:
-	WORD	$0xB3C10022	//ldgr	%f2,%r2
+	LDGR	R2, F2
 	FSUB	F4, F3
 	FMADD	F2, F4, F1
 	MOVD	$·log1pc4<>+0(SB), R2
-	WORD	$0xB3130041
+	WORD	$0xB3130041	// lcdbr %f4,%f1
 	FMOVD	0(R2), F7
 	FSUB	F3, F0
 	MOVD	$·log1pc3<>+0(SB), R2
@@ -164,7 +164,7 @@
 	FMOVD	0(R3), F2
 	WFMADB	V0, V6, V1, V0
 	MOVD	$·log1pyout<>+0(SB), R1
-	WORD	$0xB3C10065	//ldgr	%f6,%r5
+	LDGR	R5, F6
 	FMOVD	0(R1), F4
 	WFMSDB	V2, V6, V4, V2
 	MOVD	$·log1pxl2<>+0(SB), R1
diff --git a/src/math/log_s390x.s b/src/math/log_s390x.s
index 3e24ca7..7bcfdfc 100644
--- a/src/math/log_s390x.s
+++ b/src/math/log_s390x.s
@@ -63,7 +63,7 @@
 	FMOVD	x+0(FP), F0
 	MOVD	$·logrodataL21<>+0(SB), R9
 	MOVH	$0x8006, R4
-	WORD	$0xB3CD0010	//lgdr	%r1,%f0
+	LGDR	F0, R1
 	MOVD	$0x3FF0000000000000, R6
 	SRAD	$48, R1, R1
 	MOVD	$0x40F03E8000000000, R8
@@ -91,7 +91,7 @@
 	BLEU	L3
 L15:
 	FMUL	F2, F0
-	WORD	$0xB3CD0010	//lgdr	%r1,%f0
+	LGDR	F0, R1
 	SRAD	$48, R1, R1
 	SUBW	R1, R0, R2
 	SUBW	R1, R12, R3
@@ -114,7 +114,7 @@
 	MOVH	$0x7FEF, R1
 	CMPW	R5, R1
 	BGT	L1
-	WORD	$0xB3C10026	//ldgr	%f2,%r6
+	LDGR	R6, F2
 	FMUL	F2, F0
 	WORD	$0xEC4439BB	//risbg	%r4,%r4,57,128+59,3
 	BYTE	$0x03
@@ -148,14 +148,14 @@
 	WFMADB	V6, V4, V1, V4
 	FMOVD	8(R4), F1
 	WFMADB	V0, V2, V4, V2
-	WORD	$0xB3C10048	//ldgr	%f4,%r8
+	LDGR	R8, F4
 	WFMADB	V6, V2, V0, V2
 	WORD	$0xED401000	//msdb	%f1,%f4,0(%r1)
 	BYTE	$0x10
 	BYTE	$0x1F
 	MOVD	·logxl2<>+0(SB), R1
 	WORD	$0xB3130001	//lcdbr	%f0,%f1
-	WORD	$0xB3C10041	//ldgr	%f4,%r1
+	LDGR	R1, F4
 	WFMADB	V0, V4, V2, V0
 L1:
 	FMOVD	F0, ret+8(FP)
diff --git a/src/math/mod.go b/src/math/mod.go
index e1a414e..7efc018 100644
--- a/src/math/mod.go
+++ b/src/math/mod.go
@@ -24,16 +24,12 @@
 	if y == 0 || IsInf(x, 0) || IsNaN(x) || IsNaN(y) {
 		return NaN()
 	}
-	if y < 0 {
-		y = -y
-	}
+	y = Abs(y)
 
 	yfr, yexp := Frexp(y)
-	sign := false
 	r := x
 	if x < 0 {
 		r = -x
-		sign = true
 	}
 
 	for r >= y {
@@ -43,7 +39,7 @@
 		}
 		r = r - Ldexp(y, rexp-yexp)
 	}
-	if sign {
+	if x < 0 {
 		r = -r
 	}
 	return r
diff --git a/src/math/pow.go b/src/math/pow.go
index 336193b..2219a90 100644
--- a/src/math/pow.go
+++ b/src/math/pow.go
@@ -83,13 +83,7 @@
 		return 1 / Sqrt(x)
 	}
 
-	absy := y
-	flip := false
-	if absy < 0 {
-		absy = -absy
-		flip = true
-	}
-	yi, yf := Modf(absy)
+	yi, yf := Modf(Abs(y))
 	if yf != 0 && x < 0 {
 		return NaN()
 	}
@@ -147,9 +141,9 @@
 	}
 
 	// ans = a1*2**ae
-	// if flip { ans = 1 / ans }
+	// if y < 0 { ans = 1 / ans }
 	// but in the opposite order
-	if flip {
+	if y < 0 {
 		a1 = 1 / a1
 		ae = -ae
 	}
diff --git a/src/math/pow_s390x.s b/src/math/pow_s390x.s
index fd19617..754b119 100644
--- a/src/math/pow_s390x.s
+++ b/src/math/pow_s390x.s
@@ -297,7 +297,7 @@
 	FMOVD	x+0(FP), F0
 	FMOVD	y+8(FP), F2
 	MOVD	$·powrodataL51<>+0(SB), R9
-	WORD	$0xB3CD0030	//lgdr	%r3,%f0
+	LGDR	F0, R3
 	WORD	$0xC0298009	//iilf	%r2,2148095317
 	BYTE	$0x55
 	BYTE	$0x55
@@ -340,7 +340,7 @@
 	BYTE	$0x24
 	FMOVD	0(R2), F6
 	FSUBS	F1, F3
-	WORD	$0xB3C10018	//ldgr	%f1,%r8
+	LDGR	R8, F1
 	WFMSDB	V4, V1, V6, V4
 	FMOVD	152(R9), F6
 	WFMDB	V4, V4, V7
@@ -387,7 +387,7 @@
 	WFMSDB	V2, V3, V5, V3
 	VLEG	$0, 48(R9), V18
 	WFADB	V3, V5, V6
-	WORD	$0xB3CD0023	//lgdr	%r2,%f3
+	LGDR	F3, R2
 	WFMSDB	V2, V16, V6, V16
 	FMOVD	40(R9), F1
 	WFMADB	V2, V4, V16, V4
@@ -410,8 +410,8 @@
 	BYTE	$0x30
 	BYTE	$0x59
 	WFMADB	V4, V1, V3, V4
-	WORD	$0xB3CD0026	//lgdr	%r2,%f6
-	WORD	$0xB3C10015	//ldgr	%f1,%r5
+	LGDR	F6, R2
+	LDGR	R5, F1
 	SRAD	$48, R2, R2
 	FMADD	F1, F4, F1
 	RLL	$16, R2, R2
@@ -452,7 +452,7 @@
 	WORD	$0xEC1520BF	//risbgn	%r1,%r5,64-32,128+63,64+0+32
 	BYTE	$0x60
 	BYTE	$0x59
-	WORD	$0xB3CD0026	//lgdr	%r2,%f6
+	LGDR	F6, R2
 	MOVD	$powiadd<>+0(SB), R3
 	WORD	$0xEC223CBC	//risbg	%r2,%r2,60,128+60,64-60
 	BYTE	$0x04
@@ -461,7 +461,7 @@
 	WORD	$0xEC51001F	//risbgn	%r5,%r1,64-64+0,64-64+0+32-1,64-0-32
 	BYTE	$0x20
 	BYTE	$0x59
-	WORD	$0xB3C10015	//ldgr	%f1,%r5
+	LDGR	R5, F1
 	FMADD	F1, F4, F1
 	MOVD	$powxscale<>+0(SB), R1
 	WORD	$0xED121000	//mdb	%f1,0(%r2,%r1)
@@ -486,7 +486,7 @@
 	WORD	$0xC0298009	//iilf	%r2,2148095317
 	BYTE	$0x55
 	BYTE	$0x55
-	WORD	$0xB3CD0034	//lgdr	%r3,%f4
+	LGDR	F4, R3
 	WORD	$0xEC3320BF	//risbgn	%r3,%r3,64-32,128+63,64+0+32
 	BYTE	$0x60
 	BYTE	$0x59
@@ -566,11 +566,11 @@
 	BVS	L49
 L16:
 	MOVD	·pow_xnan<>+0(SB), R1
-	WORD	$0xB3C10001	//ldgr	%f0,%r1
+	LDGR	R1, F0
 	WFMDB	V4, V0, V1
 	BR	L1
 L48:
-	WORD	$0xB3CD0030	//lgdr	%r3,%f0
+	LGDR	F0, R3
 	WORD	$0xEC1320BF	//risbgn	%r1,%r3,64-32,128+63,64+0+32
 	BYTE	$0x60
 	BYTE	$0x59
diff --git a/src/math/signbit.go b/src/math/signbit.go
index 670cc1a..f6e61d6 100644
--- a/src/math/signbit.go
+++ b/src/math/signbit.go
@@ -4,7 +4,7 @@
 
 package math
 
-// Signbit returns true if x is negative or negative zero.
+// Signbit reports whether x is negative or negative zero.
 func Signbit(x float64) bool {
 	return Float64bits(x)&(1<<63) != 0
 }
diff --git a/src/math/sin.go b/src/math/sin.go
index 929cac3..cc8b136 100644
--- a/src/math/sin.go
+++ b/src/math/sin.go
@@ -118,10 +118,9 @@
 
 func cos(x float64) float64 {
 	const (
-		PI4A = 7.85398125648498535156E-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
-		PI4B = 3.77489470793079817668E-8                             // 0x3e64442d00000000,
-		PI4C = 2.69515142907905952645E-15                            // 0x3ce8469898cc5170,
-		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+		PI4A = 7.85398125648498535156E-1  // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668E-8  // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
 	)
 	// special cases
 	switch {
@@ -133,15 +132,23 @@
 	sign := false
 	x = Abs(x)
 
-	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
-	y := float64(j)      // integer part of x/(Pi/4), as float
+	var j uint64
+	var y, z float64
+	if x >= reduceThreshold {
+		j, z = trigReduce(x)
+	} else {
+		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
+		y = float64(j)           // integer part of x/(Pi/4), as float
 
-	// map zeros to origin
-	if j&1 == 1 {
-		j++
-		y++
+		// map zeros to origin
+		if j&1 == 1 {
+			j++
+			y++
+		}
+		j &= 7                               // octant modulo 2Pi radians (360 degrees)
+		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	}
-	j &= 7 // octant modulo 2Pi radians (360 degrees)
+
 	if j > 3 {
 		j -= 4
 		sign = !sign
@@ -150,7 +157,6 @@
 		sign = !sign
 	}
 
-	z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	zz := z * z
 	if j == 1 || j == 2 {
 		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
@@ -173,10 +179,9 @@
 
 func sin(x float64) float64 {
 	const (
-		PI4A = 7.85398125648498535156E-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
-		PI4B = 3.77489470793079817668E-8                             // 0x3e64442d00000000,
-		PI4C = 2.69515142907905952645E-15                            // 0x3ce8469898cc5170,
-		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+		PI4A = 7.85398125648498535156E-1  // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668E-8  // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
 	)
 	// special cases
 	switch {
@@ -193,22 +198,27 @@
 		sign = true
 	}
 
-	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
-	y := float64(j)      // integer part of x/(Pi/4), as float
+	var j uint64
+	var y, z float64
+	if x >= reduceThreshold {
+		j, z = trigReduce(x)
+	} else {
+		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
+		y = float64(j)           // integer part of x/(Pi/4), as float
 
-	// map zeros to origin
-	if j&1 == 1 {
-		j++
-		y++
+		// map zeros to origin
+		if j&1 == 1 {
+			j++
+			y++
+		}
+		j &= 7                               // octant modulo 2Pi radians (360 degrees)
+		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	}
-	j &= 7 // octant modulo 2Pi radians (360 degrees)
 	// reflect in x axis
 	if j > 3 {
 		sign = !sign
 		j -= 4
 	}
-
-	z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	zz := z * z
 	if j == 1 || j == 2 {
 		y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
diff --git a/src/math/sin_386.s b/src/math/sin_386.s
index 9d605a1..cf7679d 100644
--- a/src/math/sin_386.s
+++ b/src/math/sin_386.s
@@ -6,42 +6,8 @@
 
 // func Cos(x float64) float64
 TEXT ·Cos(SB),NOSPLIT,$0
-	FMOVD   x+0(FP), F0  // F0=x
-	FCOS                 // F0=cos(x) if -2**63 < x < 2**63
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     3(PC)        // jump if x outside range
-	FMOVDP  F0, ret+8(FP)
-	RET
-	FLDPI                // F0=Pi, F1=x
-	FADDD   F0, F0       // F0=2*Pi, F1=x
-	FXCHD   F0, F1       // F0=x, F1=2*Pi
-	FPREM1               // F0=reduced_x, F1=2*Pi
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     -3(PC)       // jump if reduction incomplete
-	FMOVDP  F0, F1       // F0=reduced_x
-	FCOS                 // F0=cos(reduced_x)
-	FMOVDP  F0, ret+8(FP)
-	RET
-	
+	JMP ·cos(SB)
+
 // func Sin(x float64) float64
 TEXT ·Sin(SB),NOSPLIT,$0
-	FMOVD   x+0(FP), F0  // F0=x
-	FSIN                 // F0=sin(x) if -2**63 < x < 2**63
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     3(PC)        // jump if x outside range
-	FMOVDP  F0, ret+8(FP)
-	RET
-	FLDPI                // F0=Pi, F1=x
-	FADDD   F0, F0       // F0=2*Pi, F1=x
-	FXCHD   F0, F1       // F0=x, F1=2*Pi
-	FPREM1               // F0=reduced_x, F1=2*Pi
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     -3(PC)       // jump if reduction incomplete
-	FMOVDP  F0, F1       // F0=reduced_x
-	FSIN                 // F0=sin(reduced_x)
-	FMOVDP  F0, ret+8(FP)
-	RET
+	JMP ·sin(SB)
diff --git a/src/math/sincos.go b/src/math/sincos.go
index 3ae193a..c002db6 100644
--- a/src/math/sincos.go
+++ b/src/math/sincos.go
@@ -2,8 +2,6 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
-// +build !386
-
 package math
 
 // Coefficients _sin[] and _cos[] are found in pkg/math/sin.go.
@@ -16,10 +14,9 @@
 //	Sincos(NaN) = NaN, NaN
 func Sincos(x float64) (sin, cos float64) {
 	const (
-		PI4A = 7.85398125648498535156E-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
-		PI4B = 3.77489470793079817668E-8                             // 0x3e64442d00000000,
-		PI4C = 2.69515142907905952645E-15                            // 0x3ce8469898cc5170,
-		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+		PI4A = 7.85398125648498535156E-1  // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668E-8  // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
 	)
 	// special cases
 	switch {
@@ -36,14 +33,21 @@
 		sinSign = true
 	}
 
-	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
-	y := float64(j)      // integer part of x/(Pi/4), as float
+	var j uint64
+	var y, z float64
+	if x >= reduceThreshold {
+		j, z = trigReduce(x)
+	} else {
+		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
+		y = float64(j)           // integer part of x/(Pi/4), as float
 
-	if j&1 == 1 { // map zeros to origin
-		j++
-		y++
+		if j&1 == 1 { // map zeros to origin
+			j++
+			y++
+		}
+		j &= 7                               // octant modulo 2Pi radians (360 degrees)
+		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	}
-	j &= 7     // octant modulo 2Pi radians (360 degrees)
 	if j > 3 { // reflect in x axis
 		j -= 4
 		sinSign, cosSign = !sinSign, !cosSign
@@ -52,7 +56,6 @@
 		cosSign = !cosSign
 	}
 
-	z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
 	zz := z * z
 	cos = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
 	sin = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
diff --git a/src/math/sincos_386.go b/src/math/sincos_386.go
deleted file mode 100644
index 38bb050..0000000
--- a/src/math/sincos_386.go
+++ /dev/null
@@ -1,13 +0,0 @@
-// 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
-
-// Sincos returns Sin(x), Cos(x).
-//
-// Special cases are:
-//	Sincos(±0) = ±0, 1
-//	Sincos(±Inf) = NaN, NaN
-//	Sincos(NaN) = NaN, NaN
-func Sincos(x float64) (sin, cos float64)
diff --git a/src/math/sincos_386.s b/src/math/sincos_386.s
deleted file mode 100644
index f700a4f..0000000
--- a/src/math/sincos_386.s
+++ /dev/null
@@ -1,28 +0,0 @@
-// Copyright 2010 The Go Authors. All rights reserved.
-// Use of this source code is governed by a BSD-style
-// license that can be found in the LICENSE file.
-
-#include "textflag.h"
-
-// func Sincos(x float64) (sin, cos float64)
-TEXT ·Sincos(SB),NOSPLIT,$0
-	FMOVD   x+0(FP), F0  // F0=x
-	FSINCOS              // F0=cos(x), F1=sin(x) if -2**63 < x < 2**63
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     4(PC)        // jump if x outside range
-	FMOVDP  F0, cos+16(FP) // F0=sin(x)
-	FMOVDP  F0, sin+8(FP)
-	RET
-	FLDPI                // F0=Pi, F1=x
-	FADDD   F0, F0       // F0=2*Pi, F1=x
-	FXCHD   F0, F1       // F0=x, F1=2*Pi
-	FPREM1               // F0=reduced_x, F1=2*Pi
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     -3(PC)       // jump if reduction incomplete
-	FMOVDP  F0, F1       // F0=reduced_x
-	FSINCOS              // F0=cos(reduced_x), F1=sin(reduced_x)
-	FMOVDP  F0, cos+16(FP) // F0=sin(reduced_x)
-	FMOVDP  F0, sin+8(FP)
-	RET
diff --git a/src/math/sinh.go b/src/math/sinh.go
index 39e7c20..573a37e 100644
--- a/src/math/sinh.go
+++ b/src/math/sinh.go
@@ -43,7 +43,7 @@
 	}
 
 	var temp float64
-	switch true {
+	switch {
 	case x > 21:
 		temp = Exp(x) * 0.5
 
diff --git a/src/math/sqrt_386.s b/src/math/sqrt_386.s
index 402d152..5a5c33a 100644
--- a/src/math/sqrt_386.s
+++ b/src/math/sqrt_386.s
@@ -4,7 +4,7 @@
 
 #include "textflag.h"
 
-// func Sqrt(x float64) float64	
+// func Sqrt(x float64) float64
 TEXT ·Sqrt(SB),NOSPLIT,$0
 	FMOVD   x+0(FP),F0
 	FSQRT
diff --git a/src/math/sqrt_arm.s b/src/math/sqrt_arm.s
index deb6712..ffc7d10 100644
--- a/src/math/sqrt_arm.s
+++ b/src/math/sqrt_arm.s
@@ -4,7 +4,7 @@
 
 #include "textflag.h"
 
-// func Sqrt(x float64) float64	
+// func Sqrt(x float64) float64
 TEXT ·Sqrt(SB),NOSPLIT,$0
 	MOVB	runtime·goarm(SB), R11
 	CMP	$5, R11
diff --git a/src/math/tan.go b/src/math/tan.go
index aa2fb37..0d5394c 100644
--- a/src/math/tan.go
+++ b/src/math/tan.go
@@ -83,10 +83,9 @@
 
 func tan(x float64) float64 {
 	const (
-		PI4A = 7.85398125648498535156E-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
-		PI4B = 3.77489470793079817668E-8                             // 0x3e64442d00000000,
-		PI4C = 2.69515142907905952645E-15                            // 0x3ce8469898cc5170,
-		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+		PI4A = 7.85398125648498535156E-1  // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668E-8  // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
 	)
 	// special cases
 	switch {
@@ -102,17 +101,22 @@
 		x = -x
 		sign = true
 	}
+	var j uint64
+	var y, z float64
+	if x >= reduceThreshold {
+		j, z = trigReduce(x)
+	} else {
+		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
+		y = float64(j)           // integer part of x/(Pi/4), as float
 
-	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
-	y := float64(j)      // integer part of x/(Pi/4), as float
+		/* map zeros and singularities to origin */
+		if j&1 == 1 {
+			j++
+			y++
+		}
 
-	/* map zeros and singularities to origin */
-	if j&1 == 1 {
-		j++
-		y++
+		z = ((x - y*PI4A) - y*PI4B) - y*PI4C
 	}
-
-	z := ((x - y*PI4A) - y*PI4B) - y*PI4C
 	zz := z * z
 
 	if zz > 1e-14 {
diff --git a/src/math/tan_386.s b/src/math/tan_386.s
index cb65a3f..4e44c26 100644
--- a/src/math/tan_386.s
+++ b/src/math/tan_386.s
@@ -6,23 +6,4 @@
 
 // func Tan(x float64) float64
 TEXT ·Tan(SB),NOSPLIT,$0
-	FMOVD   x+0(FP), F0  // F0=x
-	FPTAN                // F0=1, F1=tan(x) if -2**63 < x < 2**63
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     4(PC)        // jump if x outside range
-	FMOVDP  F0, F0       // F0=tan(x)
-	FMOVDP  F0, ret+8(FP)
-	RET
-	FLDPI                // F0=Pi, F1=x
-	FADDD   F0, F0       // F0=2*Pi, F1=x
-	FXCHD   F0, F1       // F0=x, F1=2*Pi
-	FPREM1               // F0=reduced_x, F1=2*Pi
-	FSTSW   AX           // AX=status word
-	ANDW    $0x0400, AX
-	JNE     -3(PC)       // jump if reduction incomplete
-	FMOVDP  F0, F1       // F0=reduced_x
-	FPTAN                // F0=1, F1=tan(reduced_x)
-	FMOVDP  F0, F0       // F0=tan(reduced_x)
-	FMOVDP  F0, ret+8(FP)
-	RET
+	JMP     ·tan(SB)
diff --git a/src/math/tan_s390x.s b/src/math/tan_s390x.s
index 7b05ba0..b6e2295 100644
--- a/src/math/tan_s390x.s
+++ b/src/math/tan_s390x.s
@@ -68,7 +68,7 @@
 	WFMADB	V4, V3, V2, V4
 	FMUL	F2, F2
 	VLEG	$0, 48(R5), V18
-	WORD	$0xB3CD0016	//lgdr	%r1,%f6
+	LGDR	F6, R1
 	FMOVD	40(R5), F5
 	FMOVD	32(R5), F3
 	FMADD	F1, F2, F3
@@ -82,7 +82,7 @@
 	WFLCDB	V4, V16
 	WFMADB	V2, V5, V18, V5
 	WFMADB	V1, V0, V7, V0
-	WORD	$0xA7110001	//tmll	%r1,1
+	TMLL	R1, $1
 	WFMADB	V1, V5, V3, V1
 	BNE	L12
 	WFDDB	V0, V1, V0
diff --git a/src/math/trig_reduce.go b/src/math/trig_reduce.go
new file mode 100644
index 0000000..6f8eaba
--- /dev/null
+++ b/src/math/trig_reduce.go
@@ -0,0 +1,94 @@
+// Copyright 2018 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
+
+import (
+	"math/bits"
+)
+
+// reduceThreshold is the maximum value where the reduction using Pi/4
+// in 3 float64 parts still gives accurate results.  Above this
+// threshold Payne-Hanek range reduction must be used.
+const reduceThreshold = (1 << 52) / (4 / Pi)
+
+// trigReduce implements Payne-Hanek range reduction by Pi/4
+// for x > 0.  It returns the integer part mod 8 (j) and
+// the fractional part (z) of x / (Pi/4).
+// The implementation is based on:
+// "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit"
+// K. C. Ng et al, March 24, 1992
+// The simulated multi-precision calculation of x*B uses 64-bit integer arithmetic.
+func trigReduce(x float64) (j uint64, z float64) {
+	const PI4 = Pi / 4
+	if x < PI4 {
+		return 0, x
+	}
+	// Extract out the integer and exponent such that,
+	// x = ix * 2 ** exp.
+	ix := Float64bits(x)
+	exp := int(ix>>shift&mask) - bias - shift
+	ix &^= mask << shift
+	ix |= 1 << shift
+	// Use the exponent to extract the 3 appropriate uint64 digits from mPi4,
+	// B ~ (z0, z1, z2), such that the product leading digit has the exponent -61.
+	// Note, exp >= -53 since x >= PI4 and exp < 971 for maximum float64.
+	digit, bitshift := uint(exp+61)/64, uint(exp+61)%64
+	z0 := (mPi4[digit] << bitshift) | (mPi4[digit+1] >> (64 - bitshift))
+	z1 := (mPi4[digit+1] << bitshift) | (mPi4[digit+2] >> (64 - bitshift))
+	z2 := (mPi4[digit+2] << bitshift) | (mPi4[digit+3] >> (64 - bitshift))
+	// Multiply mantissa by the digits and extract the upper two digits (hi, lo).
+	z2hi, _ := bits.Mul64(z2, ix)
+	z1hi, z1lo := bits.Mul64(z1, ix)
+	z0lo := z0 * ix
+	lo, c := bits.Add64(z1lo, z2hi, 0)
+	hi, _ := bits.Add64(z0lo, z1hi, c)
+	// The top 3 bits are j.
+	j = hi >> 61
+	// Extract the fraction and find its magnitude.
+	hi = hi<<3 | lo>>61
+	lz := uint(bits.LeadingZeros64(hi))
+	e := uint64(bias - (lz + 1))
+	// Clear implicit mantissa bit and shift into place.
+	hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1)))
+	hi >>= 64 - shift
+	// Include the exponent and convert to a float.
+	hi |= e << shift
+	z = Float64frombits(hi)
+	// Map zeros to origin.
+	if j&1 == 1 {
+		j++
+		j &= 7
+		z--
+	}
+	// Multiply the fractional part by pi/4.
+	return j, z * PI4
+}
+
+// mPi4 is the binary digits of 4/pi as a uint64 array,
+// that is, 4/pi = Sum mPi4[i]*2^(-64*i)
+// 19 64-bit digits and the leading one bit give 1217 bits
+// of precision to handle the largest possible float64 exponent.
+var mPi4 = [...]uint64{
+	0x0000000000000001,
+	0x45f306dc9c882a53,
+	0xf84eafa3ea69bb81,
+	0xb6c52b3278872083,
+	0xfca2c757bd778ac3,
+	0x6e48dc74849ba5c0,
+	0x0c925dd413a32439,
+	0xfc3bd63962534e7d,
+	0xd1046bea5d768909,
+	0xd338e04d68befc82,
+	0x7323ac7306a673e9,
+	0x3908bf177bf25076,
+	0x3ff12fffbc0b301f,
+	0xde5e2316b414da3e,
+	0xda6cfd9e4f96136e,
+	0x9e8c7ecd3cbfd45a,
+	0xea4f758fd7cbe2f6,
+	0x7a0e73ef14a525d4,
+	0xd7f6bf623f1aba10,
+	0xac06608df8f6d757,
+}
diff --git a/src/math/unsafe.go b/src/math/unsafe.go
index 5ae6742..e59f50c 100644
--- a/src/math/unsafe.go
+++ b/src/math/unsafe.go
@@ -6,16 +6,24 @@
 
 import "unsafe"
 
-// Float32bits returns the IEEE 754 binary representation of f.
+// Float32bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position.
+// Float32bits(Float32frombits(x)) == x.
 func Float32bits(f float32) uint32 { return *(*uint32)(unsafe.Pointer(&f)) }
 
-// Float32frombits returns the floating point number corresponding
-// to the IEEE 754 binary representation b.
+// Float32frombits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// Float32frombits(Float32bits(x)) == x.
 func Float32frombits(b uint32) float32 { return *(*float32)(unsafe.Pointer(&b)) }
 
-// Float64bits returns the IEEE 754 binary representation of f.
+// Float64bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position,
+// and Float64bits(Float64frombits(x)) == x.
 func Float64bits(f float64) uint64 { return *(*uint64)(unsafe.Pointer(&f)) }
 
-// Float64frombits returns the floating point number corresponding
-// the IEEE 754 binary representation b.
+// Float64frombits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// Float64frombits(Float64bits(x)) == x.
 func Float64frombits(b uint64) float64 { return *(*float64)(unsafe.Pointer(&b)) }