aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorCharles Karney <charles.karney@sri.com>2019-09-18 07:08:55 -0400
committerCharles Karney <charles.karney@sri.com>2019-09-18 07:08:55 -0400
commit93d0f53db70b581718f29601df11f4b819ef943e (patch)
treeb5dc2c1f3af4b2a8350636b75f724eaa95165224 /src
parentc080801637b5b3c1fe7462829411a6859ce581dd (diff)
downloadPROJ-93d0f53db70b581718f29601df11f4b819ef943e.tar.gz
PROJ-93d0f53db70b581718f29601df11f4b819ef943e.zip
Try removing C99 math functions from math.cpp (but leave isnan).
Diffstat (limited to 'src')
-rw-r--r--src/math.cpp138
-rw-r--r--src/proj_math.h21
2 files changed, 0 insertions, 159 deletions
diff --git a/src/math.cpp b/src/math.cpp
index 6c1b6d1b..90d35001 100644
--- a/src/math.cpp
+++ b/src/math.cpp
@@ -41,141 +41,3 @@ int pj_isnan (double x) {
/* cppcheck-suppress duplicateExpression */
return x != x;
}
-
-#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
-
-/* Define C99 compatible versions of
- * hypot
- * log1p
- * asinh
- * atanh
- * copysign
- * cbrt
- * remainder
- * remquo
- * round
- * lround
- */
-
-/* Compute hypotenuse */
-double pj_hypot(double x, double y) {
- x = fabs(x);
- y = fabs(y);
- if (x < y) {
- x /= y; /* y is nonzero */
- return y * sqrt(1 + x * x);
- } else {
- y /= (x != 0 ? x : 1);
- return x * sqrt(1 + y * y);
- }
-}
-
-/* Compute log(1+x) accurately */
-double pj_log1p(double x) {
- volatile double
- y = 1 + x,
- z = y - 1;
- /* Here's the explanation for this magic: y = 1 + z, exactly, and z
- * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
- * a good approximation to the true log(1 + x)/x. The multiplication x *
- * (log(y)/z) introduces little additional error. */
- return z == 0 ? x : x * log(y) / z;
-}
-
-/* Compute asinh(x) accurately */
-double pj_asinh(double x) {
- double y = fabs(x); /* Enforce odd parity */
- y = pj_log1p(y * (1 + y/(pj_hypot(1.0, y) + 1)));
- return x > 0 ? y : (x < 0 ? -y : x); /* asinh(-0.0) = -0.0 */
-}
-
-/* Compute atanh(x) accurately */
-double pj_atanh(double x) {
- double y = fabs(x); /* Enforce odd parity */
- y = pj_log1p(2 * y/(1 - y))/2;
- return x > 0 ? y : (x < 0 ? -y : x); /* atanh(-0.0) = -0.0 */
-}
-
-/* Implement copysign(x, y) */
-double pj_copysign(double x, double y) {
- /* 1/y trick to get the sign of -0.0 */
- return fabs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
-}
-
-/* Implement cbrt(x) */
-double pj_cbrt(double x) {
- double y = pow(fabs(x), 1/3.0); /* Return the real cube root */
- return x > 0 ? y : (x < 0 ? -y : x); /* cbrt(-0.0) = -0.0 */
-}
-
-/* Implement remainder(x, y) with ties to even */
-double pj_remainder(double x, double y) {
- double z;
- y = fabs(y); /* The result doesn't depend on the sign of y */
- z = fmod(x, y);
- if (z == 0)
- /* This shouldn't be necessary. However, before version 14 (2015),
- * Visual Studio had problems dealing with -0.0. Specifically
- * VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
- * python 2.7 on Windows 32-bit machines has the same problem. */
- z = pj_copysign(z, x);
- else if (2 * fabs(z) == y)
- z -= fmod(x, 2 * y) - z; /* Implement ties to even */
- else if (2 * fabs(z) > y)
- z += (z < 0 ? y : -y); /* Fold remaining cases to (-y/2, y/2) */
- return z;
-}
-
-/* Implement remquo(x, y, n) with n giving low 3 bits + sign of x/y */
-double pj_remquo(double x, double y, int* n) {
- double z = pj_remainder(x, y);
- if (n) {
- double
- a = pj_remainder(x, 2 * y),
- b = pj_remainder(x, 4 * y),
- c = pj_remainder(x, 8 * y);
- *n = (a > z ? 1 : (a < z ? -1 : 0));
- *n += (b > a ? 2 : (b < a ? -2 : 0));
- *n += (c > b ? 4 : (c < b ? -4 : 0));
- if (y < 0) *n *= -1;
- if (y != 0) {
- if (x/y > 0 && *n <= 0)
- *n += 8;
- else if (x/y < 0 && *n >= 0)
- *n -= 8;
- }
- }
- return z;
-}
-
-/* Implement round(x) */
-double pj_round(double x) {
- /* The handling of corner cases is copied from boost; see
- * https://github.com/boostorg/math/pull/8
- * with improvements to return -0.0 when appropriate */
- double t;
- if (0 < x && x < 0.5)
- return +0.0;
- else if (0 > x && x > -0.5)
- return -0.0;
- else if (x > 0) {
- t = ceil(x);
- return 0.5 < t - x ? t - 1 : t;
- } else if (x < 0) {
- t = floor(x);
- return 0.5 < x - t ? t + 1 : t;
- } else /* +/-0 and NaN */
- return x; /* retain sign of 0 */
-}
-
-/* Implement lround(x) */
-long pj_lround(double x) {
- /* Default value for overflow + NaN + (x == LONG_MIN) */
- long r = LONG_MIN;
- x = pj_round(x);
- if (fabs(x) < -(double)r) /* Assume (double)LONG_MIN is exact */
- r = (long)x;
- return r;
-}
-
-#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
diff --git a/src/proj_math.h b/src/proj_math.h
index 414df805..698654dd 100644
--- a/src/proj_math.h
+++ b/src/proj_math.h
@@ -61,29 +61,8 @@ extern "C" {
#endif
#endif
-double pj_hypot(double x, double y);
-double pj_log1p(double x);
-double pj_asinh(double x);
-double pj_atanh(double x);
-double pj_copysign(double x, double y);
-double pj_cbrt(double x);
-double pj_remainder(double x, double y);
-double pj_remquo(double x, double y, int* n);
-double pj_round(double x);
-long pj_lround(double x);
int PROJ_DLL pj_isnan(double x);
-#define hypot pj_hypot
-#define log1p pj_log1p
-#define asinh pj_asinh
-#define atanh pj_atanh
-#define copysign pj_copysign
-#define cbrt pj_cbrt
-#define remainder pj_remainder
-#define remquo pj_remquo
-#define round pj_round
-#define lround pj_lround
-
#ifdef isnan
#undef isnan
#endif