avoid catastrophic cancellation in divided differences (#858)

* avoid catastrophic cancellation in divided differences

* rename rteps thresh

Author: araujoms

src/Cones/epitrrelentropytri.jl

@@ -515,7 +515,7 @@ function dder3(
 end
 
 function Δ2!(Δ2::Matrix{T}, λ::Vector{T}, log_λ::Vector{T}) where {T <: Real}
-    rteps = sqrt(eps(T))
+    thresh = cbrt(eps(T))
     d = length(λ)
 
     @inbounds for j in 1:d
@@ -524,7 +524,7 @@ function Δ2!(Δ2::Matrix{T}, λ::Vector{T}, log_λ::Vector{T}) where {T <: Real
         for i in 1:(j - 1)
             λ_i = λ[i]
             λ_ij = λ_i - λ_j
-            if abs(λ_ij) < rteps
+            if abs(λ_ij) < thresh
                 Δ2[i, j] = 2 / (λ_i + λ_j)
             else
                 Δ2[i, j] = (log_λ[i] - lλ_j) / λ_ij
@@ -540,17 +540,17 @@ end
 
 function Δ3!(Δ3::Array{T, 3}, Δ2::Matrix{T}, λ::Vector{T}) where {T <: Real}
     @assert issymmetric(Δ2) # must be symmetric (wrapper is less efficient)
-    rteps = sqrt(eps(T))
+    thresh = eps(T)^(2/9)
     d = length(λ)
 
     @inbounds for k in 1:d, j in 1:k, i in 1:j
         λ_j = λ[j]
         λ_k = λ[k]
         λ_jk = λ_j - λ_k
-        if abs(λ_jk) < rteps
+        if abs(λ_jk) < thresh
             λ_i = λ[i]
             λ_ij = λ_i - λ_j
-            if abs(λ_ij) < rteps
+            if abs(λ_ij) < thresh
                 t = abs2(3 / (λ_i + λ_j + λ_k)) / -2
             else
                 t = (Δ2[i, j] - Δ2[j, k]) / λ_ij
@@ -661,7 +661,7 @@ function Δ4_ij!(
     Δ3::Array{T, 3},
     λ::Vector{T},
 ) where {T <: Real}
-    rteps = sqrt(eps(T))
+    thresh = eps(T)^(4/27)
     d = length(λ)
     λ_i = λ[i]
     λ_j = λ[j]
@@ -672,11 +672,12 @@ function Δ4_ij!(
         λ_ij = λ_i - λ_j
         λ_ik = λ_i - λ_k
         λ_il = λ_i - λ_l
-        B_ik = (abs(λ_ik) < rteps)
-        B_il = (abs(λ_il) < rteps)
+        B_ik = (abs(λ_ik) < thresh)
+        B_il = (abs(λ_il) < thresh)
 
-        if (abs(λ_ij) < rteps) && B_ik && B_il
-            t = λ_i^-3 / 3
+        if (abs(λ_ij) < thresh) && B_ik && B_il
+            λ_mean = (λ_i + λ_j + λ_k + λ_l)/4
+            t = λ_mean^-3 / 3
         elseif B_ik && B_il
             t = (Δ3[i, i, i] - Δ3[i, i, j]) / λ_ij
         elseif B_il