//---------------------------------------------------------------------------// // MIT License // // Copyright (c) 2017 StreamComputing // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this software and associated documentation files (the "Software"), to deal // in the Software without restriction, including without limitation the rights // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell // copies of the Software, and to permit persons to whom the Software is // furnished to do so, subject to the following conditions: // // The above copyright notice and this permission notice shall be included in all // copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE // SOFTWARE. //---------------------------------------------------------------------------// #ifndef OPENCL_COMPLEX_MATH #define OPENCL_COMPLEX_MATH #define CONCAT(x, y) x##y #define FNAME(name, sufix) c##name##sufix // float2 #define clrealf(complex) complex.x; #define climagf(complex) complex.y; // double2 #define clreal(complex) complex.x; #define climag(complex) complex.y; #define OPENCL_COMPLEX_MATH_FUNCS(complex_type, real_type, func_sufix, math_consts_sufix) \ complex_type CONCAT(complex, func_sufix)(real_type r, real_type i) { \ return (complex_type)(r, i); \ } \ \ complex_type FNAME(add, func_sufix)(complex_type x, complex_type y) { \ return x + y; \ } \ \ complex_type FNAME(sub, func_sufix)(complex_type x, complex_type y) { \ return x - y; \ } \ \ complex_type FNAME(add_real, func_sufix)(complex_type z, real_type r) { \ return (complex_type)(z.x + r, z.y); \ } \ \ complex_type FNAME(sub_real, func_sufix)(complex_type z, real_type r) { \ return (complex_type)(z.x - r, z.y); \ } \ \ real_type FNAME(abs, func_sufix)(complex_type z) { \ return length(z); \ } \ \ real_type FNAME(arg, func_sufix)(complex_type z) { \ return atan2(z.y, z.x); \ } \ \ complex_type FNAME(mul, func_sufix)(complex_type z1, complex_type z2) { \ real_type x1 = z1.x; \ real_type y1 = z1.y; \ real_type x2 = z2.x; \ real_type y2 = z2.y; \ return (complex_type)(x1 * x2 - y1 * y2, x1 * y2 + x2 * y1); \ } \ \ complex_type FNAME(div, func_sufix)(complex_type z1, complex_type z2) { \ real_type x1 = z1.x; \ real_type y1 = z1.y; \ real_type x2 = z2.x; \ real_type y2 = z2.y; \ real_type iabs_z2 = CONCAT(1.0, func_sufix) / FNAME(abs, func_sufix)(z2); \ return (complex_type)(((x1 * x2 * iabs_z2) + (y1 * y2 * iabs_z2)) * iabs_z2, \ ((y1 * x2 * iabs_z2) - (x1 * y2 * iabs_z2)) * iabs_z2); \ } \ \ complex_type FNAME(mul_real, func_sufix)(complex_type z, real_type r) { \ return z * r; \ } \ \ complex_type FNAME(div_real, func_sufix)(complex_type z, real_type r) { \ return z / r; \ } \ \ complex_type FNAME(conj, func_sufix)(complex_type z) { \ return (complex_type)(z.x, -z.y); \ } \ \ complex_type FNAME(proj, func_sufix)(complex_type z) { \ if (isinf(z.x) || isinf(z.y)) { \ return (complex_type)(INFINITY, (copysign(CONCAT(0.0, func_sufix), z.y))); \ } \ return z; \ } \ \ real_type FNAME(norm, func_sufix)(complex_type z) { \ /* Returns the squared magnitude of the complex number z. */ \ /* The norm calculated by this function is also known as */ \ /* field norm or absolute square. */ \ real_type x = z.x; \ real_type y = z.y; \ return x * x + y * y; \ } \ \ complex_type FNAME(polar, func_sufix)(real_type r, real_type theta) { \ /* Returns a complex number with magnitude r and phase angle theta. */ \ return (complex_type)(r * cos(theta), r * sin(theta)); \ } \ \ complex_type FNAME(exp, func_sufix)(complex_type z) { \ /* The complex exponential function e^z for z = x+i*y */ \ /* equals to e^x * cis(y), */ \ /* or, e^x * (cos(y) + i*sin(y)) */ \ real_type expx = exp(z.x); \ return (complex_type)(expx * cos(z.y), expx * sin(z.y)); \ } \ \ complex_type FNAME(log, func_sufix)(complex_type z) { \ /* log(z) = log(abs(z)) + i * arg(z) */ \ return (complex_type)(log(FNAME(abs, func_sufix)(z)), FNAME(arg, func_sufix)(z)); \ } \ \ complex_type FNAME(log10, func_sufix)(complex_type z) { \ return FNAME(log, func_sufix)(z) / log(CONCAT(10.0, func_sufix)); \ } \ \ complex_type FNAME(pow, func_sufix)(complex_type z1, complex_type z2) { \ /* (z1)^(z2) = exp(z2 * log(z1)) = cexp(mul(z2, clog(z1))) */ \ return FNAME(exp, func_sufix)(FNAME(mul, func_sufix)(z2, FNAME(log, func_sufix)(z1))); \ } \ \ complex_type FNAME(sqrt, func_sufix)(complex_type z) { \ /* */ \ real_type x = z.x; \ real_type y = z.y; \ if (x == CONCAT(0.0, func_sufix)) { \ real_type t = sqrt(fabs(y) / 2); \ return (complex_type)(t, y < CONCAT(0.0, func_sufix) ? -t : t); \ } else { \ real_type t = sqrt(2 * FNAME(abs, func_sufix)(z) + fabs(x)); \ real_type u = t / 2; \ return x > CONCAT(0.0, func_sufix) ? (complex_type)(u, y / t) \ : (complex_type)(fabs(y) / t, y < CONCAT(0.0, func_sufix) ? -u : u); \ } \ } \ \ complex_type FNAME(sin, func_sufix)(complex_type z) { \ const real_type x = z.x; \ const real_type y = z.y; \ return (complex_type)(sin(x) * cosh(y), cos(x) * sinh(y)); \ } \ \ complex_type FNAME(sinh, func_sufix)(complex_type z) { \ const real_type x = z.x; \ const real_type y = z.y; \ return (complex_type)(sinh(x) * cos(y), cosh(x) * sin(y)); \ } \ \ complex_type FNAME(cos, func_sufix)(complex_type z) { \ const real_type x = z.x; \ const real_type y = z.y; \ return (complex_type)(cos(x) * cosh(y), -sin(x) * sinh(y)); \ } \ \ complex_type FNAME(cosh, func_sufix)(complex_type z) { \ const real_type x = z.x; \ const real_type y = z.y; \ return (complex_type)(cosh(x) * cos(y), sinh(x) * sin(y)); \ } \ \ complex_type FNAME(tan, func_sufix)(complex_type z) { \ return FNAME(div, func_sufix)(FNAME(sin, func_sufix)(z), FNAME(cos, func_sufix)(z)); \ } \ \ complex_type FNAME(tanh, func_sufix)(complex_type z) { \ return FNAME(div, func_sufix)(FNAME(sinh, func_sufix)(z), FNAME(cosh, func_sufix)(z)); \ } \ \ complex_type FNAME(asinh, func_sufix)(complex_type z) { \ complex_type t = (complex_type)((z.x - z.y) * (z.x + z.y) + CONCAT(1.0, func_sufix), CONCAT(2.0, func_sufix) * z.x * z.y); \ t = FNAME(sqrt, func_sufix)(t) + z; \ return FNAME(log, func_sufix)(t); \ } \ \ complex_type FNAME(asin, func_sufix)(complex_type z) { \ complex_type t = (complex_type)(-z.y, z.x); \ t = FNAME(asinh, func_sufix)(t); \ return (complex_type)(t.y, -t.x); \ } \ \ complex_type FNAME(acosh, func_sufix)(complex_type z) { \ return CONCAT(2.0, func_sufix) * \ FNAME(log, func_sufix)(FNAME(sqrt, func_sufix)(CONCAT(0.5, func_sufix) * (z + CONCAT(1.0, func_sufix))) + \ FNAME(sqrt, func_sufix)(CONCAT(0.5, func_sufix) * (z - CONCAT(1.0, func_sufix)))); \ } \ \ complex_type FNAME(acos, func_sufix)(complex_type z) { \ complex_type t = FNAME(asin, func_sufix)(z); \ return (complex_type)(CONCAT(M_PI_2, math_consts_sufix) - t.x, -t.y); \ } \ \ complex_type FNAME(atanh, func_sufix)(complex_type z) { \ const real_type zy2 = z.y * z.y; \ real_type n = CONCAT(1.0, func_sufix) + z.x; \ real_type d = CONCAT(1.0, func_sufix) - z.x; \ n = zy2 + n * n; \ d = zy2 + d * d; \ return (complex_type)(CONCAT(0.25, func_sufix) * (log(n) - log(d)), \ CONCAT(0.5, func_sufix) * \ atan2(CONCAT(2.0, func_sufix) * z.y, CONCAT(1.0, func_sufix) - zy2 - (z.x * z.x))); \ } \ \ complex_type FNAME(atan, func_sufix)(complex_type z) { \ const real_type zx2 = z.x * z.x; \ real_type n = z.y + CONCAT(1.0, func_sufix); \ real_type d = z.y - CONCAT(1.0, func_sufix); \ n = zx2 + n * n; \ d = zx2 + d * d; \ return (complex_type)(CONCAT(0.5, func_sufix) * atan2(CONCAT(2.0, func_sufix) * z.x, CONCAT(1.0, func_sufix) - zx2 - (z.y * z.y)), \ CONCAT(0.25, func_sufix) * (log(n / d))); \ } // float complex typedef float2 cfloat; OPENCL_COMPLEX_MATH_FUNCS(float2, float, f, _F) // double complex #ifdef cl_khr_fp64 # ifdef OPENCL_COMPLEX_MATH_USE_DOUBLE # pragma OPENCL EXTENSION cl_khr_fp64 : enable typedef double2 cdouble; OPENCL_COMPLEX_MATH_FUNCS(double2, double, , ) # endif #endif #undef FNAME #undef CONCAT #endif // OPENCL_COMPLEX_MATH