#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <complex.h>
#include <unistd.h>
#include <pthread.h>
#include <assert.h>
#include "lodepng.h"

#ifndef CMPLX
#define CMPLX(r,i) ((r) + (i)*I)
#endif

const int MAXITER=200;
const double EPS=1e-1;

static double complex att_func(double complex x,double exponent){
#if 0
	static const int N=9;

	double complex r=1,y=x;
	for(int i=1;i<N;i++){
		r+=y;
		y*=x;
	}
	r+=y;
	return clog(r);
#else
	return clog((1-cpow(x,exponent+1))/(1-x));
#endif
}

struct colour{
	unsigned char r,g,b;
};

__attribute__((unused))
static struct colour make_colour(int index){
	static const struct colour clr_base[]={
		{255,100,100}, {100,255,100}, {100,100,255},
		{200,200,30}, {200,30,200}, {30,200,200},
		{100,100,100}, {200,200,200}
	};
	static const int nbase=sizeof(clr_base)/sizeof(clr_base[0]);

	return clr_base[index%nbase];
}

struct config{
	double unitsz;
	int width,height;
	double exponent;
};

struct thread_info{
	pthread_t th;
	int ystart,yend;
	unsigned char *img;
	bool do_basins;
	struct config cfg;
};

static void* thread_entry(void *arg_vp){
	struct thread_info *th_info=arg_vp;
	int ystart=th_info->ystart;
	int yend=th_info->yend;
	unsigned char *img=th_info->img;
	bool do_basins=th_info->do_basins;
	struct config cfg=th_info->cfg;

	int nbasins=0,bsz=0;
	double complex *basins=NULL;
	if(do_basins){
		bsz=16;
		basins=malloc(bsz*sizeof(double complex));
		assert(basins);
	}

	printf("Thread: y in [%d,%d)\n",ystart,yend);

	for(int iy=ystart;iy<yend;iy++){
		for(int ix=0;ix<cfg.width;ix++){
			const double complex start=CMPLX(
					(double)(ix-cfg.width/2)/cfg.unitsz,
					(double)(iy-cfg.height/2)/cfg.unitsz);
			double complex z=start;
			for(int i=0;i<MAXITER;i++){
				z=att_func(z,cfg.exponent);
			}
			if(do_basins){
				int bi;
				for(bi=0;bi<nbasins;bi++){
					if(cabs(z-basins[bi])<EPS)break;
				}
				if(bi==nbasins){
					if(nbasins==bsz){
						bsz*=2;
						basins=realloc(basins,bsz*sizeof(double complex));
						assert(basins);
					}
					fprintf(stderr,"New basin at %lf,%lf (from point %lf,%lf)\n",creal(z),cimag(z),creal(start),cimag(start));
					basins[bi]=z;
					nbasins++;
				}
				basins[bi]=(basins[bi]+z)/2;
				
				// struct colour clr=make_colour(bi);
			}

			unsigned char clrv=(int)(cabs(z)/3*255)%256;
			struct colour clr={clrv,clrv,clrv};
			img[3*(iy*cfg.width+ix)+0]=clr.r;
			img[3*(iy*cfg.width+ix)+1]=clr.g;
			img[3*(iy*cfg.width+ix)+2]=clr.b;
		}
	}

	pthread_exit(NULL);
}

int main(int argc,char **argv){
	struct config cfg={
		.unitsz=200,
		.width=500, .height=500,
		.exponent=8.5
	};

	assert(argc==1||argc==2);
	if(argc==2){
		cfg.exponent=strtod(argv[1],NULL);
	}

	unsigned char img[cfg.width*cfg.height*3];

	int nthreads=sysconf(_SC_NPROCESSORS_ONLN);
	if(nthreads==4)nthreads/=2;
	if(nthreads>=999||nthreads<1)nthreads=1;
	fprintf(stderr,"Using %d thread%s\n",nthreads,nthreads==1?"":"s");


	pthread_attr_t attr;
	if(pthread_attr_init(&attr)!=0){
		perror("pthread_attr_init");
		return 1;
	}

	struct thread_info ths[nthreads];
	for(int i=0;i<nthreads;i++){
		ths[i].ystart=i*cfg.height/nthreads;
		ths[i].yend=(i+1)*cfg.height/nthreads;
		ths[i].img=img;
		ths[i].do_basins=nthreads==1;
		ths[i].cfg=cfg;
		int ret=pthread_create(&ths[i].th,&attr,thread_entry,&ths[i]);
		if(ret!=0){
			perror("pthread_create");
			return 1;
		}
	}

	pthread_attr_destroy(&attr);

	for(int i=0;i<nthreads;i++){
		int ret=pthread_join(ths[i].th,NULL);
		if(ret!=0){
			perror("pthread_join");
			return 1;
		}
	}

	lodepng_encode_file("out.png",img,cfg.width,cfg.height,LCT_RGB,8);
}