Initial

1 files changed, 158 insertions, 0 deletions

diff --git a/attract.c b/attract.c
new file mode 100644
index 0000000..f4cb7b0
--- /dev/null
+++ b/attract.c
@@ -0,0 +1,158 @@
+#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"
+
+const int MAXITER=200;
+const double EPS=1e-1;
+
+static double complex att_func(double complex x){
+	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);
+}
+
+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{
+	const double unitsz;
+	const int width,height;
+};
+
+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);
+			}
+			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(void){
+	struct config cfg={
+		200,
+		500, 500
+	};
+
+	unsigned char img[cfg.width*cfg.height*3];
+
+	int nthreads=sysconf(_SC_NPROCESSORS_ONLN);
+	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;
+		memcpy(&ths[i].cfg,&cfg,sizeof(struct config));
+		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);
+}