UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-19 00:37:25 +01:00
Code Issues Packages Projects Releases Wiki Activity

Fix potential overflow in complex radius calculation.

Browse Source
This commit is contained in:
Jacob Young
2017-09-09 10:06:02 -04:00
committed by EmilieNumworks
parent 851c145d31
commit 7e2cc375d2
3 changed files with 42 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
liba/include/math.h
View File

@@ -10,6 +10,8 @@ LIBA_BEGIN_DECLS
#define INFINITY __builtin_inff()
#define M_E 2.71828182845904523536028747135266250
#define M_PI 3.14159265358979323846264338327950288
#define M_PI_4 0.78539816339744830961566084581987572
#define M_SQRT2 1.41421356237309504880168872420969808
#define FP_INFINITE 0x01
#define FP_NAN 0x02

32
poincare/src/complex.cpp
View File

@@ -161,10 +161,38 @@ T Complex<T>::b() const {
template <class T>
T Complex<T>::r() const {
// We want to avoid a^2 and b^2 which could both easily overflow.
// min, max = minmax(abs(a), abs(b)) (*minmax returns both arguments sorted*)
// abs(a + bi) == sqrt(a^2 + b^2)
// == sqrt(abs(a)^2 + abs(b)^2)
// == sqrt(min^2 + max^2)
// == sqrt((min^2 + max^2) * max^2/max^2)
// == sqrt((min^2 + max^2) / max^2)*sqrt(max^2)
// == sqrt(min^2/max^2 + 1) * max
// == sqrt((min/max)^2 + 1) * max
// min >= 0 &&
// max >= 0 &&
// min <= max => min/max <= 1
// => (min/max)^2 <= 1
// => (min/max)^2 + 1 <= 2
// => sqrt((min/max)^2 + 1) <= sqrt(2)
// So the calculation is guaranteed to not overflow until the final multiply.
// If (min/max)^2 underflows then min doesn't contribute anything significant
// compared to max, and the formula reduces to simply max as it should.
// We do need to be careful about the case where a == 0 && b == 0 which would
// cause a division by zero.
T min = std::fabs(m_a);
if (m_b == 0) {
return std::fabs(m_a);
return min;
}
return std::sqrt(m_a*m_a + m_b*m_b);
T max = std::fabs(m_b);
if (max < min) {
T temp = min;
min = max;
max = temp;
}
T temp = min/max;
return std::sqrt(temp*temp + 1) * max;
}
template <class T>

10
poincare/test/complex.cpp
View File

@@ -128,4 +128,14 @@ QUIZ_CASE(poincare_complex_constructor) {
b = new Complex<double>(Complex<double>::Polar(12.04159457879229548012824103, 1.4876550949));
assert(std::fabs(b->a() - 1.0) < 0.0000000001 && std::fabs(b->b()-12.0) < 0.0000000001);
delete b;
Complex<float> * c = new Complex<float>(Complex<float>::Cartesian(-2.0e20f, 2.0e20f));
assert(c->a() == -2.0e20f && c->b() == 2.0e20f);
assert(c->r() == 2.0e20f*(float)M_SQRT2 && c->th() == 3*(float)M_PI_4);
delete c;
Complex<double> * d = new Complex<double>(Complex<double>::Cartesian(1.0e155, -1.0e155));
assert(d->a() == 1.0e155 && d->b() == -1.0e155);
assert(d->r() == 1.0e155*M_SQRT2 && d->th() == -M_PI_4);
delete d;
}