From a00d5678526bfa18820d6b6840b24b1314bcc656 Mon Sep 17 00:00:00 2001 From: shaoyu Date: Mon, 9 Jun 2025 03:23:32 +0800 Subject: [PATCH 1/4] add FindBracketMono fixes https://github.com/gonum/exp/issues/70 --- root/findbracket.go | 46 ++++++++++++++++++++ root/findbracket_test.go | 93 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 root/findbracket.go create mode 100644 root/findbracket_test.go diff --git a/root/findbracket.go b/root/findbracket.go new file mode 100644 index 0000000..bd2f55c --- /dev/null +++ b/root/findbracket.go @@ -0,0 +1,46 @@ +package root + +import "math" + +// FindBracketMono finds a bracket interval [a, b] where f(a)f(b) < 0. +// f must be a monotonically increasing function. +func FindBracketMono(f func(float64) float64, guess float64) (float64, float64) { + // Make sure initial guess has the same sign as the root. + f0 := f(0) + if (guess < 0 && f0 < 0) || (guess > 0 && f0 > 0) { + guess *= -1 + } + + // r is the rate in which we adjust the interval. + var r float64 + a, fa := guess, f(guess) + if (a > 0) == (fa < 0) { + r = 2 + } else { + r = 1. / 2 + } + + // Expand bracket until x-axis is crossed. + crossed := false + b := a * r + fb := f(b) + for range 200 { + if math.Signbit(fa) != math.Signbit(fb) || fa == 0 || fb == 0 { + crossed = true + break + } + a, fa = b, fb + b *= r + fb = f(b) + } + // If unable to cross x-axis, return the largest possible bracket. + if !crossed { + if r > 1 { + return a, math.Inf(int(math.Copysign(1, b))) + } else { + return a, 0 + } + } + + return a, b +} diff --git a/root/findbracket_test.go b/root/findbracket_test.go new file mode 100644 index 0000000..7bbd90f --- /dev/null +++ b/root/findbracket_test.go @@ -0,0 +1,93 @@ +package root_test + +import ( + "fmt" + "math" + "testing" + + "gonum.org/v1/exp/root" +) + +var findBracketMonoTests = []struct { + name string + f func(float64) float64 + guess float64 +}{ + // Based on https://github.com/boostorg/math/blob/boost-1.88.0/test/test_toms748_solve.cpp + {name: "f4.4", f: func(x float64) float64 { return math.Pow(x, 4) - 0.2 }, guess: 2}, + {name: "f4.6", f: func(x float64) float64 { return math.Pow(x, 6) - 0.2 }, guess: 2}, + {name: "f4.8", f: func(x float64) float64 { return math.Pow(x, 8) - 0.2 }, guess: 2}, + {name: "f4.10", f: func(x float64) float64 { return math.Pow(x, 10) - 0.2 }, guess: 2}, + {name: "f4.12", f: func(x float64) float64 { return math.Pow(x, 12) - 0.2 }, guess: 2}, + + // Based on https://github.com/boostorg/math/blob/boost-1.88.0/test/test_root_finding_concepts.cpp + {name: "f1", f: func(x float64) float64 { return x*x*x - 27 }, guess: 27}, +} + +func TestFindBracketMono(t *testing.T) { + t.Parallel() + + for _, test := range findBracketMonoTests { + t.Run(fmt.Sprintf("%s", test.name), func(t *testing.T) { + a, b := root.FindBracketMono(test.f, test.guess) + fa, fb := test.f(a), test.f(b) + if fa*fb > 0 { + t.Fatalf("%s: invalid bracket (%f, %f)", test.name, fa, fb) + } + }) + } +} + +func TestFindBracketMonoSpecialCases(t *testing.T) { + // Bracket to positive infinity. + f := func(x float64) float64 { return math.Atan(x) - 2 } + var guess float64 = 3 + a, b := root.FindBracketMono(f, guess) + if !(a > 0 && math.IsInf(b, 1)) { + t.Errorf("FindBracketMono(Atan-2, %f): got (%f, %f) want (+a, +Inf)", guess, a, b) + } + + // Bracket to negative infinity. + f = func(x float64) float64 { return math.Atan(x) + 2 } + guess = 3 + a, b = root.FindBracketMono(f, guess) + if !(a < 0 && math.IsInf(b, -1)) { + t.Errorf("FindBracketMono(Atan+2, %f): got (%f, %f) want (-a, -Inf)", guess, a, b) + } + + // Root is a tiny positive value. + f = func(x float64) float64 { + rt := math.SmallestNonzeroFloat64 + switch { + case x > rt: + return 1 + case x < rt: + return -1 + default: + return 0 + } + } + guess = 3 + a, b = root.FindBracketMono(f, guess) + if !(a > 0 && b == 0) { + t.Errorf("FindBracketMono(f, %f): got (%g, %g) want (+a, 0)", guess, a, b) + } + + // Root is a tiny negative value. + f = func(x float64) float64 { + rt := -math.SmallestNonzeroFloat64 + switch { + case x > rt: + return 1 + case x < rt: + return -1 + default: + return 0 + } + } + guess = -3 + a, b = root.FindBracketMono(f, guess) + if !(a < 0 && b == 0) { + t.Errorf("FindBracketMono(f, %f): got (%g, %g) want (-a, 0)", guess, a, b) + } +} From 8e9407a21abfc3bec59dcd3d3dfd124c07b41f8b Mon Sep 17 00:00:00 2001 From: shaoyu Date: Mon, 9 Jun 2025 18:51:38 +0800 Subject: [PATCH 2/4] comments for maxiter --- root/findbracket.go | 18 +++++++++++------- root/findbracket_test.go | 4 ++++ 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/root/findbracket.go b/root/findbracket.go index bd2f55c..22fac59 100644 --- a/root/findbracket.go +++ b/root/findbracket.go @@ -1,10 +1,14 @@ +// Copyright ©2025 The Gonum 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 root import "math" // FindBracketMono finds a bracket interval [a, b] where f(a)f(b) < 0. // f must be a monotonically increasing function. -func FindBracketMono(f func(float64) float64, guess float64) (float64, float64) { +func FindBracketMono(f func(float64) float64, guess float64) (a, b float64) { // Make sure initial guess has the same sign as the root. f0 := f(0) if (guess < 0 && f0 < 0) || (guess > 0 && f0 > 0) { @@ -17,14 +21,16 @@ func FindBracketMono(f func(float64) float64, guess float64) (float64, float64) if (a > 0) == (fa < 0) { r = 2 } else { - r = 1. / 2 + r = 0.5 } // Expand bracket until x-axis is crossed. + // maxiter value is based on https://github.com/boostorg/math/blob/boost-1.88.0/include/boost/math/policies/policy.hpp#L130 + const maxiter = 200 crossed := false - b := a * r + b = a * r fb := f(b) - for range 200 { + for range maxiter { if math.Signbit(fa) != math.Signbit(fb) || fa == 0 || fb == 0 { crossed = true break @@ -37,10 +43,8 @@ func FindBracketMono(f func(float64) float64, guess float64) (float64, float64) if !crossed { if r > 1 { return a, math.Inf(int(math.Copysign(1, b))) - } else { - return a, 0 } + return a, 0 } - return a, b } diff --git a/root/findbracket_test.go b/root/findbracket_test.go index 7bbd90f..5355ba0 100644 --- a/root/findbracket_test.go +++ b/root/findbracket_test.go @@ -1,3 +1,7 @@ +// Copyright ©2025 The Gonum 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 root_test import ( From 862765afc625e7040f8dd2c4e3bd1e657c556482 Mon Sep 17 00:00:00 2001 From: shaoyu Date: Tue, 10 Jun 2025 13:07:15 +0800 Subject: [PATCH 3/4] change multiple assignments to single ones --- root/findbracket.go | 3 ++- root/findbracket_test.go | 55 ++++++++++++++++------------------------ 2 files changed, 24 insertions(+), 34 deletions(-) diff --git a/root/findbracket.go b/root/findbracket.go index 22fac59..b961368 100644 --- a/root/findbracket.go +++ b/root/findbracket.go @@ -17,7 +17,8 @@ func FindBracketMono(f func(float64) float64, guess float64) (a, b float64) { // r is the rate in which we adjust the interval. var r float64 - a, fa := guess, f(guess) + a = guess + fa := f(a) if (a > 0) == (fa < 0) { r = 2 } else { diff --git a/root/findbracket_test.go b/root/findbracket_test.go index 5355ba0..1c0c935 100644 --- a/root/findbracket_test.go +++ b/root/findbracket_test.go @@ -42,25 +42,15 @@ func TestFindBracketMono(t *testing.T) { } } -func TestFindBracketMonoSpecialCases(t *testing.T) { - // Bracket to positive infinity. - f := func(x float64) float64 { return math.Atan(x) - 2 } - var guess float64 = 3 - a, b := root.FindBracketMono(f, guess) - if !(a > 0 && math.IsInf(b, 1)) { - t.Errorf("FindBracketMono(Atan-2, %f): got (%f, %f) want (+a, +Inf)", guess, a, b) - } - - // Bracket to negative infinity. - f = func(x float64) float64 { return math.Atan(x) + 2 } - guess = 3 - a, b = root.FindBracketMono(f, guess) - if !(a < 0 && math.IsInf(b, -1)) { - t.Errorf("FindBracketMono(Atan+2, %f): got (%f, %f) want (-a, -Inf)", guess, a, b) - } - - // Root is a tiny positive value. - f = func(x float64) float64 { +var findBracketMonoSpecials = []struct { + name string + f func(float64) float64 + guess float64 + criteria func(a, b float64) bool +}{ + {name: "+Inf", f: func(x float64) float64 { return math.Atan(x) - 2 }, guess: 3, criteria: func(a, b float64) bool { return a > 0 && math.IsInf(b, 1) }}, + {name: "-Inf", f: func(x float64) float64 { return math.Atan(x) + 2 }, guess: 3, criteria: func(a, b float64) bool { return a < 0 && math.IsInf(b, -1) }}, + {name: "tiny positive", f: func(x float64) float64 { rt := math.SmallestNonzeroFloat64 switch { case x > rt: @@ -70,15 +60,8 @@ func TestFindBracketMonoSpecialCases(t *testing.T) { default: return 0 } - } - guess = 3 - a, b = root.FindBracketMono(f, guess) - if !(a > 0 && b == 0) { - t.Errorf("FindBracketMono(f, %f): got (%g, %g) want (+a, 0)", guess, a, b) - } - - // Root is a tiny negative value. - f = func(x float64) float64 { + }, guess: 3, criteria: func(a, b float64) bool { return a > 0 && b == 0 }}, + {name: "tiny negative", f: func(x float64) float64 { rt := -math.SmallestNonzeroFloat64 switch { case x > rt: @@ -88,10 +71,16 @@ func TestFindBracketMonoSpecialCases(t *testing.T) { default: return 0 } - } - guess = -3 - a, b = root.FindBracketMono(f, guess) - if !(a < 0 && b == 0) { - t.Errorf("FindBracketMono(f, %f): got (%g, %g) want (-a, 0)", guess, a, b) + }, guess: -3, criteria: func(a, b float64) bool { return a < 0 && b == 0 }}, +} + +func TestFindBracketMonoSpecialCases(t *testing.T) { + for _, tc := range findBracketMonoSpecials { + t.Run(tc.name, func(t *testing.T) { + a, b := root.FindBracketMono(tc.f, tc.guess) + if !tc.criteria(a, b) { + t.Fatalf("%s: wrong bracket (%g, %g)", tc.name, a, b) + } + }) } } From d21c41005e3efd3260fc2eb894f6aaa262a9728e Mon Sep 17 00:00:00 2001 From: shaoyu Date: Wed, 11 Jun 2025 19:17:45 +0800 Subject: [PATCH 4/4] ensure returned bracket satisfies a <= b --- root/findbracket.go | 12 +++++++-- root/findbracket_test.go | 57 ++++++++++++++++------------------------ 2 files changed, 32 insertions(+), 37 deletions(-) diff --git a/root/findbracket.go b/root/findbracket.go index b961368..5ba960b 100644 --- a/root/findbracket.go +++ b/root/findbracket.go @@ -8,6 +8,7 @@ import "math" // FindBracketMono finds a bracket interval [a, b] where f(a)f(b) < 0. // f must be a monotonically increasing function. +// guess is the initial guess of the bracket search. func FindBracketMono(f func(float64) float64, guess float64) (a, b float64) { // Make sure initial guess has the same sign as the root. f0 := f(0) @@ -43,9 +44,16 @@ func FindBracketMono(f func(float64) float64, guess float64) (a, b float64) { // If unable to cross x-axis, return the largest possible bracket. if !crossed { if r > 1 { - return a, math.Inf(int(math.Copysign(1, b))) + b = math.Inf(int(math.Copysign(1, b))) + } else { + b = 0 } - return a, 0 } + + // Ensure a <= b + if a > b { + a, b = b, a + } + return a, b } diff --git a/root/findbracket_test.go b/root/findbracket_test.go index 1c0c935..ff4f1e1 100644 --- a/root/findbracket_test.go +++ b/root/findbracket_test.go @@ -5,7 +5,6 @@ package root_test import ( - "fmt" "math" "testing" @@ -13,9 +12,10 @@ import ( ) var findBracketMonoTests = []struct { - name string - f func(float64) float64 - guess float64 + name string + f func(float64) float64 + guess float64 + validBracket func(a, b float64) bool }{ // Based on https://github.com/boostorg/math/blob/boost-1.88.0/test/test_toms748_solve.cpp {name: "f4.4", f: func(x float64) float64 { return math.Pow(x, 4) - 0.2 }, guess: 2}, @@ -26,30 +26,10 @@ var findBracketMonoTests = []struct { // Based on https://github.com/boostorg/math/blob/boost-1.88.0/test/test_root_finding_concepts.cpp {name: "f1", f: func(x float64) float64 { return x*x*x - 27 }, guess: 27}, -} -func TestFindBracketMono(t *testing.T) { - t.Parallel() - - for _, test := range findBracketMonoTests { - t.Run(fmt.Sprintf("%s", test.name), func(t *testing.T) { - a, b := root.FindBracketMono(test.f, test.guess) - fa, fb := test.f(a), test.f(b) - if fa*fb > 0 { - t.Fatalf("%s: invalid bracket (%f, %f)", test.name, fa, fb) - } - }) - } -} - -var findBracketMonoSpecials = []struct { - name string - f func(float64) float64 - guess float64 - criteria func(a, b float64) bool -}{ - {name: "+Inf", f: func(x float64) float64 { return math.Atan(x) - 2 }, guess: 3, criteria: func(a, b float64) bool { return a > 0 && math.IsInf(b, 1) }}, - {name: "-Inf", f: func(x float64) float64 { return math.Atan(x) + 2 }, guess: 3, criteria: func(a, b float64) bool { return a < 0 && math.IsInf(b, -1) }}, + // Special cases. + {name: "+Inf", f: func(x float64) float64 { return math.Atan(x) - 2 }, guess: 3, validBracket: func(a, b float64) bool { return a > 0 && math.IsInf(b, 1) }}, + {name: "-Inf", f: func(x float64) float64 { return math.Atan(x) + 2 }, guess: 3, validBracket: func(a, b float64) bool { return math.IsInf(a, -1) && b < 0 }}, {name: "tiny positive", f: func(x float64) float64 { rt := math.SmallestNonzeroFloat64 switch { @@ -60,7 +40,7 @@ var findBracketMonoSpecials = []struct { default: return 0 } - }, guess: 3, criteria: func(a, b float64) bool { return a > 0 && b == 0 }}, + }, guess: 3, validBracket: func(a, b float64) bool { return a == 0 && b > 0 }}, {name: "tiny negative", f: func(x float64) float64 { rt := -math.SmallestNonzeroFloat64 switch { @@ -71,15 +51,22 @@ var findBracketMonoSpecials = []struct { default: return 0 } - }, guess: -3, criteria: func(a, b float64) bool { return a < 0 && b == 0 }}, + }, guess: -3, validBracket: func(a, b float64) bool { return a < 0 && b == 0 }}, } -func TestFindBracketMonoSpecialCases(t *testing.T) { - for _, tc := range findBracketMonoSpecials { - t.Run(tc.name, func(t *testing.T) { - a, b := root.FindBracketMono(tc.f, tc.guess) - if !tc.criteria(a, b) { - t.Fatalf("%s: wrong bracket (%g, %g)", tc.name, a, b) +func TestFindBracketMono(t *testing.T) { + t.Parallel() + + for _, test := range findBracketMonoTests { + t.Run(test.name, func(t *testing.T) { + validBracket := test.validBracket + if validBracket == nil { + validBracket = func(a, b float64) bool { return test.f(a)*test.f(b) < 0 && a <= b } + } + + a, b := root.FindBracketMono(test.f, test.guess) + if !validBracket(a, b) { + t.Errorf("%s: invalid bracket (%f, %f)", test.name, a, b) } }) }