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