fractals.cpp

Generating Fractals with SSE/SSE2 You probably have heard about fractals before. They are beautiful

// (c) Kankowski Peter, 2004-2005. [e-mail]

#include "tech.h"
#define WM_MOUSEWHEEL	0x20A
#define TIMER_ZOOMIN  100
#define TIMER_ZOOMOUT 101

LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wParam,LPARAM lParam);

HWND hWnd = NULL, hDlg = NULL, hToolbar = NULL;
const char classname[] = "kankowski_fractals";
WNDCLASS wc = {CS_HREDRAW|CS_VREDRAW |CS_DBLCLKS, (WNDPROC)WndProc, 0, 0,
  NULL, NULL, NULL, NULL, NULL, classname};
WNDPROC OldProc; HACCEL hAccel;

BITMAPINFOHEADER bi = {sizeof(BITMAPINFOHEADER), 0, 0, 1, 32, BI_RGB, 0, 0, 0, 0, 0};
BITMAPINFOHEADER bis = {sizeof(BITMAPINFOHEADER), 0, 0, 1, 32, BI_RGB, 0, 0, 0, 0, 0};

HBITMAP hBMP, hsBMP;
long *bits=0, *sbits;
char buff[4096];
const double LEFTm = -2.375, RIGHTm = 0.775, TOPm = -1.1, BOTTOMm = 1.1;
const double LEFTj = -1.575, RIGHTj = 1.575, TOPj = -1.1, BOTTOMj = 1.1;
//const double LEFTm = -0.84, RIGHTm = -0.79, TOPm = 0.18, BOTTOMm = 0.21;
//const double LEFTj = -0.897, RIGHTj = -0.894, TOPj = 0.673, BOTTOMj = 0.675;
double LEFT = LEFTj, RIGHT = RIGHTj,
       TOP = TOPj, BOTTOM = BOTTOMj;
#define zoomspeed 400
int SSEmode=1, SSEsupported=0, SSE2supported=0, ms=zoomspeed;
#define ITER 64
COLORREF a[ITER+1];
LARGE_INTEGER freq;
DWORD bz; RECT clientrect = {0, 0, 0, 0};
HBITMAP holdBMP; HDC hmem;
LPDIRECTDRAW DD = 0;
LPDIRECTDRAWSURFACE PrimarySurface = 0;
LPDIRECTDRAWSURFACE OffscreenSurface = 0;
LPDIRECTDRAWCLIPPER Clipper = 0;
WPARAM WindowActive;
int DDsupported = 1, DDmode = 1, movadd = 1;

typedef HRESULT (WINAPI *DIRECTDRAWCREATE)(GUID*, LPDIRECTDRAW*, IUnknown*);
HMODULE ddraw = 0;
struct MYRGNDATA {
	RGNDATAHEADER rdh;
	RECT r[4];
};

extern "C" void PaintJuliaMandelSSE(int w, int h, float dx, float dy, long* bits,
					COLORREF* a, float left, float top, float px, float py, int piw,
					int type, int bw);
extern "C" void PaintJuliaMandelSSE2(int w, int h, double dx, double dy, long* bits,
					COLORREF* a, double left, double top, double px, double py, int piw,
					int type, int bw);
extern "C" int IsSSE();
extern "C" int IsSSE2();
extern "C" int IsMovAdd();


#define ID_ARROW 100
#define ID_ZOOM 101
#define ID_PIP 102
#define ID_MANDEL 103
#define ID_JULIA 104
#define ID_BW 105
#define ID_AXIS 106
#define ID_SSE 107
#define ID_HOME 108
#define IDC_TOOLBAR 666
char* tooltips[] = {"Table and graph of iterates for the selected point (F2)",
	"Zoom (F3)", "Picture in picture: correspondence between Mandelbrot and Julia sets (F4)",
	"Mandelbrot set (Space to toggle)", "Julia set (Space to toggle)",
	"Color on/off (F5)", "Axes on/off (F6)", "Drawing mode (F7)", "Back to initial size (F8)"};
ACCEL accel[] = {
	{FNOINVERT | FVIRTKEY, VK_F2, ID_ARROW},
	{FNOINVERT | FVIRTKEY, VK_F3, ID_ZOOM},
	{FNOINVERT | FVIRTKEY, VK_F4, ID_PIP},

	{FNOINVERT | FVIRTKEY, VK_F5, ID_BW},
	{FNOINVERT | FVIRTKEY, VK_F6, ID_AXIS},
	{FNOINVERT | FVIRTKEY, VK_F7, ID_SSE},
	{FNOINVERT | FVIRTKEY, VK_F8, ID_HOME}
};
int mode=ID_ARROW, type=ID_JULIA, bw = 0, axis = 0, fired, oldtype, beep=0;
HMENU menu, popup;
double PX = -0.12, PY = 0.74, X0 = 0, Y0 = 0;

const TBBUTTON toolbar[] = {
	{0, ID_ARROW, TBSTATE_ENABLED | TBSTATE_CHECKED, TBSTYLE_BUTTON, 0},
	{1, ID_ZOOM, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0},	
	{2, ID_PIP, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0},
	{-1, 0, 0, TBSTYLE_SEP, 0},
	{3, ID_JULIA, TBSTATE_ENABLED | TBSTATE_CHECKED, TBSTYLE_BUTTON, 0},
	{4, ID_MANDEL, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0},
	{-1, 0, 0, TBSTYLE_SEP, 0},
	{5, ID_BW, TBSTATE_ENABLED | TBSTATE_CHECKED, TBSTYLE_BUTTON, 0},
	{6, ID_AXIS, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0},
	{-1, 0, 0, TBSTYLE_SEP, 0},
	{7, ID_SSE, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0},
	{8, ID_HOME, TBSTATE_ENABLED, TBSTYLE_BUTTON, 0}
};

void PaintMandelbrotFPU(int w, int h, double dx, double dy, long* b, double left, double first, int piw) {
	int x, y, i; double x2, y2, px, py, zx, zy;
	py = first;
	if(bw) {
		for(y = 0; y < h; ++y) {
			px = left;
			for(x = 0; x < w; ++x) {
				zx = px, zy = py;
				for(i = 0; i < ITER; ++i) {
					x2 = zx * zx,  y2 = zy * zy;
					if(x2 + y2 > 4.0) {
						*b++ = 0x808080;
						goto NEXT;
					}
					zy = zx * zy * 2 + py;
					zx = x2 - y2 + px;
				}
				*b++ = 0;
NEXT:			px += dx;
			}
			b += piw / 4;
			py += dy;
		}
	}
	else {
		for(y = 0; y < h; ++y) {
			px = left;
			for(x = 0; x < w; ++x) {
				zx = px, zy = py;
				for(i = 0; i < ITER; ++i) {
					x2 = zx * zx,   y2 = zy * zy;
					if(x2 + y2 > 4.0)
						break;
					zy = zx * zy * 2 + py;
					zx = x2 - y2 + px;
				}
				*b++ = a[i];
				px += dx;
			}
			b += piw / 4;
			py += dy;
		}
	}
}

void PaintJuliaFPU(int w, int h, double dx, double dy, long* b, double left, double first, int piw) {
	int x, y, i; double x2, y2, zx2, zy2, zx, zy;
	zy2 = first;
	if(bw) {
		for(y = 0; y < h; ++y) {
			zx2 = left;
			for(x = 0; x < w; ++x) {
				zx = zx2, zy = zy2;
				for(i = 0; i < ITER; ++i) {
					x2 = zx * zx,   y2 = zy * zy;
					if(x2 + y2 > 4.0) {
						*b++ = 0x808080;
						goto NEXT;
					}
					zy = zx * zy * 2 + PY;
					zx = x2 - y2 + PX;
				}
				*b++ = 0;
NEXT:			zx2 += dx;
			}
			b += piw / 4;
			zy2 += dy;
		}
	}
	else {
		for(y = 0; y < h; ++y) {
			zx2 = left;
			for(x = 0; x < w; ++x) {
				zx = zx2, zy = zy2;
				for(i = 0; i < ITER; ++i) {
					x2 = zx * zx,   y2 = zy * zy;
					if(x2 + y2 > 4.0)
						break;
					zy = zx * zy * 2 + PY;
					zx = x2 - y2 + PX;
				}
				*b++ = a[i];
				zx2 += dx;
			}
			b += piw / 4;
			zy2 += dy;
		}
	}
}

double atod(char*& p) { // String to double
	unsigned sign = 0, v; double d = 0.0, factor = 0.1; int exp = 0; char *a = p;
	while(' ' == *a) // Skip leading spaces
		a++;
	if('+' == *a) a++;
	if('-' == *a) a++, sign = 1;
	while((v = *a - '0') <= 9) // Integer part
		a++, d = 10.0 * d + v;
	if('.' == *a || ',' == *a) { // Fractional part
		a++;
		while((v = *a - '0') <= 9)
			a++, d += v * factor, factor *= 0.1;
	}
	if((*a & 0xDF) == 'E') { // Exponent
		a++;
		if('+' == *a) a++, factor = 10;
		if('-' == *a) a++, factor = 0.1;
		while((v = *a - '0') <= 9)
			a++, exp = 10 * exp + v;
		while(exp) {
			if(exp & 1)
				d *= factor;
			exp >>= 1;
			factor *= factor; // factor = factor ^ 2
		}
	}
	p = a;
	return sign ? -d : d;
}

void __forceinline atox(char* a, double& x, double&y) { // String to complex number
	while(' ' == *a) // Skip spaces
		a++;
	if('i' == *a) { // i without number = 1 * i
		y = 1, x = 0;
		return;
	}
	x = atod(a); // Real part
	while(' ' == *a) // Skip spaces
		a++;
	if('i' == *a) { // Imaginary part goes first
		a++, y = x, x = atod(a);
		if('i' == *a)
			x = 0;
		return;
	}
	y = atod(a); // Imaginary part
	if('i' != *a)
		y = 0;
}

double ceil(double x) {
	__asm {
		FLD x
		FRNDINT
		// FISTP QWORD PTR [TEMP]
		// FILD QWORD PTR [TEMP]
	}
}

unsigned inline dtou(double x) {
	int y;
	__asm {
		FLD x
		FISTP y
		MOV eax, DWORD PTR y
	}
}

int inline dtoi(double x) {
	int y;
	__asm {
		FLD x
		FISTP y
		MOV eax, DWORD PTR y
	}
}

void dtoa(double x, char*& a) { // Double to string
	unsigned v, i; char *p = a;
	if(x < 0)
		x = -x, *p++ = '-';
	v = dtou(x); // Get integer part to v
	do {
		*p++ = char('0' + (v % 10)); // Output integer part
		v /= 10;
	} while(v);
	x = x - ceil(x); // Leave only fractional part
	if(x) {
		x *= 10;
		*p++ = '.';
		for(i = 7; i && x; i--)
			v = dtou(x), *p++ = char('0' + v), x = x - ceil(x), x *= 10;
	}
	*p = '\0';
	a = p;
}

char* xtoa(char* a, double x, double y) {
	dtoa(x, a);
	if(y >= 0)
		*a++ = '+';
	dtoa(y, a);
	set2ch(a, 'i', '\0');
	return a+1;
}

inline char* xxxtoa(int x, char* a) { // Iteration number
	*a++ = 'Z';
	if(x / 10) // If the is the second figure
		*a++ = char(x / 10 + '0');
	*a++ = char(x % 10 + '0'); // The first figure
	*a++ = '=';
	return a;
}

void EnableDisable(WPARAM id1, WPARAM id2, WPARAM id3) {
	SendMessage(hToolbar, TB_SETSTATE, id1, TBSTATE_CHECKED | TBSTATE_ENABLED);
	SendMessage(hToolbar, TB_SETSTATE, id2, TBSTATE_ENABLED);
	SendMessage(hToolbar, TB_SETSTATE, id3, TBSTATE_ENABLED);
}
void CheckUncheck(int id1, int id2, int id3) {
	CheckMenuItem(popup, id1, MF_BYCOMMAND | MF_CHECKED);
	CheckMenuItem(popup, id2, MF_BYCOMMAND | MF_UNCHECKED);
	CheckMenuItem(popup, id3, MF_BYCOMMAND | MF_UNCHECKED);
}

double flp2(double x) { // 侮痼汶栩