translate netlib l2norm algorithm to asm and remove branches
Netlib algorithm reduces overflow while calculating the l2norm of a
vector.
Translated to asm while reducing branches in NaN and Inf checks.
Overflow protection is equivalent to the Netlib standard implementation.
MODIFY
floats/floats_test.go
floats/floats_test.go
@@ -36,7 +36,6 @@ func areSlicesSame(t *testing.T, truth, comp []float64, str string) {
break
}
}
-
}
if !ok {
t.Errorf(str+". Expected %v, returned %v", truth, comp)
@@ -96,7 +95,6 @@ func TestAddTo(t *testing.T) {
if !Panics(func() { AddTo(make([]float64, 3), make([]float64, 3), make([]float64, 2)) }) {
t.Errorf("Did not panic with length mismatch")
}
-
}
func TestAddConst(t *testing.T) {
@@ -208,7 +206,6 @@ func TestCumProd(t *testing.T) {
truth = []float64{}
CumProd(emptyReceiver, emptyReceiver)
areSlicesEqual(t, truth, emptyReceiver, "Wrong cumprod returned with empty receiver")
-
}
func TestCumSum(t *testing.T) {
@@ -231,7 +228,6 @@ func TestCumSum(t *testing.T) {
truth = []float64{}
CumSum(emptyReceiver, emptyReceiver)
areSlicesEqual(t, truth, emptyReceiver, "Wrong cumsum returned with empty receiver")
-
}
func TestDistance(t *testing.T) {
@@ -270,7 +266,6 @@ func TestDistance(t *testing.T) {
if !Panics(func() { Distance([]float64{}, norms, 1) }) {
t.Errorf("Did not panic with unequal lengths")
}
-
}
func TestDiv(t *testing.T) {
@@ -398,7 +393,7 @@ func TestEqualFunc(t *testing.T) {
}
func TestEqualsRelative(t *testing.T) {
- var equalityTests = []struct {
+ equalityTests := []struct {
a, b float64
tol float64
equal bool
@@ -518,7 +513,6 @@ func TestEqualsULP(t *testing.T) {
if EqualWithinULP(1, math.NaN(), 10) {
t.Errorf("NaN returned as equal")
}
-
}
func TestEqualLengths(t *testing.T) {
@@ -699,7 +693,6 @@ func TestLogSumExp(t *testing.T) {
if math.Abs(val-truth) > EqTolerance {
t.Errorf("Wrong logsumexp for values with negative infinity")
}
-
}
func TestMaxAndIdx(t *testing.T) {
@@ -1572,7 +1565,6 @@ func TestWithin(t *testing.T) {
t.Errorf("Case %v: Idx mismatch. Want: %v, got: %v", i, test.idx, idx)
}
}
-
}
func randomSlice(l int) []float64 {

MODIFY
internal/asm/f64/l2norm.go
internal/asm/f64/l2norm.go
@@ -6,33 +6,6 @@ package f64
import "math"
-// L2NormUnitary is the level 2 norm of x.
-func L2NormUnitary(x []float64) (sum float64) {
- var scale float64
- sumSquares := 1.0
- for _, v := range x {
- if v == 0 {
- continue
- }
- absxi := math.Abs(v)
- if math.IsNaN(absxi) {
- return math.NaN()
- }
- if scale < absxi {
- s := scale / absxi
- sumSquares = 1 + sumSquares*s*s
- scale = absxi
- } else {
- s := absxi / scale
- sumSquares += s * s
- }
- }
- if math.IsInf(scale, 1) {
- return math.Inf(1)
- }
- return scale * math.Sqrt(sumSquares)
-}
-
// L2NormInc is the level 2 norm of x.
func L2NormInc(x []float64, n, incX uintptr) (sum float64) {
var scale float64

MODIFY
internal/asm/f64/l2norm_test.go
internal/asm/f64/l2norm_test.go
package f64
-import "testing"
+import (
+ "fmt"
+ "math"
+ "testing"
+)
+
+// nanwith copied from floats package
+func nanwith(payload uint64) float64 {
+ const (
+ nanBits = 0x7ff8000000000000
+ nanMask = 0xfff8000000000000
+ )
+ return math.Float64frombits(nanBits | (payload &^ nanMask))
+}
func TestL2NormUnitary(t *testing.T) {
var src_gd float64 = 1
@@ -17,6 +30,7 @@ func TestL2NormUnitary(t *testing.T) {
{want: 3.7416573867739413, x: []float64{1, 2, 3}},
{want: 3.7416573867739413, x: []float64{-1, -2, -3}},
{want: nan, x: []float64{nan}},
+ {want: nan, x: []float64{1, inf, 3, nanwith(25), 5}},
{want: 17.88854381999832, x: []float64{8, -8, 8, -8, 8}},
{want: 2.23606797749979, x: []float64{0, 1, 0, -1, 0, 1, 0, -1, 0, 1}},
} {
@@ -87,3 +101,52 @@ func TestL2DistanceUnitary(t *testing.T) {
}
}
}
+
+func BenchmarkL2NormNetlib(b *testing.B) {
+ netlib := func(x []float64) (sum float64) {
+ var scale float64
+ sumSquares := 1.0
+ for _, v := range x {
+ if v == 0 {
+ continue
+ }
+ absxi := math.Abs(v)
+ if math.IsNaN(absxi) {
+ return math.NaN()
+ }
+ if scale < absxi {
+ s := scale / absxi
+ sumSquares = 1 + sumSquares*s*s
+ scale = absxi
+ } else {
+ s := absxi / scale
+ sumSquares += s * s
+ }
+ }
+ if math.IsInf(scale, 1) {
+ return math.Inf(1)
+ }
+ return scale * math.Sqrt(sumSquares)
+ }
+
+ tests := []struct {
+ name string
+ f func(x []float64) float64
+ }{
+ {"L2NormUnitaryNetlib", netlib},
+ {"L2NormUnitary", L2NormUnitary},
+ }
+ x[0] = randomSlice(1, 1)[0] // replace the leading zero (edge case)
+ for _, test := range tests {
+ for _, ln := range []uintptr{1, 3, 10, 30, 1e2, 3e2, 1e3, 3e3, 1e4, 3e4, 1e5} {
+ b.Run(fmt.Sprintf("%s-%d", test.name, ln), func(b *testing.B) {
+ b.SetBytes(int64(64 * ln))
+ x := x[:ln]
+ b.ResetTimer()
+ for i := 0; i < b.N; i++ {
+ test.f(x)
+ }
+ })
+ }
+ }
+}

MODIFY
internal/asm/f64/stubs_amd64.go
internal/asm/f64/stubs_amd64.go
@@ -170,3 +170,29 @@ func ScalIncTo(dst []float64, incDst uintptr, alpha float64, x []float64, n, inc
// sum += x[i]
// }
func Sum(x []float64) float64
+
+// L2NormUnitary returns the L2-norm of x.
+// var scale float64
+// sumSquares := 1.0
+// for _, v := range x {
+// if v == 0 {
+// continue
+// }
+// absxi := math.Abs(v)
+// if math.IsNaN(absxi) {
+// return math.NaN()
+// }
+// if scale < absxi {
+// s := scale / absxi
+// sumSquares = 1 + sumSquares*s*s
+// scale = absxi
+// } else {
+// s := absxi / scale
+// sumSquares += s * s
+// }
+// if math.IsInf(scale, 1) {
+// return math.Inf(1)
+// }
+// }
+// return scale * math.Sqrt(sumSquares)
+func L2NormUnitary(x []float64) (sum float64)

translate netlib l2norm algorithm to asm and remove branches

Netlib algorithm reduces overflow while calculating the l2norm of a
vector.

Translated to asm while reducing branches in NaN and Inf checks.
Overflow protection is equivalent to the Netlib standard implementation.