From d6edc2c4dbd928ff929d3e95efb3cbf04fcc211b Mon Sep 17 00:00:00 2001
From: liu-du <duliujimmy@hotmail.com>
Date: Sun, 6 Nov 2022 15:26:03 +1100
Subject: [PATCH] math: improve math.Log to handle subnormal floating number
 correctly on amd64

The existing implementation for math.Log on amd64/s390x, math/log_amd64.S,
doesn't normalize subnormal IEEE 64-bit floating point number correctly.
This only occurs on amd64/s390x. The go implementation of log, which is used
on other architecture does handle subnormal floating number correctly.
This commit add normalize logic to math/log_amd64.S.

Fixes #56600
---
 src/math/all_test.go | 27 +++++++++++++++++++++++++++
 src/math/log_amd64.s | 11 +++++++++++
 2 files changed, 38 insertions(+)

diff --git a/src/math/all_test.go b/src/math/all_test.go
index 8d5e0ad439..3a9b498979 100644
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -1511,6 +1511,27 @@ var logSC = []float64{
 	NaN(),
 }
 
+var vflogBC = []float64{
+	SmallestNonzeroFloat64,
+	LargestSubnormalFloat64,
+	SmallestNormalFloat64,
+	MaxFloat64,
+	-SmallestNonzeroFloat64,
+	-LargestSubnormalFloat64,
+	-SmallestNormalFloat64,
+	-MaxFloat64,
+}
+var logBC = []float64{
+	-744.4400719213812,
+	-708.3964185322641,
+	-708.3964185322641,
+	709.782712893384,
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+}
+
 var vflogbSC = []float64{
 	Inf(-1),
 	0,
@@ -2717,6 +2738,12 @@ func TestLog(t *testing.T) {
 			t.Errorf("Log(%g) = %g, want %g", vflogSC[i], f, logSC[i])
 		}
 	}
+
+	for i := 0; i < len(vflogBC); i++ {
+		if f := Log(vflogBC[i]); !alike(logBC[i], f) {
+			t.Errorf("Log(%g) = %g, want %g", vflogBC[i], f, logBC[i])
+		}
+	}
 }
 
 func TestLogb(t *testing.T) {
diff --git a/src/math/log_amd64.s b/src/math/log_amd64.s
index d84091f23a..2ca376b213 100644
--- a/src/math/log_amd64.s
+++ b/src/math/log_amd64.s
@@ -14,6 +14,7 @@
 #define L5     1.818357216161805012e-01   // 0x3FC7466496CB03DE
 #define L6     1.531383769920937332e-01   // 0x3FC39A09D078C69F
 #define L7     1.479819860511658591e-01   // 0x3FC2F112DF3E5244
+#define Two52  4.503599627370496e15       // 0x4330000000000000 (2**52)
 #define NaN    0x7FF8000000000001
 #define NegInf 0xFFF0000000000000
 #define PosInf 0x7FF0000000000000
@@ -32,7 +33,16 @@ TEXT ·archLog(SB),NOSPLIT,$0
 	CMPQ    AX, BX
 	JLE     isInfOrNaN
 	// f1, ki := math.Frexp(x); k := float64(ki)
+	MOVQ    $0, CX
 	MOVQ    BX, X0
+	ANDQ    BX, AX
+	JNE     isNormal
+	// f, exp = normalize(f)
+	MOVSD   $Two52, X1
+	MULSD   X1, X0
+	MOVQ    X0, BX
+	MOVQ    $-52, CX
+isNormal:
 	MOVQ    $0x000FFFFFFFFFFFFF, AX
 	MOVQ    AX, X2
 	ANDPD   X0, X2
@@ -41,6 +51,7 @@ TEXT ·archLog(SB),NOSPLIT,$0
 	SHRQ    $52, BX
 	ANDL    $0x7FF, BX
 	SUBL    $0x3FE, BX
+	ADDQ    CX, BX
 	XORPS   X1, X1 // break dependency for CVTSL2SD
 	CVTSL2SD BX, X1 // x1= k, x2= f1
 	// if f1 < math.Sqrt2/2 { k -= 1; f1 *= 2 }