//------------------------------------------------------------------------------ // GB_math.h: definitions for complex types, and mathematical operators //------------------------------------------------------------------------------ // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2023, All Rights Reserved. // SPDX-License-Identifier: Apache-2.0 //------------------------------------------------------------------------------ #include "GB_casting.h" #ifndef GB_MATH_H #define GB_MATH_H //------------------------------------------------------------------------------ // integer division //------------------------------------------------------------------------------ // The GJ_idiv* definitions are used in JIT kernels only. // Integer division is done carefully so that GraphBLAS does not terminate the // user's application on divide-by-zero. To compute x/0: if x is zero, the // result is zero (like NaN). if x is negative, the result is the negative // integer with biggest magnitude (like -infinity). if x is positive, the // result is the biggest positive integer (like +infinity). inline int8_t GB_idiv_int8 (int8_t x, int8_t y) { // returns x/y when x and y are int8_t if (y == -1) { // INT32_MIN/(-1) causes floating point exception; avoid it return (-x) ; } else if (y == 0) { // zero divided by zero gives 'integer Nan' // x/0 where x is nonzero: result is integer -Inf or +Inf return ((x == 0) ? 0 : ((x < 0) ? INT8_MIN : INT8_MAX)) ; } else { // normal case for signed integer division return (x / y) ; } } #define GJ_idiv_int8_DEFN \ "int8_t GJ_idiv_int8 (int8_t x, int8_t y) \n" \ "{ \n" \ " if (y == -1) \n" \ " { \n" \ " return (-x) ; \n" \ " } \n" \ " else if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : ((x < 0) ? INT8_MIN : INT8_MAX)) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline int16_t GB_idiv_int16 (int16_t x, int16_t y) { // returns x/y when x and y are int16_t if (y == -1) { // INT32_MIN/(-1) causes floating point exception; avoid it return (-x) ; } else if (y == 0) { // zero divided by zero gives 'integer Nan' // x/0 where x is nonzero: result is integer -Inf or +Inf return ((x == 0) ? 0 : ((x < 0) ? INT16_MIN : INT16_MAX)) ; } else { // normal case for signed integer division return (x / y) ; } } #define GJ_idiv_int16_DEFN \ "int16_t GJ_idiv_int16 (int16_t x, int16_t y) \n" \ "{ \n" \ " if (y == -1) \n" \ " { \n" \ " return (-x) ; \n" \ " } \n" \ " else if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : ((x < 0) ? INT16_MIN : INT16_MAX)) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline int32_t GB_idiv_int32 (int32_t x, int32_t y) { // returns x/y when x and y are int32_t if (y == -1) { // INT32_MIN/(-1) causes floating point exception; avoid it return (-x) ; } else if (y == 0) { // zero divided by zero gives 'integer Nan' // x/0 where x is nonzero: result is integer -Inf or +Inf return ((x == 0) ? 0 : ((x < 0) ? INT32_MIN : INT32_MAX)) ; } else { // normal case for signed integer division return (x / y) ; } } #define GJ_idiv_int32_DEFN \ "int32_t GJ_idiv_int32 (int32_t x, int32_t y) \n" \ "{ \n" \ " if (y == -1) \n" \ " { \n" \ " return (-x) ; \n" \ " } \n" \ " else if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : ((x < 0) ? INT32_MIN : INT32_MAX)) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline int64_t GB_idiv_int64 (int64_t x, int64_t y) { // returns x/y when x and y are int64_t if (y == -1) { // INT32_MIN/(-1) causes floating point exception; avoid it return (-x) ; } else if (y == 0) { // zero divided by zero gives 'integer Nan' // x/0 where x is nonzero: result is integer -Inf or +Inf return ((x == 0) ? 0 : ((x < 0) ? INT64_MIN : INT64_MAX)) ; } else { // normal case for signed integer division return (x / y) ; } } #define GJ_idiv_int64_DEFN \ "int64_t GJ_idiv_int64 (int64_t x, int64_t y) \n" \ "{ \n" \ " if (y == -1) \n" \ " { \n" \ " return (-x) ; \n" \ " } \n" \ " else if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : ((x < 0) ? INT64_MIN : INT64_MAX)) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline uint8_t GB_idiv_uint8 (uint8_t x, uint8_t y) { if (y == 0) { // x/0: 0/0 is integer Nan, otherwise result is +Inf return ((x == 0) ? 0 : UINT8_MAX) ; } else { // normal case for unsigned integer division return (x / y) ; } } #define GJ_idiv_uint8_DEFN \ "uint8_t GJ_idiv_uint8 (uint8_t x, uint8_t y) \n" \ "{ \n" \ " if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : UINT8_MAX) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline uint16_t GB_idiv_uint16 (uint16_t x, uint16_t y) { if (y == 0) { // x/0: 0/0 is integer Nan, otherwise result is +Inf return ((x == 0) ? 0 : UINT16_MAX) ; } else { // normal case for unsigned integer division return (x / y) ; } } #define GJ_idiv_uint16_DEFN \ "uint16_t GJ_idiv_uint16 (uint16_t x, uint16_t y) \n" \ "{ \n" \ " if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : UINT16_MAX) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline uint32_t GB_idiv_uint32 (uint32_t x, uint32_t y) { if (y == 0) { // x/0: 0/0 is integer Nan, otherwise result is +Inf return ((x == 0) ? 0 : UINT32_MAX) ; } else { // normal case for unsigned integer division return (x / y) ; } } #define GJ_idiv_uint32_DEFN \ "uint32_t GJ_idiv_uint32 (uint32_t x, uint32_t y) \n" \ "{ \n" \ " if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : UINT32_MAX) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" inline uint64_t GB_idiv_uint64 (uint64_t x, uint64_t y) { if (y == 0) { // x/0: 0/0 is integer Nan, otherwise result is +Inf return ((x == 0) ? 0 : UINT64_MAX) ; } else { // normal case for unsigned integer division return (x / y) ; } } #define GJ_idiv_uint64_DEFN \ "uint64_t GJ_idiv_uint64 (uint64_t x, uint64_t y) \n" \ "{ \n" \ " if (y == 0) \n" \ " { \n" \ " return ((x == 0) ? 0 : UINT64_MAX) ; \n" \ " } \n" \ " else \n" \ " { \n" \ " return (x / y) ; \n" \ " } \n" \ "}" //------------------------------------------------------------------------------ // complex division //------------------------------------------------------------------------------ // The GJ_FC*_div definitions are used in JIT kernels only. // complex division is problematic. It is not supported at all on MS Visual // Studio. With other compilers, complex division exists but it has different // NaN and Inf behavior as compared with MATLAB, which causes the tests to // fail. As a result, the built-in complex division is not used, even if the // compiler supports it. // Three cases below are from ACM Algo 116, R. L. Smith, 1962. inline GxB_FC64_t GB_FC64_div (GxB_FC64_t x, GxB_FC64_t y) { double xr = GB_creal (x) ; double xi = GB_cimag (x) ; double yr = GB_creal (y) ; double yi = GB_cimag (y) ; int yr_class = fpclassify (yr) ; int yi_class = fpclassify (yi) ; if (yi_class == FP_ZERO) { // (zr,zi) = (xr,xi) / (yr,0) return (GB_CMPLX64 (xr / yr, xi / yr)) ; } else if (yr_class == FP_ZERO) { // (zr,zi) = (xr,xi) / (0,yi) = (xi,-xr) / (yi,0) return (GB_CMPLX64 (xi / yi, -xr / yi)) ; } else if (yi_class == FP_INFINITE && yr_class == FP_INFINITE) { // Using Smith's method for a very special case double r = (signbit (yr) == signbit (yi)) ? (1) : (-1) ; double d = yr + r * yi ; return (GB_CMPLX64 ((xr + xi * r) / d, (xi - xr * r) / d)) ; } else if (fabs (yr) >= fabs (yi)) { // Smith's method (1st case) double r = yi / yr ; double d = yr + r * yi ; return (GB_CMPLX64 ((xr + xi * r) / d, (xi - xr * r) / d)) ; } else { // Smith's method (2nd case) double r = yr / yi ; double d = r * yr + yi ; return (GB_CMPLX64 ((xr * r + xi) / d, (xi * r - xr) / d)) ; } } #define GJ_FC64_div_DEFN \ "GxB_FC64_t GJ_FC64_div (GxB_FC64_t x, GxB_FC64_t y) \n" \ "{ \n" \ " double xr = GB_creal (x) ; \n" \ " double xi = GB_cimag (x) ; \n" \ " double yr = GB_creal (y) ; \n" \ " double yi = GB_cimag (y) ; \n" \ " int yr_class = fpclassify (yr) ; \n" \ " int yi_class = fpclassify (yi) ; \n" \ " if (yi_class == FP_ZERO) \n" \ " { \n" \ " return (GJ_CMPLX64 (xr / yr, xi / yr)) ; \n" \ " } \n" \ " else if (yr_class == FP_ZERO) \n" \ " { \n" \ " return (GJ_CMPLX64 (xi / yi, -xr / yi)) ; \n" \ " } \n" \ " else if (yi_class == FP_INFINITE && yr_class == FP_INFINITE) \n" \ " { \n" \ " double r = (signbit (yr) == signbit (yi)) ? (1) : (-1) ; \n" \ " double d = yr + r * yi ; \n" \ " return (GJ_CMPLX64 ((xr + xi * r) / d, (xi - xr * r) / d)) ;\n" \ " } \n" \ " else if (fabs (yr) >= fabs (yi)) \n" \ " { \n" \ " double r = yi / yr ; \n" \ " double d = yr + r * yi ; \n" \ " return (GJ_CMPLX64 ((xr + xi * r) / d, (xi - xr * r) / d)) ;\n" \ " } \n" \ " else \n" \ " { \n" \ " double r = yr / yi ; \n" \ " double d = r * yr + yi ; \n" \ " return (GJ_CMPLX64 ((xr * r + xi) / d, (xi * r - xr) / d)) ;\n" \ " } \n" \ "}" inline GxB_FC32_t GB_FC32_div (GxB_FC32_t x, GxB_FC32_t y) { // single complex division: cast double complex, do the division, // and then cast back to single complex. double xr = (double) GB_crealf (x) ; double xi = (double) GB_cimagf (x) ; double yr = (double) GB_crealf (y) ; double yi = (double) GB_cimagf (y) ; GxB_FC64_t zz = GB_FC64_div (GB_CMPLX64 (xr, xi), GB_CMPLX64 (yr, yi)) ; return (GB_CMPLX32 ((float) GB_creal (zz), (float) GB_cimag (zz))) ; } #define GJ_FC32_div_DEFN \ "GxB_FC32_t GJ_FC32_div (GxB_FC32_t x, GxB_FC32_t y) \n" \ "{ \n" \ " double xr = (double) GB_crealf (x) ; \n" \ " double xi = (double) GB_cimagf (x) ; \n" \ " double yr = (double) GB_crealf (y) ; \n" \ " double yi = (double) GB_cimagf (y) ; \n" \ " GxB_FC64_t zz ; \n" \ " zz = GJ_FC64_div (GJ_CMPLX64 (xr, xi), GJ_CMPLX64 (yr, yi)) ; \n" \ " return (GJ_CMPLX32 ((float) GB_creal(zz), (float) GB_cimag(zz))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = x^y: wrappers for pow, powf, cpow, and cpowf //------------------------------------------------------------------------------ // if x or y are NaN, then z is NaN // if y is zero, then z is 1 // if (x and y are complex but with zero imaginary parts, and // (x >= 0 or if y is an integer, NaN, or Inf)), then z is real // else use the built-in C library function, z = pow (x,y) inline float GB_powf (float x, float y) { int xr_class = fpclassify (x) ; int yr_class = fpclassify (y) ; if (xr_class == FP_NAN || yr_class == FP_NAN) { // z is nan if either x or y are nan return (NAN) ; } if (yr_class == FP_ZERO) { // z is 1 if y is zero return (1) ; } // otherwise, z = powf (x,y) return (powf (x, y)) ; } #define GJ_powf_DEFN \ "float GJ_powf (float x, float y) \n" \ "{ \n" \ " int xr_class = fpclassify (x) ; \n" \ " int yr_class = fpclassify (y) ; \n" \ " if (xr_class == FP_NAN || yr_class == FP_NAN) \n" \ " { \n" \ " return (NAN) ; \n" \ " } \n" \ " if (yr_class == FP_ZERO) \n" \ " { \n" \ " return (1) ; \n" \ " } \n" \ " return (powf (x, y)) ; \n" \ "}" inline double GB_pow (double x, double y) { int xr_class = fpclassify (x) ; int yr_class = fpclassify (y) ; if (xr_class == FP_NAN || yr_class == FP_NAN) { // z is nan if either x or y are nan return (NAN) ; } if (yr_class == FP_ZERO) { // z is 1 if y is zero return (1) ; } // otherwise, z = pow (x,y) return (pow (x, y)) ; } #define GJ_pow_DEFN \ "double GJ_pow (double x, double y) \n" \ "{ \n" \ " int xr_class = fpclassify (x) ; \n" \ " int yr_class = fpclassify (y) ; \n" \ " if (xr_class == FP_NAN || yr_class == FP_NAN) \n" \ " { \n" \ " // z is nan if either x or y are nan \n" \ " return (NAN) ; \n" \ " } \n" \ " if (yr_class == FP_ZERO) \n" \ " { \n" \ " // z is 1 if y is zero \n" \ " return (1) ; \n" \ " } \n" \ " // otherwise, z = pow (x,y) \n" \ " return (pow (x, y)) ; \n" \ "}" inline GxB_FC32_t GB_FC32_pow (GxB_FC32_t x, GxB_FC32_t y) { float xr = GB_crealf (x) ; float yr = GB_crealf (y) ; int xr_class = fpclassify (xr) ; int yr_class = fpclassify (yr) ; int xi_class = fpclassify (GB_cimagf (x)) ; int yi_class = fpclassify (GB_cimagf (y)) ; if (xi_class == FP_ZERO && yi_class == FP_ZERO) { // both x and y are real; see if z should be real if (xr >= 0 || yr_class == FP_NAN || yr_class == FP_INFINITE || yr == truncf (yr)) { // z is real if x >= 0, or if y is an integer, NaN, or Inf return (GB_CMPLX32 (GB_powf (xr, yr), 0)) ; } } if (xr_class == FP_NAN || xi_class == FP_NAN || yr_class == FP_NAN || yi_class == FP_NAN) { // z is (nan,nan) if any part of x or y are nan return (GB_CMPLX32 (NAN, NAN)) ; } if (yr_class == FP_ZERO && yi_class == FP_ZERO) { // z is (1,0) if y is (0,0) return (GxB_CMPLXF (1, 0)) ; } return (GB_cpowf (x, y)) ; } #define GJ_FC32_pow_DEFN \ "GxB_FC32_t GJ_FC32_pow (GxB_FC32_t x, GxB_FC32_t y) \n" \ "{ \n" \ " float xr = GB_crealf (x) ; \n" \ " float yr = GB_crealf (y) ; \n" \ " int xr_class = fpclassify (xr) ; \n" \ " int yr_class = fpclassify (yr) ; \n" \ " int xi_class = fpclassify (GB_cimagf (x)) ; \n" \ " int yi_class = fpclassify (GB_cimagf (y)) ; \n" \ " if (xi_class == FP_ZERO && yi_class == FP_ZERO) \n" \ " { \n" \ " if (xr >= 0 || yr_class == FP_NAN || \n" \ " yr_class == FP_INFINITE || yr == truncf (yr)) \n" \ " { \n" \ " return (GJ_CMPLX32 (GJ_powf (xr, yr), 0)) ; \n" \ " } \n" \ " } \n" \ " if (xr_class == FP_NAN || xi_class == FP_NAN || \n" \ " yr_class == FP_NAN || yi_class == FP_NAN) \n" \ " { \n" \ " return (GJ_CMPLX32 (NAN, NAN)) ; \n" \ " } \n" \ " if (yr_class == FP_ZERO && yi_class == FP_ZERO) \n" \ " { \n" \ " return (GxB_CMPLXF (1, 0)) ; \n" \ " } \n" \ " return (GB_cpowf (x, y)) ; \n" \ "}" inline GxB_FC64_t GB_FC64_pow (GxB_FC64_t x, GxB_FC64_t y) { double xr = GB_creal (x) ; double yr = GB_creal (y) ; int xr_class = fpclassify (xr) ; int yr_class = fpclassify (yr) ; int xi_class = fpclassify (GB_cimag (x)) ; int yi_class = fpclassify (GB_cimag (y)) ; if (xi_class == FP_ZERO && yi_class == FP_ZERO) { // both x and y are real; see if z should be real if (xr >= 0 || yr_class == FP_NAN || yr_class == FP_INFINITE || yr == trunc (yr)) { // z is real if x >= 0, or if y is an integer, NaN, or Inf return (GB_CMPLX64 (GB_pow (xr, yr), 0)) ; } } if (xr_class == FP_NAN || xi_class == FP_NAN || yr_class == FP_NAN || yi_class == FP_NAN) { // z is (nan,nan) if any part of x or y are nan return (GB_CMPLX64 (NAN, NAN)) ; } if (yr_class == FP_ZERO && yi_class == FP_ZERO) { // z is (1,0) if y is (0,0) return (GxB_CMPLX (1, 0)) ; } return (GB_cpow (x, y)) ; } #define GJ_FC64_pow_DEFN \ "GxB_FC64_t GJ_FC64_pow (GxB_FC64_t x, GxB_FC64_t y) \n" \ "{ \n" \ " double xr = GB_creal (x) ; \n" \ " double yr = GB_creal (y) ; \n" \ " int xr_class = fpclassify (xr) ; \n" \ " int yr_class = fpclassify (yr) ; \n" \ " int xi_class = fpclassify (GB_cimag (x)) ; \n" \ " int yi_class = fpclassify (GB_cimag (y)) ; \n" \ " if (xi_class == FP_ZERO && yi_class == FP_ZERO) \n" \ " { \n" \ " if (xr >= 0 || yr_class == FP_NAN || \n" \ " yr_class == FP_INFINITE || yr == trunc (yr)) \n" \ " { \n" \ " return (GJ_CMPLX64 (GJ_pow (xr, yr), 0)) ; \n" \ " } \n" \ " } \n" \ " if (xr_class == FP_NAN || xi_class == FP_NAN || \n" \ " yr_class == FP_NAN || yi_class == FP_NAN) \n" \ " { \n" \ " return (GJ_CMPLX64 (NAN, NAN)) ; \n" \ " } \n" \ " if (yr_class == FP_ZERO && yi_class == FP_ZERO) \n" \ " { \n" \ " return (GxB_CMPLX (1, 0)) ; \n" \ " } \n" \ " return (GB_cpow (x, y)) ; \n" \ "}" inline int8_t GB_pow_int8 (int8_t x, int8_t y) { return (GB_cast_to_int8_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_int8_DEFN \ "int8_t GJ_pow_int8 (int8_t x, int8_t y) \n" \ "{ \n" \ " return (GJ_cast_to_int8 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline int16_t GB_pow_int16 (int16_t x, int16_t y) { return (GB_cast_to_int16_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_int16_DEFN \ "int16_t GJ_pow_int16 (int16_t x, int16_t y) \n" \ "{ \n" \ " return (GJ_cast_to_int16 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline int32_t GB_pow_int32 (int32_t x, int32_t y) { return (GB_cast_to_int32_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_int32_DEFN \ "int32_t GJ_pow_int32 (int32_t x, int32_t y) \n" \ "{ \n" \ " return (GJ_cast_to_int32 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline int64_t GB_pow_int64 (int64_t x, int64_t y) { return (GB_cast_to_int64_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_int64_DEFN \ "int64_t GJ_pow_int64 (int64_t x, int64_t y) \n" \ "{ \n" \ " return (GJ_cast_to_int64 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline uint8_t GB_pow_uint8 (uint8_t x, uint8_t y) { return (GB_cast_to_uint8_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_uint8_DEFN \ "int8_t GJ_pow_uint8 (int8_t x, int8_t y) \n" \ "{ \n" \ " return (GJ_cast_to_uint8 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline uint16_t GB_pow_uint16 (uint16_t x, uint16_t y) { return (GB_cast_to_uint16_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_uint16_DEFN \ "int16_t GJ_pow_uint16 (int16_t x, int16_t y) \n" \ "{ \n" \ " return (GJ_cast_to_uint16 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline uint32_t GB_pow_uint32 (uint32_t x, uint32_t y) { return (GB_cast_to_uint32_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_uint32_DEFN \ "int32_t GJ_pow_uint32 (int32_t x, int32_t y) \n" \ "{ \n" \ " return (GJ_cast_to_uint32 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" inline uint64_t GB_pow_uint64 (uint64_t x, uint64_t y) { return (GB_cast_to_uint64_t (GB_pow ((double) x, (double) y))) ; } #define GJ_pow_uint64_DEFN \ "int64_t GJ_pow_uint64 (int64_t x, int64_t y) \n" \ "{ \n" \ " return (GJ_cast_to_uint64 (GJ_pow ((double) x, (double) y))) ; \n" \ "}" //------------------------------------------------------------------------------ // frexp for float and double //------------------------------------------------------------------------------ inline float GB_frexpxf (float x) { // ignore the exponent and just return the mantissa int exp_ignored ; return (frexpf (x, &exp_ignored)) ; } #define GJ_frexpxf_DEFN \ "float GJ_frexpxf (float x) \n" \ "{ \n" \ " int exp_ignored ; \n" \ " return (frexpf (x, &exp_ignored)) ; \n" \ "}" inline float GB_frexpef (float x) { // ignore the mantissa and just return the exponent int exp ; (void) frexpf (x, &exp) ; return ((float) exp) ; } #define GJ_frexpef_DEFN \ "float GJ_frexpef (float x) \n" \ "{ \n" \ " int exp ; \n" \ " (void) frexpf (x, &exp) ; \n" \ " return ((float) exp) ; \n" \ "}" inline double GB_frexpx (double x) { // ignore the exponent and just return the mantissa int exp_ignored ; return (frexp (x, &exp_ignored)) ; } #define GJ_frexpx_DEFN \ "double GJ_frexpx (double x) \n" \ "{ \n" \ " int exp_ignored ; \n" \ " return (frexp (x, &exp_ignored)) ; \n" \ "}" inline double GB_frexpe (double x) { // ignore the mantissa and just return the exponent int exp ; (void) frexp (x, &exp) ; return ((double) exp) ; } #define GJ_frexpe_DEFN \ "double GJ_frexpe (double x) \n" \ "{ \n" \ " int exp ; \n" \ " (void) frexp (x, &exp) ; \n" \ " return ((double) exp) ; \n" \ "}" //------------------------------------------------------------------------------ // signum functions //------------------------------------------------------------------------------ inline float GB_signumf (float x) { if (isnan (x)) return (x) ; return ((float) ((x < 0) ? (-1) : ((x > 0) ? 1 : 0))) ; } #define GJ_signumf_DEFN \ "float GJ_signumf (float x) \n" \ "{ \n" \ " if (isnan (x)) return (x) ; \n" \ " return ((float) ((x < 0) ? (-1) : ((x > 0) ? 1 : 0))) ; \n" \ "}" inline double GB_signum (double x) { if (isnan (x)) return (x) ; return ((double) ((x < 0) ? (-1) : ((x > 0) ? 1 : 0))) ; } #define GJ_signum_DEFN \ "double GJ_signum (double x) \n" \ "{ \n" \ " if (isnan (x)) return (x) ; \n" \ " return ((double) ((x < 0) ? (-1) : ((x > 0) ? 1 : 0))) ; \n" \ "}" inline GxB_FC32_t GB_csignumf (GxB_FC32_t x) { if (GB_crealf (x) == 0 && GB_cimagf (x) == 0) { return (GxB_CMPLXF (0,0)) ; } float y = GB_cabsf (x) ; return (GB_CMPLX32 (GB_crealf (x) / y, GB_cimagf (x) / y)) ; } #define GJ_csignumf_DEFN \ "GxB_FC32_t GJ_csignumf (GxB_FC32_t x) \n" \ "{ \n" \ " if (GB_crealf (x) == 0 && GB_cimagf (x) == 0) \n" \ " { \n" \ " return (GxB_CMPLXF (0,0)) ; \n" \ " } \n" \ " float y = GB_cabsf (x) ; \n" \ " return (GJ_CMPLX32 (GB_crealf (x) / y, GB_cimagf (x) / y)) ; \n" \ "}" inline GxB_FC64_t GB_csignum (GxB_FC64_t x) { if (GB_creal (x) == 0 && GB_cimag (x) == 0) { return (GxB_CMPLX (0,0)) ; } double y = GB_cabs (x) ; return (GB_CMPLX64 (GB_creal (x) / y, GB_cimag (x) / y)) ; } #define GJ_csignum_DEFN \ "GxB_FC64_t GJ_csignum (GxB_FC64_t x) \n" \ "{ \n" \ " if (GB_creal (x) == 0 && GB_cimag (x) == 0) \n" \ " { \n" \ " return (GxB_CMPLX (0,0)) ; \n" \ " } \n" \ " double y = GB_cabs (x) ; \n" \ " return (GJ_CMPLX64 (GB_creal (x) / y, GB_cimag (x) / y)) ; \n" \ "}" //------------------------------------------------------------------------------ // complex functions //------------------------------------------------------------------------------ // The C11 math.h header defines the ceil, floor, round, trunc, // exp2, expm1, log10, log1pm, or log2 functions for float and double, // but the corresponding functions do not appear in the C11 complex.h. // These functions are used instead, for float complex and double complex. //------------------------------------------------------------------------------ // z = ceil (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_cceilf (GxB_FC32_t x) { return (GB_CMPLX32 (ceilf (GB_crealf (x)), ceilf (GB_cimagf (x)))) ; } #define GJ_cceilf_DEFN \ "GxB_FC32_t GJ_cceilf (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_CMPLX32 (ceilf (GB_crealf (x)), ceilf (GB_cimagf (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = ceil (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_cceil (GxB_FC64_t x) { return (GB_CMPLX64 (ceil (GB_creal (x)), ceil (GB_cimag (x)))) ; } #define GJ_cceil_DEFN \ "GxB_FC64_t GJ_cceil (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_CMPLX64 (ceil (GB_creal (x)), ceil (GB_cimag (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = floor (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_cfloorf (GxB_FC32_t x) { return (GB_CMPLX32 (floorf (GB_crealf (x)), floorf (GB_cimagf (x)))) ; } #define GJ_cfloorf_DEFN \ "GxB_FC32_t GJ_cfloorf (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_CMPLX32 (floorf (GB_crealf (x)), floorf (GB_cimagf (x)))) ;\n" \ "}" //------------------------------------------------------------------------------ // z = floor (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_cfloor (GxB_FC64_t x) { return (GB_CMPLX64 (floor (GB_creal (x)), floor (GB_cimag (x)))) ; } #define GJ_cfloor_DEFN \ "GxB_FC64_t GJ_cfloor (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_CMPLX64 (floor (GB_creal (x)), floor (GB_cimag (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = round (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_croundf (GxB_FC32_t x) { return (GB_CMPLX32 (roundf (GB_crealf (x)), roundf (GB_cimagf (x)))) ; } #define GJ_croundf_DEFN \ "GxB_FC32_t GJ_croundf (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_CMPLX32 (roundf (GB_crealf (x)), roundf (GB_cimagf (x)))) ;\n" \ "}" //------------------------------------------------------------------------------ // z = round (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_cround (GxB_FC64_t x) { return (GB_CMPLX64 (round (GB_creal (x)), round (GB_cimag (x)))) ; } #define GJ_cround_DEFN \ "GxB_FC64_t GJ_cround (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_CMPLX64 (round (GB_creal (x)), round (GB_cimag (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = trunc (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_ctruncf (GxB_FC32_t x) { return (GB_CMPLX32 (truncf (GB_crealf (x)), truncf (GB_cimagf (x)))) ; } #define GJ_ctruncf_DEFN \ "GxB_FC32_t GJ_ctruncf (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_CMPLX32 (truncf (GB_crealf (x)), truncf (GB_cimagf (x)))) ;\n" \ "}" //------------------------------------------------------------------------------ // z = trunc (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_ctrunc (GxB_FC64_t x) { return (GB_CMPLX64 (trunc (GB_creal (x)), trunc (GB_cimag (x)))) ; } #define GJ_ctrunc_DEFN \ "GxB_FC64_t GJ_ctrunc (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_CMPLX64 (trunc (GB_creal (x)), trunc (GB_cimag (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = exp2 (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_cexp2f (GxB_FC32_t x) { if (fpclassify (GB_cimagf (x)) == FP_ZERO) { // x is real, use exp2f return (GB_CMPLX32 (exp2f (GB_crealf (x)), 0)) ; } return (GB_FC32_pow (GxB_CMPLXF (2,0), x)) ; // z = 2^x } #define GJ_cexp2f_DEFN \ "GxB_FC32_t GJ_cexp2f (GxB_FC32_t x) \n" \ "{ \n" \ " if (fpclassify (GB_cimagf (x)) == FP_ZERO) \n" \ " { \n" \ " return (GJ_CMPLX32 (exp2f (GB_crealf (x)), 0)) ; \n" \ " } \n" \ " return (GJ_FC32_pow (GxB_CMPLXF (2,0), x)) ; \n" \ "}" //------------------------------------------------------------------------------ // z = exp2 (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_cexp2 (GxB_FC64_t x) { if (fpclassify (GB_cimag (x)) == FP_ZERO) { // x is real, use exp2 return (GB_CMPLX64 (exp2 (GB_creal (x)), 0)) ; } return (GB_FC64_pow (GxB_CMPLX (2,0), x)) ; // z = 2^x } #define GJ_cexp2_DEFN \ "GxB_FC64_t GJ_cexp2 (GxB_FC64_t x) \n" \ "{ \n" \ " if (fpclassify (GB_cimag (x)) == FP_ZERO) \n" \ " { \n" \ " return (GJ_CMPLX64 (exp2 (GB_creal (x)), 0)) ; \n" \ " } \n" \ " return (GJ_FC64_pow (GxB_CMPLX (2,0), x)) ; \n" \ "}" //------------------------------------------------------------------------------ // z = expm1 (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_cexpm1 (GxB_FC64_t x) { // FUTURE: GB_cexpm1 is not accurate // z = cexp (x) - 1 GxB_FC64_t z = GB_cexp (x) ; return (GB_CMPLX64 (GB_creal (z) - 1, GB_cimag (z))) ; } #define GJ_cexpm1_DEFN \ "GxB_FC64_t GJ_cexpm1 (GxB_FC64_t x) \n" \ "{ \n" \ " GxB_FC64_t z = GB_cexp (x) ; \n" \ " return (GJ_CMPLX64 (GB_creal (z) - 1, GB_cimag (z))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = expm1 (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_cexpm1f (GxB_FC32_t x) { // typecast to double and use GB_cexpm1 GxB_FC64_t z = GB_CMPLX64 ((double) GB_crealf (x), (double) GB_cimagf (x)) ; z = GB_cexpm1 (z) ; return (GB_CMPLX32 ((float) GB_creal (z), (float) GB_cimag (z))) ; } #define GJ_cexpm1f_DEFN \ "GxB_FC32_t GJ_cexpm1f (GxB_FC32_t x) \n" \ "{ \n" \ " GxB_FC64_t z = GJ_CMPLX64 ((double) GB_crealf (x), \n" \ " (double) GB_cimagf (x)) ; \n" \ " z = GJ_cexpm1 (z) ; \n" \ " return (GJ_CMPLX32 ((float) GB_creal (z), \n" \ " (float) GB_cimag (z))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log1p (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_clog1p (GxB_FC64_t x) { // FUTURE: GB_clog1p is not accurate // z = clog (1+x) return (GB_clog (GB_CMPLX64 (GB_creal (x) + 1, GB_cimag (x)))) ; } #define GJ_clog1p_DEFN \ "GxB_FC64_t GJ_clog1p (GxB_FC64_t x) \n" \ "{ \n" \ " return (GB_clog (GJ_CMPLX64 (GB_creal (x) + 1, GB_cimag (x)))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log1p (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_clog1pf (GxB_FC32_t x) { // typecast to double and use GB_clog1p GxB_FC64_t z = GB_CMPLX64 ((double) GB_crealf (x), (double) GB_cimagf (x)) ; z = GB_clog1p (z) ; return (GB_CMPLX32 ((float) GB_creal (z), (float) GB_cimag (z))) ; } #define GJ_clog1pf_DEFN \ "GxB_FC32_t GJ_clog1pf (GxB_FC32_t x) \n" \ "{ \n" \ " GxB_FC64_t z = GJ_CMPLX64 ((double) GB_crealf (x), \n" \ " (double) GB_cimagf (x)) ; \n" \ " z = GJ_clog1p (z) ; \n" \ " return (GJ_CMPLX32 ((float) GB_creal (z), \n" \ " (float) GB_cimag (z))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log10 (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_clog10f (GxB_FC32_t x) { // z = log (x) / log (10) return (GB_FC32_div (GB_clogf (x), GxB_CMPLXF (2.3025851f, 0))) ; } #define GJ_clog10f_DEFN \ "GxB_FC32_t GJ_clog10f (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_FC32_div (GB_clogf (x), GxB_CMPLXF (2.3025851f, 0))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log10 (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_clog10 (GxB_FC64_t x) { // z = log (x) / log (10) return (GB_FC64_div (GB_clog (x), GxB_CMPLX (2.302585092994045901, 0))) ; } #define GJ_clog10_DEFN \ "GxB_FC64_t GJ_clog10 (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_FC64_div (GB_clog (x), \n" \ " GxB_CMPLX (2.302585092994045901, 0))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log2 (x) for float complex //------------------------------------------------------------------------------ inline GxB_FC32_t GB_clog2f (GxB_FC32_t x) { // z = log (x) / log (2) return (GB_FC32_div (GB_clogf (x), GxB_CMPLXF (0.69314718f, 0))) ; } #define GJ_clog2f_DEFN \ "GxB_FC32_t GJ_clog2f (GxB_FC32_t x) \n" \ "{ \n" \ " return (GJ_FC32_div (GB_clogf (x), GxB_CMPLXF (0.69314718f, 0))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = log2 (x) for double complex //------------------------------------------------------------------------------ inline GxB_FC64_t GB_clog2 (GxB_FC64_t x) { // z = log (x) / log (2) return (GB_FC64_div (GB_clog (x), GxB_CMPLX (0.693147180559945286, 0))) ; } #define GJ_clog2_DEFN \ "GxB_FC64_t GJ_clog2 (GxB_FC64_t x) \n" \ "{ \n" \ " return (GJ_FC64_div (GB_clog (x), \n" \ " GxB_CMPLX (0.693147180559945286, 0))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = isinf (x) for float complex //------------------------------------------------------------------------------ inline bool GB_cisinff (GxB_FC32_t x) { return (isinf (GB_crealf (x)) || isinf (GB_cimagf (x))) ; } #define GJ_cisinff_DEFN \ "bool GJ_cisinff (GxB_FC32_t x) \n" \ "{ \n" \ " return (isinf (GB_crealf (x)) || isinf (GB_cimagf (x))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = isinf (x) for double complex //------------------------------------------------------------------------------ inline bool GB_cisinf (GxB_FC64_t x) { return (isinf (GB_creal (x)) || isinf (GB_cimag (x))) ; } #define GJ_cisinf_DEFN \ "bool GJ_cisinf (GxB_FC64_t x) \n" \ "{ \n" \ " return (isinf (GB_creal (x)) || isinf (GB_cimag (x))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = isnan (x) for float complex //------------------------------------------------------------------------------ inline bool GB_cisnanf (GxB_FC32_t x) { return (isnan (GB_crealf (x)) || isnan (GB_cimagf (x))) ; } #define GJ_cisnanf_DEFN \ "bool GJ_cisnanf (GxB_FC32_t x) \n" \ "{ \n" \ " return (isnan (GB_crealf (x)) || isnan (GB_cimagf (x))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = isnan (x) for double complex //------------------------------------------------------------------------------ inline bool GB_cisnan (GxB_FC64_t x) { return (isnan (GB_creal (x)) || isnan (GB_cimag (x))) ; } #define GJ_cisnan_DEFN \ "bool GJ_cisnan (GxB_FC64_t x) \n" \ "{ \n" \ " return (isnan (GB_creal (x)) || isnan (GB_cimag (x))) ; \n" \ "}" //------------------------------------------------------------------------------ // z = isfinite (x) for float complex //------------------------------------------------------------------------------ inline bool GB_cisfinitef (GxB_FC32_t x) { return (isfinite (GB_crealf (x)) && isfinite (GB_cimagf (x))) ; } #define GJ_cisfinitef_DEFN \ "bool GJ_cisfinitef (GxB_FC32_t x) \n" \ "{ \n" \ " return (isfinite (GB_crealf (x)) && isfinite (GB_cimagf (x))) ;\n" \ "}" //------------------------------------------------------------------------------ // z = isfinite (x) for double complex //------------------------------------------------------------------------------ inline bool GB_cisfinite (GxB_FC64_t x) { return (isfinite (GB_creal (x)) && isfinite (GB_cimag (x))) ; } #define GJ_cisfinite_DEFN \ "bool GJ_cisfinite (GxB_FC64_t x) \n" \ "{ \n" \ " return (isfinite (GB_creal (x)) && isfinite (GB_cimag (x))) ; \n" \ "}" #endif