File: common\soi1.c
1 /*
2 * soi.c -- SOI
3 *
4 * Simultaneous Orbit Iteration Image Generation Method. Computes
5 * rectangular regions by tracking the orbits of only a few key points.
6 *
7 * Copyright (c) 1994-1997 Michael R. Ganss. All Rights Reserved.
8 *
9 * This file is distributed under the same conditions as
10 * AlmondBread. For further information see
11 * <URL:http://www.cs.tu-berlin.de/~rms/AlmondBread>.
12 *
13 */
14 #include <time.h>
15 #include <string.h>
16 #include <malloc.h>
17 #include "port.h"
18 #include "prototyp.h"
19
20 #define DBLS double
21 #define FABS(x) fabs(x)
22 #define FREXP(x,y) frexp(x,y)
23
24 #define TRUE 1
25 #define FALSE 0
26 #define EVERY 15
27 #define BASIN_COLOR 0
28
29 extern int rhombus_stack[10];
30 extern int rhombus_depth;
31 extern int max_rhombus_depth;
32 extern int minstackavail;
33 extern int minstack; /* need this much stack to recurse */
34 static DBLS twidth;
35 static DBLS equal;
36
37 #if 0
59 #endif
60
61 static long iteration(register DBLS cr, register DBLS ci,
62 register DBLS re, register DBLS im,
63 long start)
64 {
65 old.x = re;
66 old.y = im;
67 tempsqrx = sqr(old.x);
68 tempsqry = sqr(old.y);
69 floatparm = &init;
70 floatparm->x = cr;
71 floatparm->y = ci;
72 while(ORBITCALC()==0 && start < maxit)
73 start++;
74 if (start >= maxit)
75 start = BASIN_COLOR;
76 return(start);
77 }
78 #if 0
86 #endif
87
88 static void puthline(int x1,int y1,int x2,int color)
89 {
90 int x;
91 for(x=x1;x<=x2;x++)
92 (*plot)(x,y1,color);
93 }
94
95 static void putbox(int x1, int y1, int x2, int y2, int color)
96 {
97 for(; y1<=y2; y1++)
98 puthline(x1,y1,x2,color);
99 }
100
101 /* maximum side length beyond which we start regular scanning instead of
102 subdividing */
103 #define SCAN 16
104
105 /* pixel interleave used in scanning */
106 #define INTERLEAVE 4
107
108 /* compute the value of the interpolation polynomial at (x,y) */
109 #define GET_REAL(x,y) \
110 interpolate(cim1,midi,cim2,\
111 interpolate(cre1,midr,cre2,zre1,zre5,zre2,x),\
112 interpolate(cre1,midr,cre2,zre6,zre9,zre7,x),\
113 interpolate(cre1,midr,cre2,zre3,zre8,zre4,x),y)
114 #define GET_IMAG(x,y) \
115 interpolate(cre1,midr,cre2,\
116 interpolate(cim1,midi,cim2,zim1,zim6,zim3,y),\
117 interpolate(cim1,midi,cim2,zim5,zim9,zim8,y),\
118 interpolate(cim1,midi,cim2,zim2,zim7,zim4,y),x)
119
120 /* compute the value of the interpolation polynomial at (x,y)
121 from saved values before interpolation failed to stay within tolerance */
122 #define GET_SAVED_REAL(x,y) \
123 interpolate(cim1,midi,cim2,\
124 interpolate(cre1,midr,cre2,sr1,sr5,sr2,x),\
125 interpolate(cre1,midr,cre2,sr6,sr9,sr7,x),\
126 interpolate(cre1,midr,cre2,sr3,sr8,sr4,x),y)
127 #define GET_SAVED_IMAG(x,y) \
128 interpolate(cre1,midr,cre2,\
129 interpolate(cim1,midi,cim2,si1,si6,si3,y),\
130 interpolate(cim1,midi,cim2,si5,si9,si8,y),\
131 interpolate(cim1,midi,cim2,si2,si7,si4,y),x)
132
133 /* compute the value of the interpolation polynomial at (x,y)
134 during scanning. Here, key values do not change, so we can precompute
135 coefficients in one direction and simply evaluate the polynomial
136 during scanning. */
137 #define GET_SCAN_REAL(x,y) \
138 interpolate(cim1,midi,cim2,\
139 EVALUATE(cre1,midr,br10,br11,br12,x),\
140 EVALUATE(cre1,midr,br20,br21,br22,x),\
141 EVALUATE(cre1,midr,br30,br31,br32,x),y)
142 #define GET_SCAN_IMAG(x,y) \
143 interpolate(cre1,midr,cre2,\
144 EVALUATE(cim1,midi,bi10,bi11,bi12,y),\
145 EVALUATE(cim1,midi,bi20,bi21,bi22,y),\
146 EVALUATE(cim1,midi,bi30,bi31,bi32,y),x)
147
148 /* compute coefficients of Newton polynomial (b0,..,b2) from
149 (x0,w0),..,(x2,w2). */
150 #define INTERPOLATE(x0,x1,x2,w0,w1,w2,b0,b1,b2) \
151 b0=w0;\
152 b1=(w1-w0)/(x1-x0);\
153 b2=((w2-w1)/(x2-x1)-b1)/(x2-x0)
154
155 /* evaluate Newton polynomial given by (x0,b0),(x1,b1) at x:=t */
156 #define EVALUATE(x0,x1,b0,b1,b2,t) \
157 ((b2*(t-x1)+b1)*(t-x0)+b0)
158
159 /* Newton Interpolation.
160 It computes the value of the interpolation polynomial given by
161 (x0,w0)..(x2,w2) at x:=t */
162 static DBLS interpolate(DBLS x0, DBLS x1, DBLS x2,
163 DBLS w0, DBLS w1, DBLS w2,
164 DBLS t)
165 {
166 register DBLS b0=w0,b1=w1,b2=w2,b;
167
168 /*b0=(r0*b1-r1*b0)/(x1-x0);
169 b1=(r1*b2-r2*b1)/(x2-x1);
170 b0=(r0*b1-r2*b0)/(x2-x0);
171
172 return (DBLS)b0;*/
173 b=(b1-b0)/(x1-x0);
174 return (DBLS)((((b2-b1)/(x2-x1)-b)/(x2-x0))*(t-x1)+b)*(t-x0)+b0;
175 /*
176 if(t<x1)
177 return w0+((t-x0)/(x1-x0))*(w1-w0);
178 else
179 return w1+((t-x1)/(x2-x1))*(w2-w1);*/
180 }
181 #if (_MSC_VER >= 700)
182 #pragma code_seg ("soi3_text") /* place following in an overlay */
183 #endif
184
185 /* SOICompute - Perform simultaneous orbit iteration for a given rectangle
186
187 Input: cre1..cim2 : values defining the four corners of the rectangle
188 x1..y2 : corresponding pixel values
189 zre1..zim9 : intermediate iterated values of the key points (key values)
190
191 (cre1,cim1) (cre2,cim1)
192 (zre1,zim1) (zre5,zim5) (zre2,zim2)
193 +------------+------------+
194 | | |
195 | | |
196 (zre6,zim6) (zre9,zim9) (zre7,zim7)
197 | | |
198 | | |
199 +------------+------------+
200 (zre3,zim3) (zre8,zim8) (zre4,zim4)
201 (cre1,cim2) (cre2,cim2)
202
203 iter : current number of iterations
204 */
205 static DBLS zre1, zim1, zre2, zim2, zre3, zim3, zre4, zim4, zre5, zim5,
206 zre6, zim6, zre7, zim7, zre8, zim8, zre9, zim9;
207 /*
208 The purpose of this macro is to reduce the number of parameters of the
209 function rhombus(), since this is a recursive function, and stack space
210 under DOS is extremely limited.
211 */
212
213 #define RHOMBUS(CRE1,CRE2,CIM1,CIM2,X1,X2,Y1,Y2,ZRE1,ZIM1,ZRE2,ZIM2,ZRE3,ZIM3,\
214 ZRE4, ZIM4, ZRE5, ZIM5,ZRE6, ZIM6, ZRE7, ZIM7, ZRE8, ZIM8, ZRE9, ZIM9,ITER) \
215 zre1=(ZRE1);zim1=(ZIM1);\
216 zre2=(ZRE2);zim2=(ZIM2);\
217 zre3=(ZRE3);zim3=(ZIM3);\
218 zre4=(ZRE4);zim4=(ZIM4);\
219 zre5=(ZRE5);zim5=(ZIM5);\
220 zre6=(ZRE6);zim6=(ZIM6);\
221 zre7=(ZRE7);zim7=(ZIM7);\
222 zre8=(ZRE8);zim8=(ZIM8);\
223 zre9=(ZRE9);zim9=(ZIM9);\
224 status=rhombus((CRE1),(CRE2),(CIM1),(CIM2),(X1),(X2),(Y1),(Y2),(ITER))
225
226 static int rhombus(DBLS cre1, DBLS cre2, DBLS cim1, DBLS cim2,
227 int x1, int x2, int y1, int y2, long iter)
228 {
229 /* The following variables do not need their values saved */
230 /* used in scanning */
231 static long far savecolor, color, helpcolor;
232 static int far x,y,z,savex;
233
234 #if 0
235 static DBLS far re,im,restep,imstep,interstep,helpre;
236 static DBLS far zre,zim;
237 /* interpolation coefficients */
238 static DBLS far br10,br11,br12,br20,br21,br22,br30,br31,br32;
239 static DBLS far bi10,bi11,bi12,bi20,bi21,bi22,bi30,bi31,bi32;
240 /* ratio of interpolated test point to iterated one */
241 static DBLS far l1,l2;
242 /* squares of key values */
243 static DBLS far rq1,iq1;
244 static DBLS far rq2,iq2;
245 static DBLS far rq3,iq3;
246 static DBLS far rq4,iq4;
247 static DBLS far rq5,iq5;
248 static DBLS far rq6,iq6;
249 static DBLS far rq7,iq7;
250 static DBLS far rq8,iq8;
251 static DBLS far rq9,iq9;
252
253 /* test points */
254 static DBLS far cr1,cr2;
255 static DBLS far ci1,ci2;
256 static DBLS far tzr1,tzi1,tzr2,tzi2,tzr3,tzi3,tzr4,tzi4;
257 static DBLS far trq1,tiq1,trq2,tiq2,trq3,tiq3,trq4,tiq4; 258 #else
259 #define re mem_static[ 0]
260 #define im mem_static[ 1]
261 #define restep mem_static[ 2]
262 #define imstep mem_static[ 3]
263 #define interstep mem_static[ 4]
264 #define helpre mem_static[ 5]
265 #define zre mem_static[ 6]
266 #define zim mem_static[ 7]
267 #define br10 mem_static[ 8]
268 #define br11 mem_static[ 9]
269 #define br12 mem_static[10]
270 #define br20 mem_static[11]
271 #define br21 mem_static[12]
272 #define br22 mem_static[13]
273 #define br30 mem_static[14]
274 #define br31 mem_static[15]
275 #define br32 mem_static[16]
276 #define bi10 mem_static[17]
277 #define bi11 mem_static[18]
278 #define bi12 mem_static[19]
279 #define bi20 mem_static[20]
280 #define bi21 mem_static[21]
281 #define bi22 mem_static[22]
282 #define bi30 mem_static[23]
283 #define bi31 mem_static[24]
284 #define bi32 mem_static[25]
285 #define l1 mem_static[26]
286 #define l2 mem_static[27]
287 #define rq1 mem_static[28]
288 #define iq1 mem_static[29]
289 #define rq2 mem_static[30]
290 #define iq2 mem_static[31]
291 #define rq3 mem_static[32]
292 #define iq3 mem_static[33]
293 #define rq4 mem_static[34]
294 #define iq4 mem_static[35]
295 #define rq5 mem_static[36]
296 #define iq5 mem_static[37]
297 #define rq6 mem_static[38]
298 #define iq6 mem_static[39]
299 #define rq7 mem_static[40]
300 #define iq7 mem_static[41]
301 #define rq8 mem_static[42]
302 #define iq8 mem_static[43]
303 #define rq9 mem_static[44]
304 #define iq9 mem_static[45]
305 #define cr1 mem_static[46]
306 #define cr2 mem_static[47]
307 #define ci1 mem_static[48]
308 #define ci2 mem_static[49]
309 #define tzr1 mem_static[50]
310 #define tzi1 mem_static[51]
311 #define tzr2 mem_static[52]
312 #define tzi2 mem_static[53]
313 #define tzr3 mem_static[54]
314 #define tzi3 mem_static[55]
315 #define tzr4 mem_static[56]
316 #define tzi4 mem_static[57]
317 #define trq1 mem_static[58]
318 #define tiq1 mem_static[59]
319 #define trq2 mem_static[60]
320 #define tiq2 mem_static[61]
321 #define trq3 mem_static[62]
322 #define tiq3 mem_static[63]
323 #define trq4 mem_static[64]
324 #define tiq4 mem_static[65]
325
326 #endif
327 /* number of iterations before SOI iteration cycle */
328 static long far before;
329 static int far avail;
330
331 /* the variables below need to have local copies for recursive calls */
332 int far *mem_int;
333 DBLS far *mem;
334 DBLS far *mem_static;
335 /* center of rectangle */
336 DBLS midr=(cre1+cre2)/2,midi=(cim1+cim2)/2;
337
338 #if 0
339 /* saved values of key values */
340 DBLS sr1,si1,sr2,si2,sr3,si3,sr4,si4;
341 DBLS sr5,si5,sr6,si6,sr7,si7,sr8,si8,sr9,si9;
342 /* key values for subsequent rectangles */
343 DBLS re10,re11,re12,re13,re14,re15,re16,re17,re18,re19,re20,re21;
344 DBLS im10,im11,im12,im13,im14,im15,im16,im17,im18,im19,im20,im21;
345 DBLS re91,re92,re93,re94,im91,im92,im93,im94; 346 #else
347 #define esc1 mem_int[0]
348 #define esc2 mem_int[1]
349 #define esc3 mem_int[2]
350 #define esc4 mem_int[3]
351 #define esc5 mem_int[4]
352 #define esc6 mem_int[5]
353 #define esc7 mem_int[6]
354 #define esc8 mem_int[7]
355 #define esc9 mem_int[8]
356 #define tesc1 mem_int[9]
357 #define tesc2 mem_int[10]
358 #define tesc3 mem_int[11]
359 #define tesc4 mem_int[12]
360
361 #define sr1 mem[ 0]
362 #define si1 mem[ 1]
363 #define sr2 mem[ 2]
364 #define si2 mem[ 3]
365 #define sr3 mem[ 4]
366 #define si3 mem[ 5]
367 #define sr4 mem[ 6]
368 #define si4 mem[ 7]
369 #define sr5 mem[ 8]
370 #define si5 mem[ 9]
371 #define sr6 mem[10]
372 #define si6 mem[11]
373 #define sr7 mem[12]
374 #define si7 mem[13]
375 #define sr8 mem[14]
376 #define si8 mem[15]
377 #define sr9 mem[16]
378 #define si9 mem[17]
379 #define re10 mem[18]
380 #define re11 mem[19]
381 #define re12 mem[20]
382 #define re13 mem[21]
383 #define re14 mem[22]
384 #define re15 mem[23]
385 #define re16 mem[24]
386 #define re17 mem[25]
387 #define re18 mem[26]
388 #define re19 mem[27]
389 #define re20 mem[28]
390 #define re21 mem[29]
391 #define im10 mem[30]
392 #define im11 mem[31]
393 #define im12 mem[32]
394 #define im13 mem[33]
395 #define im14 mem[34]
396 #define im15 mem[35]
397 #define im16 mem[36]
398 #define im17 mem[37]
399 #define im18 mem[38]
400 #define im19 mem[39]
401 #define im20 mem[40]
402 #define im21 mem[41]
403 #define re91 mem[42]
404 #define re92 mem[43]
405 #define re93 mem[44]
406 #define re94 mem[45]
407 #define im91 mem[46]
408 #define im92 mem[47]
409 #define im93 mem[48]
410 #define im94 mem[49]
411 #endif
412
413 int status = 0;
414 rhombus_depth++;
415
416 #if 1
417 /* what we go through under DOS to deal with memory! We re-use
418 the sizeofstring array (8k). The first 660 bytes is for
419 static variables, then we make our own "stack" with copies
420 for each recursive call of rhombus() for the rest.
421 */
422 mem_int = (int far *)sizeofstring;
423 mem_static = (DBLS far *)(mem_int + 13);
424 mem = mem_static+ 66 + 50*rhombus_depth;
425 #endif
426
427 if((avail = stackavail()) < minstackavail)
428 minstackavail = avail;
429 if(rhombus_depth > max_rhombus_depth)
430 max_rhombus_depth = rhombus_depth;
431 rhombus_stack[rhombus_depth] = avail;
432
433 if(keypressed())
434 {
435 status = 1;
436 goto rhombus_done;
437 }
438 if(iter>maxit)
439 {
440 putbox(x1,y1,x2,y2,0);
441 status = 0;
442 goto rhombus_done;
443 }
444
445 if((y2-y1<=SCAN) || (avail < minstack))
446 {
447 /* finish up the image by scanning the rectangle */
448 scan:
449 INTERPOLATE(cre1,midr,cre2,zre1,zre5,zre2,br10,br11,br12);
450 INTERPOLATE(cre1,midr,cre2,zre6,zre9,zre7,br20,br21,br22);
451 INTERPOLATE(cre1,midr,cre2,zre3,zre8,zre4,br30,br31,br32);
452
453 INTERPOLATE(cim1,midi,cim2,zim1,zim6,zim3,bi10,bi11,bi12);
454 INTERPOLATE(cim1,midi,cim2,zim5,zim9,zim8,bi20,bi21,bi22);
455 INTERPOLATE(cim1,midi,cim2,zim2,zim7,zim4,bi30,bi31,bi32);
456
457 restep=(cre2-cre1)/(x2-x1);
458 imstep=(cim2-cim1)/(y2-y1);
459 interstep=INTERLEAVE*restep;
460
461 for(y=y1, im=cim1; y<y2; y++, im+=imstep)
462 {
463 if(keypressed())
464 {
465 status = 1;
466 goto rhombus_done;
467 }
468 zre=GET_SCAN_REAL(cre1,im);
469 zim=GET_SCAN_IMAG(cre1,im);
470 savecolor=iteration(cre1,im,zre,zim,iter);
471 if(savecolor < 0)
472 {
473 status = 1;
474 goto rhombus_done;
475 }
476 savex=x1;
477 for(x=x1+INTERLEAVE, re=cre1+interstep; x<x2;
478 x+=INTERLEAVE, re+=interstep)
479 {
480 zre=GET_SCAN_REAL(re,im);
481 zim=GET_SCAN_IMAG(re,im);
482
483 color=iteration(re,im,zre,zim,iter);
484 if(color < 0)
485 {
486 status = 1;
487 goto rhombus_done;
488 }
489 else if(color==savecolor)
490 continue;
491
492 for (z=x-1, helpre=re-restep; z>x-INTERLEAVE; z--,helpre-=restep)
493 {
494 zre=GET_SCAN_REAL(helpre,im);
495 zim=GET_SCAN_IMAG(helpre,im);
496 helpcolor=iteration(helpre,im,zre,zim,iter);
497 if(helpcolor < 0)
498 {
499 status = 1;
500 goto rhombus_done;
501 }
502 else if(helpcolor==savecolor)
503 break;
504 (*plot)(z,y,(int)(helpcolor&255));
505 }
506
507 if(savex<z)
508 puthline(savex, y, z, (int)(savecolor&255));
509 else
510 (*plot)(savex, y, (int)(savecolor&255));
511
512 savex = x;
513 savecolor = color;
514 }
515
516 for (z=x2-1, helpre=cre2-restep; z>savex; z--,helpre-=restep)
517 {
518 zre=GET_SCAN_REAL(helpre,im);
519 zim=GET_SCAN_IMAG(helpre,im);
520 helpcolor=iteration(helpre,im,zre,zim,iter);
521 if(helpcolor < 0)
522 {
523 status = 1;
524 goto rhombus_done;
525 }
526 else if(helpcolor==savecolor)
527 break;
528
529 (*plot)(z,y,(int)(helpcolor&255));
530 }
531
532 if(savex<z)
533 puthline(savex, y, z, (int)(savecolor&255));
534 else
535 (*plot)(savex, y, (int)(savecolor&255));
536 }
537 status = 0;
538 goto rhombus_done;
539 }
540
541 rq1=zre1*zre1; iq1=zim1*zim1;
542 rq2=zre2*zre2; iq2=zim2*zim2;
543 rq3=zre3*zre3; iq3=zim3*zim3;
544 rq4=zre4*zre4; iq4=zim4*zim4;
545 rq5=zre5*zre5; iq5=zim5*zim5;
546 rq6=zre6*zre6; iq6=zim6*zim6;
547 rq7=zre7*zre7; iq7=zim7*zim7;
548 rq8=zre8*zre8; iq8=zim8*zim8;
549 rq9=zre9*zre9; iq9=zim9*zim9;
550
551 cr1=0.75*cre1+0.25*cre2; cr2=0.25*cre1+0.75*cre2;
552 ci1=0.75*cim1+0.25*cim2; ci2=0.25*cim1+0.75*cim2;
553
554 tzr1=GET_REAL(cr1,ci1);
555 tzi1=GET_IMAG(cr1,ci1);
556
557 tzr2=GET_REAL(cr2,ci1);
558 tzi2=GET_IMAG(cr2,ci1);
559
560 tzr3=GET_REAL(cr1,ci2);
561 tzi3=GET_IMAG(cr1,ci2);
562
563 tzr4=GET_REAL(cr2,ci2);
564 tzi4=GET_IMAG(cr2,ci2);
565
566 trq1=tzr1*tzr1;
567 tiq1=tzi1*tzi1;
568
569 trq2=tzr2*tzr2;
570 tiq2=tzi2*tzi2;
571
572 trq3=tzr3*tzr3;
573 tiq3=tzi3*tzi3;
574
575 trq4=tzr4*tzr4;
576 tiq4=tzi4*tzi4;
577
578 before=iter;
579
580 for(;;)
581 {
582 sr1=zre1; si1=zim1;
583 sr2=zre2; si2=zim2;
584 sr3=zre3; si3=zim3;
585 sr4=zre4; si4=zim4;
586 sr5=zre5; si5=zim5;
587 sr6=zre6; si6=zim6;
588 sr7=zre7; si7=zim7;
589 sr8=zre8; si8=zim8;
590 sr9=zre9; si9=zim9;
591
592
593 #define SOI_ORBIT1(zr,rq,zi,iq,cr,ci,esc) \
594 tempsqrx = rq;\
595 tempsqry = iq;\
596 old.x = zr;\
597 old.y = zi;\
598 floatparm->x = cr;\
599 floatparm->y = ci;\
600 esc = ORBITCALC();\
601 rq = tempsqrx;\
602 iq = tempsqry;\
603 zr = new.x;\
604 zi = new.y
605
606 #define SOI_ORBIT(zr,rq,zi,iq,cr,ci,esc) \
607 zi=(zi+zi)*zr+ci;\
608 zr=rq-iq+cr;\
609 rq=zr*zr;\
610 iq=zi*zi;\
611 esc = ((rq+iq)>16.0)?1:0
612
613 /* iterate key values */
614 SOI_ORBIT(zre1,rq1,zim1,iq1,cre1,cim1,esc1);
615 /*
616 zim1=(zim1+zim1)*zre1+cim1;
617 zre1=rq1-iq1+cre1;
618 rq1=zre1*zre1;
619 iq1=zim1*zim1;
620 */
621 SOI_ORBIT(zre2,rq2,zim2,iq2,cre2,cim1,esc2);
622 /*
623 zim2=(zim2+zim2)*zre2+cim1;
624 zre2=rq2-iq2+cre2;
625 rq2=zre2*zre2;
626 iq2=zim2*zim2;
627 */
628 SOI_ORBIT(zre3,rq3,zim3,iq3,cre1,cim2,esc3);
629 /*
630 zim3=(zim3+zim3)*zre3+cim2;
631 zre3=rq3-iq3+cre1;
632 rq3=zre3*zre3;
633 iq3=zim3*zim3;
634 */
635 SOI_ORBIT(zre4,rq4,zim4,iq4,cre2,cim2,esc4);
636 /*
637 zim4=(zim4+zim4)*zre4+cim2;
638 zre4=rq4-iq4+cre2;
639 rq4=zre4*zre4;
640 iq4=zim4*zim4;
641 */
642 SOI_ORBIT(zre5,rq5,zim5,iq5,midr,cim1,esc5);
643 /*
644 zim5=(zim5+zim5)*zre5+cim1;
645 zre5=rq5-iq5+midr;
646 rq5=zre5*zre5;
647 iq5=zim5*zim5;
648 */
649 SOI_ORBIT(zre6,rq6,zim6,iq6,cre1,midi,esc6);
650 /*
651 zim6=(zim6+zim6)*zre6+midi;
652 zre6=rq6-iq6+cre1;
653 rq6=zre6*zre6;
654 iq6=zim6*zim6;
655 */
656 SOI_ORBIT(zre7,rq7,zim7,iq7,cre2,midi,esc7);
657 /*
658 zim7=(zim7+zim7)*zre7+midi;
659 zre7=rq7-iq7+cre2;
660 rq7=zre7*zre7;
661 iq7=zim7*zim7;
662 */
663 SOI_ORBIT(zre8,rq8,zim8,iq8,midr,cim2,esc8);
664 /*
665 zim8=(zim8+zim8)*zre8+cim2;
666 zre8=rq8-iq8+midr;
667 rq8=zre8*zre8;
668 iq8=zim8*zim8;
669 */
670 SOI_ORBIT(zre9,rq9,zim9,iq9,midr,midi,esc9);
671 /*
672 zim9=(zim9+zim9)*zre9+midi;
673 zre9=rq9-iq9+midr;
674 rq9=zre9*zre9;
675 iq9=zim9*zim9;
676 */
677 /* iterate test point */
678 SOI_ORBIT(tzr1,trq1,tzi1,tiq1,cr1,ci1,tesc1);
679 /*
680 tzi1=(tzi1+tzi1)*tzr1+ci1;
681 tzr1=trq1-tiq1+cr1;
682 trq1=tzr1*tzr1;
683 tiq1=tzi1*tzi1;
684 */
685
686 SOI_ORBIT(tzr2,trq2,tzi2,tiq2,cr2,ci1,tesc2);
687 /*
688 tzi2=(tzi2+tzi2)*tzr2+ci1;
689 tzr2=trq2-tiq2+cr2;
690 trq2=tzr2*tzr2;
691 tiq2=tzi2*tzi2;
692 */
693 SOI_ORBIT(tzr3,trq3,tzi3,tiq3,cr1,ci2,tesc3);
694 /*
695 tzi3=(tzi3+tzi3)*tzr3+ci2;
696 tzr3=trq3-tiq3+cr1;
697 trq3=tzr3*tzr3;
698 tiq3=tzi3*tzi3;
699 */
700 SOI_ORBIT(tzr4,trq4,tzi4,tiq4,cr2,ci2,tesc4);
701 /*
702 tzi4=(tzi4+tzi4)*tzr4+ci2;
703 tzr4=trq4-tiq4+cr2;
704 trq4=tzr4*tzr4;
705 tiq4=tzi4*tzi4;
706 */
707 iter++;
708
709 /* if one of the iterated values bails out, subdivide */
710 /*
711 if((rq1+iq1)>16.0||
712 (rq2+iq2)>16.0||
713 (rq3+iq3)>16.0||
714 (rq4+iq4)>16.0||
715 (rq5+iq5)>16.0||
716 (rq6+iq6)>16.0||
717 (rq7+iq7)>16.0||
718 (rq8+iq8)>16.0||
719 (rq9+iq9)>16.0||
720 (trq1+tiq1)>16.0||
721 (trq2+tiq2)>16.0||
722 (trq3+tiq3)>16.0||
723 (trq4+tiq4)>16.0)
724 break;
725 */
726 if(esc1||esc2||esc3||esc4||esc5||esc6||esc7||esc8||esc9||
727 tesc1||tesc2||tesc3||tesc4)
728 break;
729
730 /* if maximum number of iterations is reached, the whole rectangle
731 can be assumed part of M. This is of course best case behavior
732 of SOI, we seldomly get there */
733 if(iter>maxit)
734 {
735 putbox(x1,y1,x2,y2,0);
736 status = 0;
737 goto rhombus_done;
738 }
739
740 /* now for all test points, check whether they exceed the
741 allowed tolerance. if so, subdivide */
742 l1=GET_REAL(cr1,ci1);
743 l1=(tzr1==0.0)?
744 (l1==0.0)?1.0:1000.0:
745 l1/tzr1;
746 if(FABS(1.0-l1)>twidth)
747 break;
748
749 l2=GET_IMAG(cr1,ci1);
750 l2=(tzi1==0.0)?
751 (l2==0.0)?1.0:1000.0:
752 l2/tzi1;
753 if(FABS(1.0-l2)>twidth)
754 break;
755
756 l1=GET_REAL(cr2,ci1);
757 l1=(tzr2==0.0)?
758 (l1==0.0)?1.0:1000.0:
759 l1/tzr2;
760 if(FABS(1.0-l1)>twidth)
761 break;
762
763 l2=GET_IMAG(cr2,ci1);
764 l2=(tzi2==0.0)?
765 (l2==0.0)?1.0:1000.0:
766 l2/tzi2;
767 if(FABS(1.0-l2)>twidth)
768 break;
769
770 l1=GET_REAL(cr1,ci2);
771 l1=(tzr3==0.0)?
772 (l1==0.0)?1.0:1000.0:
773 l1/tzr3;
774 if(FABS(1.0-l1)>twidth)
775 break;
776
777 l2=GET_IMAG(cr1,ci2);
778 l2=(tzi3==0.0)?
779 (l2==0.0)?1.0:1000.0:
780 l2/tzi3;
781 if(FABS(1.0-l2)>twidth)
782 break;
783
784 l1=GET_REAL(cr2,ci2);
785 l1=(tzr4==0.0)?
786 (l1==0.0)?1.0:1000.0:
787 l1/tzr4;
788 if(FABS(1.0-l1)>twidth)
789 break;
790
791 l2=GET_IMAG(cr2,ci2);
792 l2=(tzi4==0.0)?
793 (l2==0.0)?1.0:1000.0:
794 l2/tzi4;
795 if(FABS(1.0-l2)>twidth)
796 break;
797 }
798
799 iter--;
800
801 /* this is a little heuristic I tried to improve performance. */
802 if(iter-before<10)
803 {
804 zre1=sr1; zim1=si1;
805 zre2=sr2; zim2=si2;
806 zre3=sr3; zim3=si3;
807 zre4=sr4; zim4=si4;
808 zre5=sr5; zim5=si5;
809 zre6=sr6; zim6=si6;
810 zre7=sr7; zim7=si7;
811 zre8=sr8; zim8=si8;
812 zre9=sr9; zim9=si9;
813 goto scan;
814 }
815
816 /* compute key values for subsequent rectangles */
817
818 re10=interpolate(cre1,midr,cre2,sr1,sr5,sr2,cr1);
819 im10=interpolate(cre1,midr,cre2,si1,si5,si2,cr1);
820
821 re11=interpolate(cre1,midr,cre2,sr1,sr5,sr2,cr2);
822 im11=interpolate(cre1,midr,cre2,si1,si5,si2,cr2);
823
824 re20=interpolate(cre1,midr,cre2,sr3,sr8,sr4,cr1);
825 im20=interpolate(cre1,midr,cre2,si3,si8,si4,cr1);
826
827 re21=interpolate(cre1,midr,cre2,sr3,sr8,sr4,cr2);
828 im21=interpolate(cre1,midr,cre2,si3,si8,si4,cr2);
829
830 re15=interpolate(cre1,midr,cre2,sr6,sr9,sr7,cr1);
831 im15=interpolate(cre1,midr,cre2,si6,si9,si7,cr1);
832
833 re16=interpolate(cre1,midr,cre2,sr6,sr9,sr7,cr2);
834 im16=interpolate(cre1,midr,cre2,si6,si9,si7,cr2);
835
836 re12=interpolate(cim1,midi,cim2,sr1,sr6,sr3,ci1);
837 im12=interpolate(cim1,midi,cim2,si1,si6,si3,ci1);
838
839 re14=interpolate(cim1,midi,cim2,sr2,sr7,sr4,ci1);
840 im14=interpolate(cim1,midi,cim2,si2,si7,si4,ci1);
841
842 re17=interpolate(cim1,midi,cim2,sr1,sr6,sr3,ci2);
843 im17=interpolate(cim1,midi,cim2,si1,si6,si3,ci2);
844
845 re19=interpolate(cim1,midi,cim2,sr2,sr7,sr4,ci2);
846 im19=interpolate(cim1,midi,cim2,si2,si7,si4,ci2);
847
848 re13=interpolate(cim1,midi,cim2,sr5,sr9,sr8,ci1);
849 im13=interpolate(cim1,midi,cim2,si5,si9,si8,ci1);
850
851 re18=interpolate(cim1,midi,cim2,sr5,sr9,sr8,ci2);
852 im18=interpolate(cim1,midi,cim2,si5,si9,si8,ci2);
853
854 re91=GET_SAVED_REAL(cr1,ci1);
855 re92=GET_SAVED_REAL(cr2,ci1);
856 re93=GET_SAVED_REAL(cr1,ci2);
857 re94=GET_SAVED_REAL(cr2,ci2);
858
859 im91=GET_SAVED_IMAG(cr1,ci1);
860 im92=GET_SAVED_IMAG(cr2,ci1);
861 im93=GET_SAVED_IMAG(cr1,ci2);
862 im94=GET_SAVED_IMAG(cr2,ci2);
863
864 RHOMBUS(cre1,midr,cim1,midi,x1,((x1+x2)>>1),y1,((y1+y2)>>1),
865 sr1,si1,
866 sr5,si5,
867 sr6,si6,
868 sr9,si9,
869 re10,im10,
870 re12,im12,
871 re13,im13,
872 re15,im15,
873 re91,im91,
874 iter);
875 RHOMBUS(midr,cre2,cim1,midi,(x1+x2)>>1,x2,y1,(y1+y2)>>1,
876 sr5,si5,
877 sr2,si2,
878 sr9,si9,
879 sr7,si7,
880 re11,im11,
881 re13,im13,
882 re14,im14,
883 re16,im16,
884 re92,im92,
885 iter);
886 RHOMBUS(cre1,midr,midi,cim2,x1,(x1+x2)>>1,(y1+y2)>>1,y2,
887 sr6,si6,
888 sr9,si9,
889 sr3,si3,
890 sr8,si8,
891 re15,im15,
892 re17,im17,
893 re18,im18,
894 re20,im20,
895 re93,im93,
896 iter);
897 RHOMBUS(midr,cre2,midi,cim2,(x1+x2)>>1,x2,(y1+y2)>>1,y2,
898 sr9,si9,
899 sr7,si7,
900 sr8,si8,
901 sr4,si4,
902 re16,im16,
903 re18,im18,
904 re19,im19,
905 re21,im21,
906 re94,im94,
907 iter);
908 rhombus_done:
909 rhombus_depth--;
910 return(status);
911 }
912
913 #if (_MSC_VER >= 700)
914 #pragma code_seg ()
915 #endif
916
917 void soi(void)
918 {
919 int status;
920 DBLS tolerance=0.1;
921 DBLS stepx, stepy;
922 DBLS xxminl, xxmaxl, yyminl, yymaxl;
923 minstackavail = 30000;
924 rhombus_depth = -1;
925 max_rhombus_depth = 0;
926 if(bf_math)
927 {
928 xxminl = (DBLS)bftofloat(bfxmin);
929 yyminl = (DBLS)bftofloat(bfymin);
930 xxmaxl = (DBLS)bftofloat(bfxmax);
931 yymaxl = (DBLS)bftofloat(bfymax);
932 }
933 else
934 {
935 xxminl = xxmin;
936 yyminl = yymin;
937 xxmaxl = xxmax;
938 yymaxl = yymax;
939 }
940 twidth=tolerance/(xdots-1);
941 stepx = (xxmaxl - xxminl) / xdots;
942 stepy = (yyminl - yymaxl) / ydots;
943 equal = (stepx < stepy ? stepx : stepy);
944
945 RHOMBUS(xxminl,xxmaxl,yymaxl,yyminl,
946 0,xdots,0,ydots,
947 xxminl,yymaxl,
948 xxmaxl,yymaxl,
949 xxminl,yyminl,
950 xxmaxl,yyminl,
951 (xxmaxl+xxminl)/2,yymaxl,
952 xxminl,(yymaxl+yyminl)/2,
953 xxmaxl,(yymaxl+yyminl)/2,
954 (xxmaxl+xxminl)/2,yyminl,
955 (xxminl+xxmaxl)/2,(yymaxl+yyminl)/2,
956 1);
957 }
958