math: fix math.log10() for -exclude @vlib/math/*.c.v (fix #23136) (#23140)

vlib/math/log.v

@@ -1,5 +1,10 @@
 module math
 
+const two54 = f64(1.80143985094819840000e+16)
+const ivln10 = f64(4.34294481903251816668e-01)
+const log10_2hi = f64(3.01029995663611771306e-01)
+const log10_2lo = f64(3.69423907715893078616e-13)
+
 // log_n returns log base b of x
 pub fn log_n(x f64, b f64) f64 {
 	y := log(x)
@@ -9,8 +14,43 @@ pub fn log_n(x f64, b f64) f64 {
 
 // log10 returns the decimal logarithm of x.
 // The special cases are the same as for log.
+// log10(10**N) = N  for N=0,1,...,22.
 pub fn log10(x f64) f64 {
-	return log(x) * (1.0 / ln10)
+	// https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_log10.c
+
+	mut x_ := x
+	mut hx := i64(f64_bits(x_))
+	mut k := i32(0)
+	if hx < i64(0x0010000000000000) {
+		// x < 2**-1022
+		if hx & 0x7fffffffffffffff == 0 {
+			return inf(-1) // log(+-0)=-inf
+		}
+		if hx < 0 {
+			return (x_ - x_) / (x_ - x_) // log(-#) = NaN
+		}
+		k = k - 54
+		x_ *= two54 // subnormal number, scale up x
+		hx = i64(f64_bits(x_))
+	}
+
+	// scale up resulted in a NaN number
+	if hx >= u64(0x7ff0000000000000) {
+		return x_ + x_
+	}
+
+	k = k + i32((u64((hx >> 52) - 1023)))
+	i := i32((u64(k) & 0x8000000000000000) >> 63)
+	hx = (hx & 0x000fffffffffffff) | (u64(0x3ff - i) << 52)
+	y := f64(k + i)
+	/*
+	if FIX_INT_FP_CONVERT_ZERO && y == 0.0 {
+		y = 0.0
+	}
+	*/
+	x_ = f64_from_bits(u64(hx))
+	z := y * log10_2lo + ivln10 * log(x_)
+	return z + y * log10_2hi
 }
 
 // log2 returns the binary logarithm of x.

vlib/math/log_test.v

@@ -9,9 +9,8 @@ fn test_log2_base() {
 }
 
 fn test_log10_base() {
-	// Note: `assert math.log10(10.0) == 1.0` currently fails, when the pure V implementation is used with
-	// `./v -exclude @vlib/math/*.c.v vlib/math/log_test.v`
-	assert math.veryclose(math.log10(10.0), 1.0)
+	assert math.log10(10.0) == 1.0
+	assert math.log10(0.00000000000000001) == -17.0
 }
 
 fn test_log1p_base() {

vlib/math/math_test.v

@@ -108,10 +108,10 @@ const log_ = [f64(1.605231462693062999102599e+00), 2.0462560018708770653153909e+
 const logb_ = [f64(2.0000000000000000e+00), 2.0000000000000000e+00, -2.0000000000000000e+00,
 	2.0000000000000000e+00, 3.0000000000000000e+00, 1.0000000000000000e+00, 2.0000000000000000e+00,
 	1.0000000000000000e+00, 0.0000000000000000e+00, 3.0000000000000000e+00]
-const log10_ = [f64(6.9714316642508290997617083e-01), 8.886776901739320576279124e-01,
-	-5.5770832400658929815908236e-01, 6.998900476822994346229723e-01, 9.8391002850684232013281033e-01,
-	4.6633031029295153334285302e-01, 7.1842557117242328821552533e-01, 4.3583479968917773161304553e-01,
-	2.6133617905227038228626834e-01, 9.3881606348649405716214241e-01]
+const log10_ = [f64(6.9714316642508289412205613e-01), 8.8867769017393205555066515e-01,
+	-5.5770832400658929815908236e-01, 6.998900476822994346229723e-01, 9.8391002850684228242528206e-01,
+	4.6633031029295152203317798e-01, 7.1842557117242322739514293e-01, 4.3583479968917770985825655e-01,
+	2.6133617905227035649318168e-01, 9.3881606348649405716214241e-01]
 const log1p_ = [f64(4.8590257759797794104158205e-02), 7.4540265965225865330849141e-02,
 	-2.7726407903942672823234024e-03, -5.1404917651627649094953380e-02,
 	9.1998280672258624681335010e-02, 2.8843762576593352865894824e-02, 5.0969534581863707268992645e-02,