From 8703c0f2cf0e2441b9868b7f8845702f065df0e9 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 24 Sep 2015 17:04:13 -0400 Subject: [PATCH] Fixed an issue reported by @gassmoeller. This fixes an issue in the real_to_unit mapping code where, in the quadratic formula, it was possible to reach a division by zero issue for certain inputs. --- source/fe/mapping_q_generic.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 42fea442f7..7f20043225 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -259,8 +259,10 @@ namespace internal eta1 = -c/b; eta2 = -c/b; } - // special case #2: if c is very small: - else if (std::abs(c/b) < 1e-12) + // special case #2: if c is very small or the square root of the + // discriminant is nearly b. + else if (std::abs(c) < 1e-12*std::abs(b) + || std::abs(std::sqrt(discriminant) - b) <= 1e-14*std::abs(b)) { eta1 = (-b - std::sqrt(discriminant)) / (2*a); eta2 = (-b + std::sqrt(discriminant)) / (2*a); -- 2.39.5