From b66b6aee38fee21e3789a249e14e32a25f696efe Mon Sep 17 00:00:00 2001 From: josh Date: Thu, 29 Nov 2007 03:48:43 +0000 Subject: [PATCH] initial commit, pthread support missing git-svn-id: svn://anubis/misc/newton@17 bd8a9e45-a331-0410-811e-c64571078777 --- Makefile | 17 +++++ complex.c | 69 +++++++++++++++++ complex.h | 25 +++++++ newton.c | 217 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 328 insertions(+) create mode 100644 Makefile create mode 100644 complex.c create mode 100644 complex.h create mode 100644 newton.c diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..7c2969c --- /dev/null +++ b/Makefile @@ -0,0 +1,17 @@ + +CC := gcc +TARGET := newton +OBJS := complex.o newton.o +CFLAGS := -O2 `sdl-config --cflags` +LDFLAGS := `sdl-config --libs` + +all: $(TARGET) + +$(TARGET): $(OBJS) + $(CC) -o $@ $^ $(LDFLAGS) + +%.o: %.c + $(CC) -o $@ -c $^ $(CFLAGS) + +clean: + -rm -f $(TARGET) $(OBJS) *~ diff --git a/complex.c b/complex.c new file mode 100644 index 0000000..dc3376d --- /dev/null +++ b/complex.c @@ -0,0 +1,69 @@ +/* Josh Holtrop + * 10/12/05 + * Complex number functions + */ + +#include "complex.h" + +void complex_add(complex_t *c1, complex_t *c2, complex_t *result) +{ + result->a = c1->a + c2->a; + result->b = c1->b + c2->b; +} + +void complex_adds(complex_t *c1, double scalar, complex_t *result) +{ + result->a = c1->a + scalar; + result->b = c1->b; +} + +void complex_sub(complex_t *c1, complex_t *c2, complex_t *result) +{ + result->a = c1->a - c2->a; + result->b = c1->b - c2->b; +} + +void complex_subs(complex_t *c1, double scalar, complex_t *result) +{ + result->a = c1->a - scalar; + result->b = c1->b; +} + +void complex_mul(complex_t *c1, complex_t *c2, complex_t *result) +{ + double a = c1->a * c2->a - c1->b * c2->b; + result->b = c1->a * c2->b + c2->a * c1->b; + result->a = a; +} + +void complex_div(complex_t *c1, complex_t *c2, complex_t *result) +{ + double a = c1->a, b = c1->b, c = c2->a, d = c2->b; + double ra = (a * c + b * d) / (c * c + d * d); + result->b = (c * b - a * d) / (c * c + d * d); + result->a = ra; +} + +void complex_muls(complex_t *c, double scalar, complex_t *result) +{ + result->a = c->a * scalar; + result->b = c->b * scalar; +} + +void complex_divs(complex_t *c, double scalar, complex_t *result) +{ + result->a = c->a / scalar; + result->b = c->b / scalar; +} + +void complex_pow(complex_t *c, int n, complex_t *result) +{ + /* Doesn't handle c^0 */ + result->a = c->a; + result->b = c->b; + for ( ; n > 1; n--) + { + complex_mul(result, c, result); + } +} + diff --git a/complex.h b/complex.h new file mode 100644 index 0000000..7ea971e --- /dev/null +++ b/complex.h @@ -0,0 +1,25 @@ +/* Josh Holtrop + * 10/12/05 + * Complex number functions + */ + +#ifndef __COMPLEX_H__ +#define __COMPLEX_H__ __COMPLEX_H__ + +typedef struct +{ + double a; + double b; +} complex_t; + +void complex_add(complex_t *c1, complex_t *c2, complex_t *result); +void complex_adds(complex_t *c1, double scalar, complex_t *result); +void complex_sub(complex_t *c1, complex_t *c2, complex_t *result); +void complex_subs(complex_t *c1, double scalar, complex_t *result); +void complex_mul(complex_t *c1, complex_t *c2, complex_t *result); +void complex_div(complex_t *c1, complex_t *c2, complex_t *result); +void complex_muls(complex_t *c, double scalar, complex_t *result); +void complex_divs(complex_t *c, double scalar, complex_t *result); +void complex_pow(complex_t *c, int n, complex_t *result); + +#endif diff --git a/newton.c b/newton.c new file mode 100644 index 0000000..b42644a --- /dev/null +++ b/newton.c @@ -0,0 +1,217 @@ +/* + * Josh Holtrop + * 2007-11-28 + * Newton Fractal Renderer using MPI and SDL + */ + +#include +#include +#include +#include "complex.h" + +#define LEEWAY 0.001 +#define MAX_ITERATIONS 32 +#define DEFAULT_WIN_WIDTH 800 +#define DEFAULT_WIN_HEIGHT 600 +#define PROGNAME "Josh's MPI/pthread newton fractal renderer" + +/* Prototypes */ +SDL_Surface * sdl_init(unsigned int width, unsigned int height); +void calculateRow(int winWidth, int winHeight, int row, double viewWidth, + double viewHeight, double x_center, double y_center, Uint32 * rowVals); +Uint32 computePoint(double x, double y); +void draw(SDL_Surface * screen, int winWidth, int winHeight, + double x_center, double y_center, double viewWidth, double viewHeight); +void main_loop(SDL_Surface * screen, int winWidth, int winHeight, + double x_center, double y_center, double viewWidth, double viewHeight); + +int main(int argc, char * argv[]) +{ + SDL_Surface * screen; + int winWidth = DEFAULT_WIN_WIDTH; + int winHeight = DEFAULT_WIN_HEIGHT; + double viewWidth = 2.0, x_center = 0.0, y_center = 0.0; + double viewHeight = ((double)winHeight/(double)winWidth) * viewWidth; + + if ( !(screen = sdl_init(winWidth, winHeight)) ) + { + fprintf(stderr, "SDL initialization error!\n"); + return -1; + } + + main_loop(screen, winWidth, winHeight, x_center, y_center, + viewWidth, viewHeight); + + return 0; +} + +void main_loop(SDL_Surface * screen, int winWidth, int winHeight, + double x_center, double y_center, double viewWidth, double viewHeight) +{ + SDL_Event event; + for (;;) + { + draw(screen, winWidth, winHeight, x_center, y_center, + viewWidth, viewHeight); + while(SDL_WaitEvent(&event)) + if (event.type == SDL_QUIT || + event.type == SDL_MOUSEBUTTONDOWN) + break; + switch (event.type) + { + case SDL_QUIT: + return; + case SDL_MOUSEBUTTONDOWN: + { + int button = event.button.button; + int x = event.motion.x; + int y = event.motion.y; + switch (button) + { + case 1: /* re-center on click point, zoom in */ + x_center += (viewWidth / winWidth) * (x - (winWidth >> 1)); + y_center += (viewHeight / winHeight) * ((winHeight >> 1) - y); + viewWidth /= 2.0; + viewHeight /= 2.0; + break; + case 2: /* rectangular selection to zoom into */ + { + int ox = x, oy = y; + while (SDL_WaitEvent(&event)) + if (event.type == SDL_MOUSEBUTTONUP) + break; + x = event.motion.x; + y = event.motion.y; + if (x < ox) + { int t = ox; ox = x; x = ox; } + if (y < oy) + { int t = oy; oy = y; y = oy; } + x_center += (viewWidth / winWidth) * + (ox + ((x-ox) >> 1) - (winWidth >> 1)); + y_center += (viewHeight / winHeight) * + ((winHeight >> 1) - (oy + ((y-oy) >> 1))); + viewWidth = (viewWidth / winWidth) * (x-ox+1); + viewHeight = (viewHeight / winHeight) * (y-oy+1); + } + break; + case 3: /* reset view */ + viewWidth = 2.0; x_center = 0.0; y_center = 0.0; + viewHeight = ((double)winHeight/(double)winWidth) * viewWidth; + break; + case 4: /* zoom in */ + viewWidth /= 2.0; + viewHeight /= 2.0; + break; + case 5: /* zoom out */ + viewWidth *= 2.0; + viewHeight *= 2.0; + break; + } + } + } + } +} + +void draw(SDL_Surface * screen, int winWidth, int winHeight, + double x_center, double y_center, double viewWidth, double viewHeight) +{ + Uint32 * pixels = (Uint32 *) screen->pixels; + int ix, iy; + for (iy = 0; iy < winHeight; iy++) + { + calculateRow(winWidth, winHeight, iy, viewWidth, viewHeight, + x_center, y_center, pixels); + pixels += winWidth; + } + SDL_Flip(screen); +} + +SDL_Surface * sdl_init(unsigned int width, unsigned int height) +{ + SDL_Surface * screen; + + if (SDL_Init(SDL_INIT_VIDEO)) + { + fprintf(stderr, "Failed to initialize SDL!\n"); + return NULL; + } + + atexit(SDL_Quit); + + if (!(screen = SDL_SetVideoMode(width, height, 32, + SDL_DOUBLEBUF | SDL_HWSURFACE))) + { + fprintf(stderr, "Failed to set video mode!\n"); + return NULL; + } + + SDL_WM_SetCaption(PROGNAME, PROGNAME); + + return screen; +} + +void calculateRow(int winWidth, int winHeight, int row, double viewWidth, + double viewHeight, double x_center, double y_center, Uint32 * rowVals) +{ + int i; + double xspacing = viewWidth / winWidth; + double yspacing = viewHeight / winHeight; + double y = ((winHeight >> 1) - row) * yspacing + y_center; + double x = x_center - (viewWidth / 2); + for (i = 0; i < winWidth; i++, x += xspacing) + *rowVals++ = computePoint(x, y); +} + +Uint32 computePoint(double x, double y) +{ + int n; + complex_t rootGuess = {y, x}; + complex_t t1, t2; + for (n = 0; n < MAX_ITERATIONS; n++) + { + /* percentage of color to illuminate based on current iteration */ + double pIlum = 1.0 - (n * 0.03); + /* + * These if statements check to see if the complex number (the root + * guess) is within LEEWAY distance of a real root. If so, a unique + * color is returned reflecting which root it is close to and how many + * iterations it took for the root guess to get to that root + */ + if ((fabs(rootGuess.a - 1) < LEEWAY) + && (fabs(rootGuess.b) < LEEWAY)) + return ((Uint32)(0xFF * pIlum)) << 16; + if ((fabs(rootGuess.a + 1) < LEEWAY) + && (fabs(rootGuess.b) < LEEWAY)) + return ((Uint32)(0xFF * pIlum) << 16) + ((Uint32)(0xFF * pIlum) << 8); + if ((fabs(rootGuess.a - .5) < LEEWAY) + && (fabs(rootGuess.b - sqrt(3)/2) < LEEWAY)) + return (Uint32)(0xFF * pIlum) << 8; + if ((fabs(rootGuess.a + .5) < LEEWAY) + && (fabs(rootGuess.b - sqrt(3)/2) < LEEWAY)) + return (Uint32)(0x88 * pIlum); + if ((fabs(rootGuess.a - .5) < LEEWAY) + && (fabs(rootGuess.b + sqrt(3)/2) < LEEWAY)) + return (Uint32)(0xFF * pIlum); + if ((fabs(rootGuess.a + .5) < LEEWAY) + && (fabs(rootGuess.b + sqrt(3)/2) < LEEWAY)) + return ((Uint32)(0xFF * pIlum) << 16) + (Uint32)(0xFF * pIlum); + + /* This expression evaluates the next complex number to be tested + * for being close to a root. It uses Newton's method for finding + * roots of equations according to the following recursive equation: + * x_n+1 = x_n - y_n / dy_n + * --> x_n+1 = x_n - (x_n^6 - 1) / 6x_n^5 + * More information can be found at: + * http://www.willamette.edu/~sekino/fractal/fractal.htm + * http://www.math.iastate.edu/danwell/Fexplain/newt1.html + * and Thomas' CALCULUS 10th edition by Finney, Weir, & Giordano pg. 302 + */ + complex_pow(&rootGuess, 6, &t1); + complex_subs(&t1, 1, &t1); + complex_pow(&rootGuess, 5, &t2); + complex_muls(&t2, 6, &t2); + complex_div(&t1, &t2, &t1); + complex_sub(&rootGuess, &t1, &rootGuess); + } + return 0; +}