1 /* This file was taken from XaoS-the fast portable realtime interactive
2    fractal zoomer. but it is simplified for BB. You may get complette
3    sources at XaoS homepage (http://www.paru.cas.cz/~hubicka/XaoS
4  */
5 #ifdef _plan9_
6 #include <u.h>
7 #include <stdio.h>
8 #include <libc.h>
9 #else
10 #include <math.h>
11 #include <string.h>
12 #endif
13 #include "zoom.h"
14 #include "complex.h"
15 /*This code is really ugly-unstrictalized, unextensible, works just
16    for rectangular views. I didn't discovered better way how to get
17    higher speed. For example in case I will implement zre zim iterating
18    function as seperate static one should drastically slow down on my
19    hp425e. Also gotos improves speed-does not force useless if..
20    I can not use normal calculating routines, because lots of thinks
21    must be done in internal loop so internal loop must be inside
22    main loop, because shares lots of variables otherwise register
23    alocation is imposible
24  */
25 
26 #undef SAG			/*solid anti-guessing */
27 #define NOT_CALCULATED 0
28 #define INSET 254
29 #define INPROCESS 255
30 #define RMIN -range
31 #define RMAX range
32 #define IMIN -range
33 #define IMAX range
34 #define QMAX 1000
35 #ifdef STATISTICS
36 int iters2, guessed2, unguessed2, total2, frames2;
37 #endif
do_julia(zoom_context * context,register number_t pre,register number_t pim)38 void do_julia(zoom_context * context, register number_t pre, register number_t pim)
39 {
40     int i, i1, i2, j, x, y, scanline = context->scanline;
41     int iter;
42     number_t rp = 0, ip = 0 /*, t */ ;
43     register number_t zre, zim, im, xdelta, ydelta, br, tmp, range,
44      rangep;
45     number_t xstep, ystep;
46     unsigned char *queue[QMAX];
47     unsigned char **qptr;
48     unsigned char *addr, *addr1 = context->vbuff;
49     const int mode = context->currentformula->juliamode;
50 #ifdef STATISTICS
51     int guessed = 0, unguessed = 0, iters = 0;
52     frames2++;
53 #endif
54     if (mode != 6)
55 	range = 2.0;
56     else
57 	range = 4.0;
58     rangep = range * range;
59 
60     for (i = 0; i < context->height; i++) {
61 	addr = addr1 + i * context->scanline;
62 	memset((char *) addr, NOT_CALCULATED, context->width);
63     }
64     xdelta = context->width / (RMAX - RMIN);
65     ydelta = context->height / (IMAX - IMIN);
66     xstep = (RMAX - RMIN) / context->width;
67     ystep = (IMAX - IMIN) / context->height;
68     for (i = 0; i < context->height; i++) {
69 	im = IMIN + (i + 0.5) * ystep;
70 	x = sqrt(rangep - im * im) * xdelta + 0.5;
71 	addr = addr1 + i * scanline;
72 	y = context->width / 2 - x;
73 	for (j = 0; j < y; j++) {
74 	    addr[j] = 1;
75 	}
76 	y = context->width;
77 	for (j = context->width / 2 + x; j < y; j++) {
78 	    addr[j] = 1;
79 	}
80     }
81     for (i2 = 0; i2 < 2; i2++)
82 	for (i1 = 0; i1 < context->height; i1++) {
83 	    if (i1 % 2)
84 		i = context->height / 2 - i1 / 2;
85 	    else
86 		i = context->height / 2 + i1 / 2 + 1;
87 	    if (i >= context->height)
88 		continue;
89 	    im = IMIN + (i + 0.5) * ystep;
90 	    for (j = (i + i2) & 1; j < context->width; j += 2) {
91 		STAT(total2++;
92 		    )
93 		    addr = addr1 + j + i * scanline;
94 		if (*addr != NOT_CALCULATED)
95 		    continue;
96 		x = j;
97 		y = i;
98 		if (x > 0 && x < context->width &&
99 		    *(addr + 1) != NOT_CALCULATED && *(addr + 1) == *(addr - 1) && y > 0 && y < context->height && *(addr + 1) ==
100 		*(addr + scanline) && *(addr + 1) == *(addr - scanline)) {
101 		    *addr = *(addr + 1);
102 		    continue;
103 		}
104 		zim = im;
105 		zre = RMIN + (j + 0.5) * xstep;
106 		iter = 0;
107 		qptr = queue;
108 		ip = (zim * zim);
109 		rp = (zre * zre);
110 		if (rp + ip > rangep)
111 		    goto outset;
112 		while (1) {
113 		    if (*addr != NOT_CALCULATED
114 #ifdef SAG
115 			&& (*addr == INPROCESS ||
116 			    (!(x > 0 && x < context->width
117 			       && y > 0 && y < context->height &&
118 			       *(addr + 1) != NOT_CALCULATED && (
119 								    (*(addr + 1) != *(addr - 1) && *(addr - 1) != NOT_CALCULATED) ||
120 								    (*(addr + 1) != *(addr + scanline) && *(addr + scanline) != NOT_CALCULATED) ||
121 								    (*(addr + 1) != *(addr - scanline) && *(addr - scanline) != NOT_CALCULATED)))))
122 #endif
123 			) {
124 			if (*addr == INPROCESS || *addr == INSET) {
125 			    *qptr = addr;
126 			    qptr++;
127 			    STAT(guessed++;
128 				)
129 				goto inset;
130 			}
131 			STAT(guessed++;
132 			    )
133 			    iter = *addr;
134 			goto outset;
135 		    }
136 #ifdef STATISTICS
137 		    if (*addr != NOT_CALCULATED)
138 			unguessed++;
139 #endif
140 		    if (*addr != INPROCESS) {
141 			*qptr = addr;
142 			qptr++;
143 			*addr = INPROCESS;
144 			if (qptr >= queue + QMAX)
145 			    goto inset;
146 		    }
147 		    STAT(iters++;
148 			)
149 			switch (mode) {		/*technique called cut and paste */
150 		    case 0:
151 			zim = (zim * zre) * 2 + pim;
152 			zre = rp - ip + pre;
153 			break;
154 		    case 1:
155 			rp = zre * (rp - 3 * ip);
156 			zim = zim * (3 * zre * zre - ip) + pim;
157 			zre = rp + pre;
158 			break;
159 #ifdef _RED_
160 		    case 2:
161 			rp = rp * rp - 6 * rp * ip + ip * ip + pre;
162 			zim = 4 * zre * zre * zre * zim - 4 * zre * ip * zim + pim;
163 			zre = rp;
164 			break;
165 		    case 3:
166 			c_pow4(zre, zim, rp, ip);
167 			c_mul(zre, zim, rp, ip, t, zim);
168 			zre = t + pre;
169 			zim += pim;
170 			break;
171 		    case 4:
172 			c_pow3(zre, zim, rp, ip);
173 			c_pow3(rp, ip, t, zim);
174 			zre = t + pre;
175 			zim += pim;
176 			break;
177 		    case 5:
178 			if (zre >= 0) {
179 			    c_mul(zre - 1, zim, pre, pim, rp, zim);
180 			}
181 			else {
182 			    c_mul(zre + 1, zim, pre, pim, rp, zim);
183 			}
184 			zre = rp;
185 			break;
186 		    case 6:
187 			br = zre + zre + pre - 2;
188 			tmp = zre * zim;
189 			zre = rp - ip + pre - 1;
190 			ip = zim + zim + pim;
191 			zim = tmp + tmp + pim;
192 			tmp = br * br + ip * ip;
193 			rp = (zre * br + zim * ip) / tmp;
194 			ip = (zim * br - zre * ip) / tmp;
195 			zre = (rp + ip) * (rp - ip);
196 			zim = rp * ip;
197 			zim += zim;
198 			rp = zre - 1;
199 			if (rp * rp + zim * zim < 0.001)
200 			    goto outset;
201 			break;
202 #endif
203 		    }
204 		    ip = (zim * zim);
205 		    rp = (zre * zre);
206 		    if (rp + ip > 4)
207 			goto outset;
208 		    x = (zre - RMIN) * xdelta;
209 		    y = (zim - IMIN) * ydelta;
210 		    addr = addr1 + x + y * scanline;
211 		    if (x > 0 && x < context->width &&
212 			*(addr + 1) != NOT_CALCULATED && *(addr + 1) == *(addr - 1) && y > 0 && y < context->height && *(addr + 1) ==
213 			*(addr + scanline) && *(addr + 1) == *(addr - scanline)) {
214 			*addr = *(addr + 1);
215 		    }
216 		}
217 		addr = addr1 + x + y * scanline;
218 		if (*addr == INPROCESS || *addr == INSET) {
219 		  inset:
220 		    while (qptr > queue) {
221 			qptr--;
222 			**qptr = INSET;
223 		    }
224 		}
225 		else {
226 		  outset:
227 		    x = context->maxiter;
228 		    y = context->num_colors;
229 		    while (qptr > queue) {
230 			qptr--;
231 			**qptr = (iter % y) + 1;
232 			iter++;
233 		    }
234 		}
235 	    }
236 	}
237     for (i = 0; i < context->height; i++) {
238 	unsigned char *colors = context->colors;
239 	addr = context->vbuff + i * context->scanline;
240 	addr1 = addr + context->width;
241 	x = context->num_colors;
242 	while (addr < addr1) {
243 	    *addr = colors[*addr % x];
244 	    addr++;
245 	}
246     }
247 #ifdef STATISTICS
248     printf("iters %i guessed %i unguessed %i\n", iters, guessed, unguessed);
249     iters2 += iters;
250     guessed2 += guessed;
251     unguessed2 += unguessed;
252     printf("frames %i total %i iters %i guessed %i unguessed %i\n", frames2, total2, iters2, guessed2, unguessed2);
253 #endif
254 }
255