Skip to content

Instantly share code, notes, and snippets.

@knotman90 knotman90/cl_complex.cl
Last active Feb 3, 2017

Embed
What would you like to do?
/**
Author: Davide Spataro
email:davide.spataro@unical.it
www.davidespataro.it
*/
#define DOUBLE_PRECISION (1)
typedef double2 cl_double_complex;
typedef float2 cl_float_complex;
#ifdef DOUBLE_PRECISION
typedef cl_double_complex cl_complex;
typedef double TYPE;
#else
typedef cl_float_complex cl_complex;
typedef float TYPE;
#endif
inline TYPE cl_complex_real_part(const cl_complex* n){
return n->x;
}
inline TYPE cl_complex_imaginary_part(const cl_complex* n){
return n->y;
}
/*
* Returns modulus of complex number (its length):
*/
inline TYPE cl_complex_modulus(const cl_complex* n){
return (sqrt( (n->x*n->x)+(n->y*n->y) ));
}
inline cl_complex cl_complex_add(const cl_complex* a, const cl_complex* b){
return (cl_complex)( a->x + b->x, a->y + b->y );
}
inline cl_complex cl_complex_multiply(const cl_complex* a, const cl_complex* b){
return (cl_complex)(a->x*b->x - a->y*b->y, a->x*b->y + a->y*b->x);
}
/*
* Computes the integer Integer power of a complex number.
*
* */
inline cl_complex cl_complex_ipow(const cl_complex* base , int exp ){
cl_complex res;
res.x=res.y=1;
while(exp){
if(exp & 1)
res=cl_complex_multiply(&res,base);
exp>>=1;
res = cl_complex_multiply(&res,&res);
}
return res;
}
inline cl_complex cl_complex_divide(const cl_complex* a, const cl_complex* b){
const TYPE dividend = (b->x*b->x + b->y*b->y);
return (cl_complex)((a->x*b->x + a->y*b->y)/dividend , (a->y*b->x - a->x*b->y)/dividend);
}
/*
* Get the argument of a complex number (its angle):
* http://en.wikipedia.org/wiki/Complex_number#Absolute_value_and_argument
*/
inline TYPE cl_complex_argument(const cl_complex* a){
if(a->x > 0){
return atan(a->y / a->x);
}else if(a->x < 0 && a->y >= 0){
return atan(a->y / a->x) + M_PI;
}else if(a->x < 0 && a->y < 0){
return atan(a->y / a->x) - M_PI;
}else if(a->x == 0 && a->y > 0){
return M_PI/2;
}else if(a->x == 0 && a->y < 0){
return -M_PI/2;
}else{
return 0;
}
}
/*
* Returns the Square root of complex number.
* Although a complex number has two square roots,
* only one of them -the principal square root- is computed.
* see wikipedia:http://en.wikipedia.org/wiki/Square_root#Principal_square_root_of_a_complex_number
*/
inline cl_complex cl_complex_sqrt(const cl_complex* n){
const TYPE sm = sqrt(cl_complex_modulus(n));
const TYPE a2 = cl_complex_argument(n)/2;
const TYPE ca2 = cos(a2);
const TYPE sa2 = sin(a2);
return (cl_complex)(sm * ca2 , sm * sa2);
}
/*
* Computes the exponential function for the complex number z = x+iy
* as: exp(z) = e^x(cos(y) + isin(y))
* See: https://en.wikipedia.org/wiki/Exponential_function#Complex_plane
* */
inline cl_complex_exp(const cl_complex* n){
const TYPE e = exp(n->x);
return (cl_complex)(e*cos(n->y) , e*sin(n->y));
}
/*
* Computes the logatirm of the complex number z= x+iy
* x+iy = re^{i\theta}
* log(x+iy) = log(re^{i\theta} = log(r)+log(e^{i\theta}
* where r is the module of |z| = sqrt(x^2+y^2)
* log(z) = log(|z|) + iarg(z) = log(sqrt(x^2+y^2) + i atan(y/b)
*
* */
inline cl_complex_log(const cl_complex* z){
return (cl_complex)( log(cl_complex_modulus(z)) , cl_complex_argument(z));
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.