diff --git a/src/genz-malik.jl b/src/genz-malik.jl index f2ca0c5..252a6df 100644 --- a/src/genz-malik.jl +++ b/src/genz-malik.jl @@ -114,9 +114,9 @@ error `E` (via the given `norm`), and the suggested coordinate `k` ∈ `1:n` to subdivide next. """ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) where {n,T} - c = (a+b)*T(0.5) - Δ = (b-a)*T(0.5) - V = abs(prod(Δ)) + c = T(0.5).*(a.+b) + Δ = T(0.5).*abs.(b.-a) + V = prod(Δ) f₁ = f(c) @@ -124,6 +124,7 @@ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) wher f₃ = zero(f₁) kdivide = 1 maxdivdiff = zero(norm(f₁)) + ε = eps(typeof(maxdivdiff)) twelvef₁ = 12f₁ for i = 1:n p₂ = Δ .* g.p[1][i] @@ -135,9 +136,11 @@ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) wher # fourth divided difference: f₃ᵢ-2f₁ - 7*(f₂ᵢ-2f₁), # where 7 = (λ₃/λ₂)^2 [see van Dooren and de Ridder] divdiff = norm(f₃ᵢ + twelvef₁ - 7*f₂ᵢ) - if divdiff > maxdivdiff + if (δ = divdiff - maxdivdiff) > ε kdivide = i maxdivdiff = divdiff + elseif abs(δ) < ε && Δ[i] > Δ[kdivide] + kdivide = i end end