318 lines
10 KiB
C
318 lines
10 KiB
C
//---------------------------------------------------------------------------//
|
|
// 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
|