cl_complex.cl

/**
Author: Davide Spataro
email:[email protected]
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));
}

Hi. I have seen your video. looks nice. I was wandering how c i changing in this video so made this program:

#include <stdio.h>      /* Standard Library of Input and Output */
#include <complex.h>    /* Standard Library of Complex Numbers */


/*



  gcc j.c -lm -Wall
  ./a.out > j.txt
  gnuplot plot ./j.txt
 
 https://gist.github.com/knotman90/9b4e47dc04e7dbfefe4a
 https://www.youtube.com/watch?v=0kId_t6AqjI
 
  Host function that prepares data array and passes it to the CUDA kernel.
 */
//# c = 1j              # dentrite fractal
//# c = -0.123 + 0.745j # douady's rabbit fractal
//# c = -0.750 + 0j     # san marco fractal
//# c = -0.391 - 0.587j # siegel disk fractal
//# c = -0.7 - 0.3j     # NEAT cauliflower thingy
//# c = -0.75 - 0.2j    # galaxies
//# c = -0.75 + 0.15j   # groovy
//# c = -0.7 + 0.35j    # frost
int main(void) {

	
	
	double incre =0.00000003;
	double inci =-0.00009;
	double startre=-0.75;
	double starti=0.09;
	//double zoom=1.0;
	complex double c;
	
	
	printf("# x y \n"); // header for txt file with data for gnuplot
	
	
	for(int i=0;i<2500;i++){
		//zoom+=0.001;
		
		c =  startre+i*incre + (starti+i*inci)*I ;
		// 		sprintf(fileName, "julia_%f_%f.png", startre+i*incre,starti+i*inci);
		printf("%.16f %+.16f \n", creal(c), cimag(c));
		
	}



	return 0;
}

Is this OK ?