1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1996-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if ! defined (octave_lo_mappers_h) 27 #define octave_lo_mappers_h 1 28 29 #include "octave-config.h" 30 31 #include <cmath> 32 33 #include <limits> 34 35 #include "lo-ieee.h" 36 #include "oct-cmplx.h" 37 #include "oct-inttypes-fwd.h" 38 39 namespace octave 40 { 41 namespace math 42 { 43 extern OCTAVE_API bool isna (double x); 44 extern OCTAVE_API bool isna (float x); 45 extern OCTAVE_API bool isna (const Complex& x); 46 extern OCTAVE_API bool isna (const FloatComplex& x); 47 48 extern OCTAVE_API bool is_NaN_or_NA (const Complex& x); 49 extern OCTAVE_API bool is_NaN_or_NA (const FloatComplex& x); 50 copysign(double x,double y)51 inline double copysign (double x, double y) { return std::copysign (x, y); } copysign(float x,float y)52 inline float copysign (float x, float y) { return std::copysignf (x, y); } 53 signbit(double x)54 inline double signbit (double x) { return std::signbit (x); } signbit(float x)55 inline float signbit (float x) { return std::signbit (x); } 56 57 // Test for negative sign. 58 extern OCTAVE_API bool negative_sign (double x); 59 extern OCTAVE_API bool negative_sign (float x); 60 61 // Test for positive sign. positive_sign(double x)62 inline bool positive_sign (double x) { return ! negative_sign (x); } positive_sign(float x)63 inline bool positive_sign (float x) { return ! negative_sign (x); } 64 65 extern OCTAVE_API Complex acos (const Complex& x); 66 extern OCTAVE_API FloatComplex acos (const FloatComplex& x); 67 68 extern OCTAVE_API Complex asin (const Complex& x); 69 extern OCTAVE_API FloatComplex asin (const FloatComplex& x); 70 atan(const Complex & x)71 inline Complex atan (const Complex& x) { return std::atan (x); } atan(const FloatComplex & x)72 inline FloatComplex atan (const FloatComplex& x) { return std::atan (x); } 73 74 // The C++ standard would normally return a std::complex value for conj 75 // even when the input is fully real. Octave overrides this. conj(double x)76 inline double conj (double x) { return x; } conj(float x)77 inline float conj (float x) { return x; } 78 79 template <typename T> 80 std::complex<T> conj(const std::complex<T> & x)81 conj (const std::complex<T>& x) 82 { 83 return std::conj (x); 84 } 85 log2(double x)86 inline double log2 (double x) { return std::log2 (x); } log2(float x)87 inline float log2 (float x) { return std::log2f (x); } 88 89 extern OCTAVE_API Complex log2 (const Complex& x); 90 extern OCTAVE_API FloatComplex log2 (const FloatComplex& x); 91 92 extern OCTAVE_API double log2 (double x, int& exp); 93 extern OCTAVE_API float log2 (float x, int& exp); 94 95 extern OCTAVE_API Complex log2 (const Complex& x, int& exp); 96 extern OCTAVE_API FloatComplex log2 (const FloatComplex& x, int& exp); 97 exp2(double x)98 inline double exp2 (double x) { return std::exp2 (x); } exp2(float x)99 inline float exp2 (float x) { return std::exp2f (x); } 100 101 template <typename T> 102 std::complex<T> ceil(const std::complex<T> & x)103 ceil (const std::complex<T>& x) 104 { 105 return std::complex<T> (std::ceil (std::real (x)), 106 std::ceil (std::imag (x))); 107 } 108 109 template <typename T> 110 std::complex<T> trunc(const std::complex<T> & x)111 trunc (const std::complex<T>& x) 112 { 113 return std::complex<T> (std::trunc (std::real (x)), 114 std::trunc (std::imag (x))); 115 } 116 117 // Provide alias for trunc under the more familiar name of fix. fix(double x)118 inline double fix (double x) { return std::trunc (x); } fix(float x)119 inline float fix (float x) { return std::trunc (x); } 120 121 template <typename T> 122 std::complex<T> fix(const std::complex<T> & x)123 fix (const std::complex<T>& x) 124 { 125 return trunc (x); 126 } 127 128 template <typename T> 129 std::complex<T> floor(const std::complex<T> & x)130 floor (const std::complex<T>& x) 131 { 132 return std::complex<T> (std::floor (std::real (x)), 133 std::floor (std::imag (x))); 134 } 135 round(double x)136 inline double round (double x) { return std::round (x); } round(float x)137 inline float round (float x) { return std::roundf (x); } 138 139 template <typename T> 140 std::complex<T> round(const std::complex<T> & x)141 round (const std::complex<T>& x) 142 { 143 return std::complex<T> (round (std::real (x)), round (std::imag (x))); 144 } 145 146 inline double roundb(double x)147 roundb (double x) 148 { 149 double t = round (x); 150 151 if (fabs (x - t) == 0.5) 152 t = 2 * std::trunc (0.5 * t); 153 154 return t; 155 } 156 157 inline float roundb(float x)158 roundb (float x) 159 { 160 float t = round (x); 161 162 if (fabsf (x - t) == 0.5f) 163 t = 2 * std::trunc (0.5f * t); 164 165 return t; 166 } 167 168 template <typename T> 169 std::complex<T> roundb(const std::complex<T> & x)170 roundb (const std::complex<T>& x) 171 { 172 return std::complex<T> (roundb (std::real (x)), roundb (std::imag (x))); 173 } 174 175 extern OCTAVE_API double frexp (double x, int *expptr); 176 extern OCTAVE_API float frexp (float x, int *expptr); 177 isnan(bool)178 inline bool isnan (bool) { return false; } isnan(char)179 inline bool isnan (char) { return false; } 180 isnan(double x)181 inline bool isnan (double x) { return std::isnan (x); } isnan(float x)182 inline bool isnan (float x) { return std::isnan (x); } 183 184 // FIXME: Do we need the isnan overload for complex? 185 template <typename T> 186 bool isnan(const std::complex<T> & x)187 isnan (const std::complex<T>& x) 188 { 189 return (isnan (std::real (x)) || isnan (std::imag (x))); 190 } 191 isfinite(double x)192 inline bool isfinite (double x) { return std::isfinite (x); } isfinite(float x)193 inline bool isfinite (float x) { return std::isfinite (x); } 194 195 // FIXME: Do we need isfinite overload for complex? 196 template <typename T> 197 bool isfinite(const std::complex<T> & x)198 isfinite (const std::complex<T>& x) 199 { 200 return (isfinite (std::real (x)) && isfinite (std::imag (x))); 201 } 202 isinf(double x)203 inline bool isinf (double x) { return std::isinf (x); } isinf(float x)204 inline bool isinf (float x) { return std::isinf (x); } 205 206 // FIXME: Do we need isinf overload for complex? 207 template <typename T> 208 bool isinf(const std::complex<T> & x)209 isinf (const std::complex<T>& x) 210 { 211 return (isinf (std::real (x)) || isinf (std::imag (x))); 212 } 213 214 // Some useful tests, that are commonly repeated. 215 // Test for a finite integer. 216 217 // FIXME: Benchmark whether trunc might be faster than round here. isinteger(double x)218 inline bool isinteger (double x) { return isfinite (x) && x == round (x); } isinteger(float x)219 inline bool isinteger (float x) { return isfinite (x) && x == round (x); } 220 221 inline double signum(double x)222 signum (double x) 223 { 224 double tmp = 0.0; 225 226 if (x < 0.0) 227 tmp = -1.0; 228 else if (x > 0.0) 229 tmp = 1.0; 230 231 return isnan (x) ? numeric_limits<double>::NaN () : tmp; 232 } 233 234 inline float signum(float x)235 signum (float x) 236 { 237 float tmp = 0.0f; 238 239 if (x < 0.0f) 240 tmp = -1.0f; 241 else if (x > 0.0f) 242 tmp = 1.0f; 243 244 return isnan (x) ? numeric_limits<float>::NaN () : tmp; 245 } 246 247 template <typename T> 248 std::complex<T> signum(const std::complex<T> & x)249 signum (const std::complex<T>& x) 250 { 251 T tmp = abs (x); 252 253 return tmp == 0 ? 0.0 : x / tmp; 254 } 255 256 // Convert X to the nearest integer value. Should not pass NaN to 257 // this function. 258 259 // For integer types? Hmm. Need to be sure T is an integer type... 260 template <typename T> 261 T x_nint(T x)262 x_nint (T x) 263 { 264 return x; 265 } 266 267 template <> x_nint(double x)268 inline double x_nint (double x) 269 { 270 return (isfinite (x) ? std::floor (x + 0.5) : x); 271 } 272 273 template <> x_nint(float x)274 inline float x_nint (float x) 275 { 276 return (isfinite (x) ? std::floor (x + 0.5f) : x); 277 } 278 279 extern OCTAVE_API octave_idx_type nint_big (double x); 280 extern OCTAVE_API octave_idx_type nint_big (float x); 281 282 extern OCTAVE_API int nint (double x); 283 extern OCTAVE_API int nint (float x); 284 285 template <typename T> 286 T mod(T x,T y)287 mod (T x, T y) 288 { 289 T retval; 290 291 if (y == 0) 292 retval = x; 293 else 294 { 295 T q = x / y; 296 297 if (x_nint (y) != y 298 && (std::abs ((q - x_nint (q)) / x_nint (q)) 299 < std::numeric_limits<T>::epsilon ())) 300 retval = 0; 301 else 302 { 303 T n = std::floor (q); 304 305 // Prevent use of extra precision. 306 volatile T tmp = y * n; 307 308 retval = x - tmp; 309 } 310 } 311 312 if (x != y && y != 0) 313 retval = copysign (retval, y); 314 315 return retval; 316 } 317 318 template <typename T> 319 T rem(T x,T y)320 rem (T x, T y) 321 { 322 T retval; 323 324 if (y == 0) 325 retval = numeric_limits<T>::NaN (); 326 else 327 { 328 T q = x / y; 329 330 if (x_nint (y) != y 331 && (std::abs ((q - x_nint (q)) / x_nint (q)) 332 < std::numeric_limits<T>::epsilon ())) 333 retval = 0; 334 else 335 { 336 T n = std::trunc (q); 337 338 // Prevent use of extra precision. 339 volatile T tmp = y * n; 340 341 retval = x - tmp; 342 } 343 } 344 345 if (x != y && y != 0) 346 retval = copysign (retval, x); 347 348 return retval; 349 } 350 351 // Generic min, max definitions 352 template <typename T> 353 T min(T x,T y)354 min (T x, T y) 355 { 356 return x <= y ? x : y; 357 } 358 359 template <typename T> 360 T max(T x,T y)361 max (T x, T y) 362 { 363 return x >= y ? x : y; 364 } 365 366 // This form is favorable. GCC will translate (x <= y ? x : y) without a 367 // jump, hence the only conditional jump involved will be the first 368 // (isnan), infrequent and hence friendly to branch prediction. 369 370 inline double min(double x,double y)371 min (double x, double y) 372 { 373 return isnan (y) ? x : (x <= y ? x : y); 374 } 375 376 inline double max(double x,double y)377 max (double x, double y) 378 { 379 return isnan (y) ? x : (x >= y ? x : y); 380 } 381 382 inline float min(float x,float y)383 min (float x, float y) 384 { 385 return isnan (y) ? x : (x <= y ? x : y); 386 } 387 388 inline float max(float x,float y)389 max (float x, float y) 390 { 391 return isnan (y) ? x : (x >= y ? x : y); 392 } 393 394 inline std::complex<double> min(const std::complex<double> & x,const std::complex<double> & y)395 min (const std::complex<double>& x, const std::complex<double>& y) 396 { 397 return abs (x) <= abs (y) ? x : (isnan (x) ? x : y); 398 } 399 400 inline std::complex<float> min(const std::complex<float> & x,const std::complex<float> & y)401 min (const std::complex<float>& x, const std::complex<float>& y) 402 { 403 return abs (x) <= abs (y) ? x : (isnan (x) ? x : y); 404 } 405 406 inline std::complex<double> max(const std::complex<double> & x,const std::complex<double> & y)407 max (const std::complex<double>& x, const std::complex<double>& y) 408 { 409 return abs (x) >= abs (y) ? x : (isnan (x) ? x : y); 410 } 411 412 inline std::complex<float> max(const std::complex<float> & x,const std::complex<float> & y)413 max (const std::complex<float>& x, const std::complex<float>& y) 414 { 415 return abs (x) >= abs (y) ? x : (isnan (x) ? x : y); 416 } 417 418 template <typename T> 419 inline octave_int<T> min(const octave_int<T> & x,const octave_int<T> & y)420 min (const octave_int<T>& x, const octave_int<T>& y) 421 { 422 return xmin (x, y); 423 } 424 425 template <typename T> 426 inline octave_int<T> max(const octave_int<T> & x,const octave_int<T> & y)427 max (const octave_int<T>& x, const octave_int<T>& y) 428 { 429 return xmax (x, y); 430 } 431 432 // These map reals to Complex. 433 434 extern OCTAVE_API Complex rc_acos (double); 435 extern OCTAVE_API FloatComplex rc_acos (float); 436 437 extern OCTAVE_API Complex rc_acosh (double); 438 extern OCTAVE_API FloatComplex rc_acosh (float); 439 440 extern OCTAVE_API Complex rc_asin (double); 441 extern OCTAVE_API FloatComplex rc_asin (float); 442 443 extern OCTAVE_API Complex rc_atanh (double); 444 extern OCTAVE_API FloatComplex rc_atanh (float); 445 446 extern OCTAVE_API Complex rc_log (double); 447 extern OCTAVE_API FloatComplex rc_log (float); 448 449 extern OCTAVE_API Complex rc_log2 (double); 450 extern OCTAVE_API FloatComplex rc_log2 (float); 451 452 extern OCTAVE_API Complex rc_log10 (double); 453 extern OCTAVE_API FloatComplex rc_log10 (float); 454 455 extern OCTAVE_API Complex rc_sqrt (double); 456 extern OCTAVE_API FloatComplex rc_sqrt (float); 457 } 458 } 459 460 #endif 461