File: common\fracsubr.c
1 /*
2 FRACSUBR.C contains subroutines which belong primarily to CALCFRAC.C and
3 FRACTALS.C, i.e. which are non-fractal-specific fractal engine subroutines.
4 */
5
6 #ifndef USE_VARARGS
7 #include <stdarg.h>
8 #else
9 #include <varargs.h> 10 #endif
11
12 #ifndef XFRACT
13 #include <sys/timeb.h>
14 #endif
15 #include <sys/types.h>
16 #include <time.h>
17 /* see Fractint.c for a description of the "include" hierarchy */
18 #include "port.h"
19 #include "prototyp.h"
20 #include "fractype.h"
21
22
23 /* routines in this module */
24
25 static long _fastcall fudgetolong(double d);
26 static double _fastcall fudgetodouble(long l);
27 static void _fastcall adjust_to_limits(double);
28 static void _fastcall smallest_add(double *);
29 static int _fastcall ratio_bad(double,double);
30 static void _fastcall plotdorbit(double,double,int);
31 static int _fastcall combine_worklist(void);
32
33 static void _fastcall adjust_to_limitsbf(double);
34 static void _fastcall smallest_add_bf(bf_t);
35 int resume_len; /* length of resume info */
36 static int resume_offset; /* offset in resume info gets */
37 int taborhelp; /* kludge for sound and tab or help key press */
38
39 #define FUDGEFACTOR 29 /* fudge all values up by 2**this */
40 #define FUDGEFACTOR2 24 /* (or maybe this) */
41
42 void set_grid_pointers()
43 {
44 dx0 = MK_FP(extraseg,0);
45 dy1 = (dx1 = (dy0 = dx0 + xdots) + ydots) + ydots;
46 lx0 = (long far *) dx0;
47 ly1 = (lx1 = (ly0 = lx0 + xdots) + ydots) + ydots;
48 set_pixel_calc_functions();
49 }
50
51 void fill_dx_array(void)
52 {
53 int i;
54 if(use_grid)
55 {
56 dx0[0] = xxmin; /* fill up the x, y grids */
57 dy0[0] = yymax;
58 dx1[0] = dy1[0] = 0;
59 for (i = 1; i < xdots; i++ ) {
60 dx0[i] = (double)(dx0[0] + i*delxx);
61 dy1[i] = (double)(dy1[0] - i*delyy2);
62 }
63 for (i = 1; i < ydots; i++ ) {
64 dy0[i] = (double)(dy0[0] - i*delyy);
65 dx1[i] = (double)(dx1[0] + i*delxx2);
66 }
67 }
68 }
69 void fill_lx_array(void)
70 {
71 int i;
72 /* note that lx1 & ly1 values can overflow into sign bit; since */
73 /* they're used only to add to lx0/ly0, 2s comp straightens it out */
74 if(use_grid)
75 {
76 lx0[0] = xmin; /* fill up the x, y grids */
77 ly0[0] = ymax;
78 lx1[0] = ly1[0] = 0;
79 for (i = 1; i < xdots; i++ ) {
80 lx0[i] = lx0[i-1] + delx;
81 ly1[i] = ly1[i-1] - dely2;
82 }
83 for (i = 1; i < ydots; i++ ) {
84 ly0[i] = ly0[i-1] - dely;
85 lx1[i] = lx1[i-1] + delx2;
86 }
87 }
88 }
89
90 void fractal_floattobf(void)
91 {
92 int i;
93 init_bf_dec(getprecdbl(CURRENTREZ));
94 floattobf(bfxmin,xxmin);
95 floattobf(bfxmax,xxmax);
96 floattobf(bfymin,yymin);
97 floattobf(bfymax,yymax);
98 floattobf(bfx3rd,xx3rd);
99 floattobf(bfy3rd,yy3rd);
100
101 for (i = 0; i < MAXPARAMS; i++)
102 if(typehasparm(fractype,i,NULL))
103 floattobf(bfparms[i],param[i]);
104 calc_status = 0;
105 }
106
107
108 #ifdef _MSC_VER
109 #if _MSC_VER == 800
110 /* MSC8 doesn't correctly calculate the address of certain arrays here */
111 #pragma optimize( "", off )
112 #endif
113 #endif
114
115 int use_grid;
116
117 void calcfracinit(void) /* initialize a *pile* of stuff for fractal calculation */
118 {
119 int tries = 0;
120 int i, gotprec;
121 long xytemp;
122 double ftemp;
123 coloriter=oldcoloriter = 0L;
124 for(i=0;i<10;i++)
125 rhombus_stack[i] = 0;
126
127 /* set up grid array compactly leaving space at end */
128 /* space req for grid is 2(xdots+ydots)*sizeof(long or double) */
129 /* space available in extraseg is 65536 Bytes */
130 xytemp = xdots + ydots;
131 if( ((usr_floatflag == 0) && (xytemp * sizeof(long) > 32768)) ||
132 ((usr_floatflag == 1) && (xytemp * sizeof(double) > 32768)) ||
133 debugflag == 3800)
134 {
135 use_grid=0;
136 floatflag = usr_floatflag = 1;
137 }
138 else
139 use_grid=1;
140
141 set_grid_pointers();
142
143 if(!(curfractalspecific->flags & BF_MATH))
144 {
145 int tofloat;
146 if((tofloat=curfractalspecific->tofloat) == NOFRACTAL)
147 bf_math = 0;
148 else if(!(fractalspecific[tofloat].flags & BF_MATH))
149 bf_math = 0;
150 else if(bf_math)
151 {
152 curfractalspecific = &fractalspecific[tofloat];
153 fractype = tofloat;
154 }
155 }
156
157 /* switch back to double when zooming out if using arbitrary precision */
158 if(bf_math)
159 {
160 gotprec=getprecbf(CURRENTREZ);
161 if((gotprec <= DBL_DIG+1 && debugflag != 3200) || math_tol[1] >= 1.0)
162 {
163 bfcornerstofloat();
164 bf_math = 0;
165 }
166 else
167 init_bf_dec(gotprec);
168 }
169 else if((fractype==MANDEL || fractype==MANDELFP) && debugflag==3200)
170 {
171 fractype=MANDELFP;
172 curfractalspecific = &fractalspecific[MANDELFP];
173 fractal_floattobf();
174 usr_floatflag = 1;
175 }
176 else if((fractype==JULIA || fractype==JULIAFP) && debugflag==3200)
177 {
178 fractype=JULIAFP;
179 curfractalspecific = &fractalspecific[JULIAFP];
180 fractal_floattobf();
181 usr_floatflag = 1;
182 }
183 else if((fractype==LMANDELZPOWER || fractype==FPMANDELZPOWER) && debugflag==3200)
184 {
185 fractype=FPMANDELZPOWER;
186 curfractalspecific = &fractalspecific[FPMANDELZPOWER];
187 fractal_floattobf();
188 usr_floatflag = 1;
189 }
190 else if((fractype==LJULIAZPOWER || fractype==FPJULIAZPOWER) && debugflag==3200)
191 {
192 fractype=FPJULIAZPOWER;
193 curfractalspecific = &fractalspecific[FPJULIAZPOWER];
194 fractal_floattobf();
195 usr_floatflag = 1;
196 }
197 else
198 free_bf_vars();
199 if(bf_math)
200 floatflag=1;
201 else
202 floatflag = usr_floatflag;
203 if (calc_status == 2) { /* on resume, ensure floatflag correct */
204 if (curfractalspecific->isinteger)
205 floatflag = 0;
206 else
207 floatflag = 1;
208 }
209 /* if floating pt only, set floatflag for TAB screen */
210 if (!curfractalspecific->isinteger && curfractalspecific->tofloat == NOFRACTAL)
211 floatflag = 1;
212 if (usr_stdcalcmode == 's') {
213 if (fractype == MANDEL || fractype == MANDELFP)
214 floatflag = 1;
215 else
216 usr_stdcalcmode = '1';
217 }
218
219 init_restart:
220
221 /* the following variables may be forced to a different setting due to
222 calc routine constraints; usr_xxx is what the user last said is wanted,
223 xxx is what we actually do in the current situation */
224 stdcalcmode = usr_stdcalcmode;
225 periodicitycheck = usr_periodicitycheck;
226 distest = usr_distest;
227 biomorph = usr_biomorph;
228
229 potflag = 0;
230 if (potparam[0] != 0.0
231 && colors >= 64
232 && (curfractalspecific->calctype == StandardFractal
233 || curfractalspecific->calctype == calcmand
234 || curfractalspecific->calctype == calcmandfp)) {
235 potflag = 1;
236 distest = usr_distest = 0; /* can't do distest too */
237 }
238
239 if (distest)
240 floatflag = 1; /* force floating point for dist est */
241
242 if (floatflag) { /* ensure type matches floatflag */
243 if (curfractalspecific->isinteger != 0
244 && curfractalspecific->tofloat != NOFRACTAL)
245 fractype = curfractalspecific->tofloat;
246 }
247 else {
248 if (curfractalspecific->isinteger == 0
249 && curfractalspecific->tofloat != NOFRACTAL)
250 fractype = curfractalspecific->tofloat;
251 }
252 /* match Julibrot with integer mode of orbit */
253 if(fractype == JULIBROTFP && fractalspecific[neworbittype].isinteger)
254 {
255 int i;
256 if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
257 neworbittype = i;
258 else
259 fractype = JULIBROT;
260 }
261 else if(fractype == JULIBROT && fractalspecific[neworbittype].isinteger==0)
262 {
263 int i;
264 if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
265 neworbittype = i;
266 else
267 fractype = JULIBROTFP;
268 }
269
270 curfractalspecific = &fractalspecific[fractype];
271
272 integerfractal = curfractalspecific->isinteger;
273
274 /* if (fractype == JULIBROT)
275 rqlim = 4;
276 else */ if (potflag && potparam[2] != 0.0)
277 rqlim = potparam[2];
278 /* else if (decomp[0] > 0 && decomp[1] > 0)
279 rqlim = (double)decomp[1]; */
280 else if (bailout) /* user input bailout */
281 rqlim = bailout;
282 else if (biomorph != -1) /* biomorph benefits from larger bailout */
283 rqlim = 100;
284 else
285 rqlim = curfractalspecific->orbit_bailout;
286 if (integerfractal) /* the bailout limit mustn't be too high here */
287 if (rqlim > 127.0) rqlim = 127.0;
288
289 if ((curfractalspecific->flags&NOROTATE) != 0) {
290 /* ensure min<max and unrotated rectangle */
291 if (xxmin > xxmax) { ftemp = xxmax; xxmax = xxmin; xxmin = ftemp; }
292 if (yymin > yymax) { ftemp = yymax; yymax = yymin; yymin = ftemp; }
293 xx3rd = xxmin; yy3rd = yymin;
294 }
295
296 /* set up bitshift for integer math */
297 bitshift = FUDGEFACTOR2; /* by default, the smaller shift */
298 if (integerfractal > 1) /* use specific override from table */
299 bitshift = integerfractal;
300 if (integerfractal == 0) { /* float? */
301 if ((i = curfractalspecific->tofloat) != NOFRACTAL) /* -> int? */
302 {
303 if (fractalspecific[i].isinteger > 1) /* specific shift? */
304 bitshift = fractalspecific[i].isinteger;
305 }
306 else
307 bitshift = 16; /* to allow larger corners */
308 }
309 /* We want this code if we're using the assembler calcmand */
310 if (fractype == MANDEL || fractype == JULIA) { /* adust shift bits if.. */
311 if (potflag == 0 /* not using potential */
312 && (param[0] > -2.0 && param[0] < 2.0) /* parameters not too large */
313 && (param[1] > -2.0 && param[1] < 2.0)
314 && !invert /* and not inverting */
315 && biomorph == -1 /* and not biomorphing */
316 && rqlim <= 4.0 /* and bailout not too high */
317 && (outside > -2 || outside < -6) /* and no funny outside stuff */
318 && debugflag != 1234 /* and not debugging */
319 && closeprox <= 2.0 /* and closeprox not too large */
320 && bailoutest == Mod) /* and bailout test = mod */
321 bitshift = FUDGEFACTOR; /* use the larger bitshift */
322 }
323
324 fudge = 1L << bitshift;
325
326 l_at_rad = fudge/32768L;
327 f_at_rad = 1.0/32768L;
328
329 /* now setup arrays of real coordinates corresponding to each pixel */
330 if(bf_math)
331 adjust_to_limitsbf(1.0); /* make sure all corners in valid range */
332 else
333 {
334 adjust_to_limits(1.0); /* make sure all corners in valid range */
335 delxx = (LDBL)(xxmax - xx3rd) / (LDBL)dxsize; /* calculate stepsizes */
336 delyy = (LDBL)(yymax - yy3rd) / (LDBL)dysize;
337 delxx2 = (LDBL)(xx3rd - xxmin) / (LDBL)dysize;
338 delyy2 = (LDBL)(yy3rd - yymin) / (LDBL)dxsize;
339 fill_dx_array();
340 }
341
342 if(fractype != CELLULAR && fractype != ANT) /* fudgetolong fails w >10 digits in double */
343 {
344 creal = fudgetolong(param[0]); /* integer equivs for it all */
345 cimag = fudgetolong(param[1]);
346 xmin = fudgetolong(xxmin);
347 xmax = fudgetolong(xxmax);
348 x3rd = fudgetolong(xx3rd);
349 ymin = fudgetolong(yymin);
350 ymax = fudgetolong(yymax);
351 y3rd = fudgetolong(yy3rd);
352 delx = fudgetolong((double)delxx);
353 dely = fudgetolong((double)delyy);
354 delx2 = fudgetolong((double)delxx2);
355 dely2 = fudgetolong((double)delyy2);
356 }
357
358 /* skip this if plasma to avoid 3d problems */
359 /* skip if bf_math to avoid extraseg conflict with dx0 arrays */
360 /* skip if ifs, ifs3d, or lsystem to avoid crash when mathtolerance */
361 /* is set. These types don't auto switch between float and integer math */
362 if (fractype != PLASMA && bf_math == 0
363 && fractype != IFS && fractype != IFS3D && fractype != LSYSTEM)
364 {
365 if (integerfractal && !invert && use_grid)
366 {
367 if ( (delx == 0 && delxx != 0.0)
368 || (delx2 == 0 && delxx2 != 0.0)
369 || (dely == 0 && delyy != 0.0)
370 || (dely2 == 0 && delyy2 != 0.0) )
371 goto expand_retry;
372
373 fill_lx_array(); /* fill up the x,y grids */
374 /* past max res? check corners within 10% of expected */
375 if ( ratio_bad((double)lx0[xdots-1]-xmin,(double)xmax-x3rd)
376 || ratio_bad((double)ly0[ydots-1]-ymax,(double)y3rd-ymax)
377 || ratio_bad((double)lx1[(ydots>>1)-1],((double)x3rd-xmin)/2)
378 || ratio_bad((double)ly1[(xdots>>1)-1],((double)ymin-y3rd)/2) )
379 {
380 expand_retry:
381 if (integerfractal /* integer fractal type? */
382 && curfractalspecific->tofloat != NOFRACTAL)
383 floatflag = 1; /* switch to floating pt */
384 else
385 adjust_to_limits(2.0); /* double the size */
386 if (calc_status == 2) /* due to restore of an old file? */
387 calc_status = 0; /* whatever, it isn't resumable */
388 goto init_restart;
389 } /* end if ratio bad */
390
391 /* re-set corners to match reality */
392 xmax = lx0[xdots-1] + lx1[ydots-1];
393 ymin = ly0[ydots-1] + ly1[xdots-1];
394 x3rd = xmin + lx1[ydots-1];
395 y3rd = ly0[ydots-1];
396 xxmin = fudgetodouble(xmin);
397 xxmax = fudgetodouble(xmax);
398 xx3rd = fudgetodouble(x3rd);
399 yymin = fudgetodouble(ymin);
400 yymax = fudgetodouble(ymax);
401 yy3rd = fudgetodouble(y3rd);
402 } /* end if (integerfractal && !invert && use_grid) */
403 else
404 {
405 double dx0,dy0,dx1,dy1;
406 /* set up dx0 and dy0 analogs of lx0 and ly0 */
407 /* put fractal parameters in doubles */
408 dx0 = xxmin; /* fill up the x, y grids */
409 dy0 = yymax;
410 dx1 = dy1 = 0;
411 /* this way of defining the dx and dy arrays is not the most
412 accurate, but it is kept because it is used to determine
413 the limit of resolution */
414 for (i = 1; i < xdots; i++ ) {
415 dx0 = (double)(dx0 + (double)delxx);
416 dy1 = (double)(dy1 - (double)delyy2);
417 }
418 for (i = 1; i < ydots; i++ ) {
419 dy0 = (double)(dy0 - (double)delyy);
420 dx1 = (double)(dx1 + (double)delxx2);
421 }
422 if(bf_math == 0) /* redundant test, leave for now */
423 {
424 double testx_try, testx_exact;
425 double testy_try, testy_exact;
426 /* Following is the old logic for detecting failure of double
427 precision. It has two advantages: it is independent of the
428 representation of numbers, and it is sensitive to resolution
429 (allows depper zooms at lower resolution. However it fails
430 for rotations of exactly 90 degrees, so we added a safety net
431 by using the magnification. */
432 if(++tries < 2) /* for safety */
433 {
434 static FCODE err[] = {"precision-detection error"};
435 if(tries > 1) stopmsg(0, err);
436 /* Previously there were four tests of distortions in the
437 zoom box used to detect precision limitations. In some
438 cases of rotated/skewed zoom boxs, this causes the algorithm
439 to bail out to arbitrary precision too soon. The logic
440 now only tests the larger of the two deltas in an attempt
441 to repair this bug. This should never make the transition
442 to arbitrary precision sooner, but always later.*/
443 if(fabs(xxmax-xx3rd) > fabs(xx3rd-xxmin))
444 {
445 testx_exact = xxmax-xx3rd;
446 testx_try = dx0-xxmin;
447 }
448 else
449 {
450 testx_exact = xx3rd-xxmin;
451 testx_try = dx1;
452 }
453 if(fabs(yy3rd-yymax) > fabs(yymin-yy3rd))
454 {
455 testy_exact = yy3rd-yymax;
456 testy_try = dy0-yymax;
457 }
458 else
459 {
460 testy_exact = yymin-yy3rd;
461 testy_try = dy1;
462 }
463 if(ratio_bad(testx_try,testx_exact) ||
464 ratio_bad(testy_try,testy_exact))
465 {
466 if(curfractalspecific->flags & BF_MATH)
467 {
468 fractal_floattobf();
469 goto init_restart;
470 }
471 goto expand_retry;
472 } /* end if ratio_bad etc. */
473 } /* end if tries < 2 */
474 } /* end if bf_math == 0 */
475
476 /* if long double available, this is more accurate */
477 fill_dx_array(); /* fill up the x, y grids */
478
479 /* re-set corners to match reality */
480 xxmax = (double)(xxmin + (xdots-1)*delxx + (ydots-1)*delxx2);
481 yymin = (double)(yymax - (ydots-1)*delyy - (xdots-1)*delyy2);
482 xx3rd = (double)(xxmin + (ydots-1)*delxx2);
483 yy3rd = (double)(yymax - (ydots-1)*delyy);
484
485 } /* end else */
486 } /* end if not plasma */
487
488 /* for periodicity close-enough, and for unity: */
489 /* min(max(delx,delx2),max(dely,dely2) */
490 ddelmin = fabs((double)delxx);
491 if (fabs((double)delxx2) > ddelmin)
492 ddelmin = fabs((double)delxx2);
493 if (fabs((double)delyy) > fabs((double)delyy2))
494 {
495 if (fabs((double)delyy) < ddelmin)
496 ddelmin = fabs((double)delyy);
497 }
498 else if (fabs((double)delyy2) < ddelmin)
499 ddelmin = fabs((double)delyy2);
500 delmin = fudgetolong(ddelmin);
501
502 /* calculate factors which plot real values to screen co-ords */
503 /* calcfrac.c plot_orbit routines have comments about this */
504 ftemp = (double)((0.0-delyy2) * delxx2 * dxsize * dysize
505 - (xxmax-xx3rd) * (yy3rd-yymax));
506 if(ftemp != 0)
507 {
508 plotmx1 = (double)(delxx2 * dxsize * dysize / ftemp);
509 plotmx2 = (yy3rd-yymax) * dxsize / ftemp;
510 plotmy1 = (double)((0.0-delyy2) * dxsize * dysize / ftemp);
511 plotmy2 = (xxmax-xx3rd) * dysize / ftemp;
512 }
513 if(bf_math == 0)
514 free_bf_vars();
515 else
516 {
517 /* zap all of extraseg except high area to flush out bugs */
518 /* in production version this code can be deleted */
519 char far *extra;
520 extra = (char far *)MK_FP(extraseg,0);
521 far_memset(extra,0,(unsigned int)(0x10000l-(bflength+2)*22U));
522 }
523 }
524
525 #ifdef _MSC_VER
526 #if _MSC_VER == 800
527 #pragma optimize( "", on ) /* restore optimization options */
528 #endif
529 #endif
530
531 static long _fastcall fudgetolong(double d)
532 {
533 if ((d *= fudge) > 0) d += 0.5;
534 else d -= 0.5;
535 return (long)d;
536 }
537
538 static double _fastcall fudgetodouble(long l)
539 {
540 char buf[30];
541 double d;
542 sprintf(buf,"%.9g",(double)l / fudge);
543 #ifndef XFRACT
544 sscanf(buf,"%lg",&d);
547 #endif
548 return d;
549 }
550
551 void adjust_cornerbf(void)
552 {
553 /* make edges very near vert/horiz exact, to ditch rounding errs and */
554 /* to avoid problems when delta per axis makes too large a ratio */
555 double ftemp;
556 double Xmagfactor, Rotation, Skew;
557 LDBL Magnification;
558
559 bf_t bftemp, bftemp2;
560 bf_t btmp1;
561 int saved; saved = save_stack();
562 bftemp = alloc_stack(rbflength+2);
563 bftemp2 = alloc_stack(rbflength+2);
564 btmp1 = alloc_stack(rbflength+2);
565
566 /* While we're at it, let's adjust the Xmagfactor as well */
567 /* use bftemp, bftemp2 as bfXctr, bfYctr */
568 cvtcentermagbf(bftemp, bftemp2, &Magnification, &Xmagfactor, &Rotation, &Skew);
569 ftemp = fabs(Xmagfactor);
570 if (ftemp != 1 && ftemp >= (1-aspectdrift) && ftemp <= (1+aspectdrift))
571 {
572 Xmagfactor = sign(Xmagfactor);
573 cvtcornersbf(bftemp, bftemp2, Magnification, Xmagfactor, Rotation, Skew);
574 }
575
576 /* ftemp=fabs(xx3rd-xxmin); */
577 abs_a_bf(sub_bf(bftemp,bfx3rd,bfxmin));
578
579 /* ftemp2=fabs(xxmax-xx3rd);*/
580 abs_a_bf(sub_bf(bftemp2,bfxmax,bfx3rd));
581
582 /* if( (ftemp=fabs(xx3rd-xxmin)) < (ftemp2=fabs(xxmax-xx3rd)) ) */
583 if(cmp_bf(bftemp,bftemp2) < 0)
584 {
585 /* if (ftemp*10000 < ftemp2 && yy3rd != yymax) */
586 if (cmp_bf(mult_bf_int(btmp1,bftemp,10000),bftemp2) < 0
587 && cmp_bf(bfy3rd,bfymax) != 0 )
588 /* xx3rd = xxmin; */
589 copy_bf(bfx3rd, bfxmin);
590 }
591
592 /* else if (ftemp2*10000 < ftemp && yy3rd != yymin) */
593 if (cmp_bf(mult_bf_int(btmp1,bftemp2,10000),bftemp) < 0
594 && cmp_bf(bfy3rd,bfymin) != 0 )
595 /* xx3rd = xxmax; */
596 copy_bf(bfx3rd, bfxmax);
597
598 /* ftemp=fabs(yy3rd-yymin); */
599 abs_a_bf(sub_bf(bftemp,bfy3rd,bfymin));
600
601 /* ftemp2=fabs(yymax-yy3rd); */
602 abs_a_bf(sub_bf(bftemp2,bfymax,bfy3rd));
603
604 /* if( (ftemp=fabs(yy3rd-yymin)) < (ftemp2=fabs(yymax-yy3rd)) ) */
605 if(cmp_bf(bftemp,bftemp2) < 0)
606 {
607 /* if (ftemp*10000 < ftemp2 && xx3rd != xxmax) */
608 if (cmp_bf(mult_bf_int(btmp1,bftemp,10000),bftemp2) < 0
609 && cmp_bf(bfx3rd,bfxmax) != 0 )
610 /* yy3rd = yymin; */
611 copy_bf(bfy3rd, bfymin);
612 }
613
614 /* else if (ftemp2*10000 < ftemp && xx3rd != xxmin) */
615 if (cmp_bf(mult_bf_int(btmp1,bftemp2,10000),bftemp) < 0
616 && cmp_bf(bfx3rd,bfxmin) != 0 )
617 /* yy3rd = yymax; */
618 copy_bf(bfy3rd, bfymax);
619
620
621 restore_stack(saved);
622 }
623
624 void adjust_corner(void)
625 {
626 /* make edges very near vert/horiz exact, to ditch rounding errs and */
627 /* to avoid problems when delta per axis makes too large a ratio */
628 double ftemp,ftemp2;
629 double Xctr, Yctr, Xmagfactor, Rotation, Skew;
630 LDBL Magnification;
631
632 if(!integerfractal)
633 {
634 /* While we're at it, let's adjust the Xmagfactor as well */
635 cvtcentermag(&Xctr, &Yctr, &Magnification, &Xmagfactor, &Rotation, &Skew);
636 ftemp = fabs(Xmagfactor);
637 if (ftemp != 1 && ftemp >= (1-aspectdrift) && ftemp <= (1+aspectdrift))
638 {
639 Xmagfactor = sign(Xmagfactor);
640 cvtcorners(Xctr, Yctr, Magnification, Xmagfactor, Rotation, Skew);
641 }
642 }
643
644 if( (ftemp=fabs(xx3rd-xxmin)) < (ftemp2=fabs(xxmax-xx3rd)) ) {
645 if (ftemp*10000 < ftemp2 && yy3rd != yymax)
646 xx3rd = xxmin;
647 }
648
649 if (ftemp2*10000 < ftemp && yy3rd != yymin)
650 xx3rd = xxmax;
651
652 if( (ftemp=fabs(yy3rd-yymin)) < (ftemp2=fabs(yymax-yy3rd)) ) {
653 if (ftemp*10000 < ftemp2 && xx3rd != xxmax)
654 yy3rd = yymin;
655 }
656
657 if (ftemp2*10000 < ftemp && xx3rd != xxmin)
658 yy3rd = yymax;
659
660 }
661
662 static void _fastcall adjust_to_limitsbf(double expand)
663 {
664 LDBL limit;
665 bf_t bcornerx[4],bcornery[4];
666 bf_t blowx,bhighx,blowy,bhighy,blimit,bftemp;
667 bf_t bcenterx,bcentery,badjx,badjy,btmp1,btmp2;
668 bf_t bexpand;
669 int i;
670 int saved; saved = save_stack();
671 bcornerx[0] = alloc_stack(rbflength+2);
672 bcornerx[1] = alloc_stack(rbflength+2);
673 bcornerx[2] = alloc_stack(rbflength+2);
674 bcornerx[3] = alloc_stack(rbflength+2);
675 bcornery[0] = alloc_stack(rbflength+2);
676 bcornery[1] = alloc_stack(rbflength+2);
677 bcornery[2] = alloc_stack(rbflength+2);
678 bcornery[3] = alloc_stack(rbflength+2);
679 blowx = alloc_stack(rbflength+2);
680 bhighx = alloc_stack(rbflength+2);
681 blowy = alloc_stack(rbflength+2);
682 bhighy = alloc_stack(rbflength+2);
683 blimit = alloc_stack(rbflength+2);
684 bftemp = alloc_stack(rbflength+2);
685 bcenterx = alloc_stack(rbflength+2);
686 bcentery = alloc_stack(rbflength+2);
687 badjx = alloc_stack(rbflength+2);
688 badjy = alloc_stack(rbflength+2);
689 btmp1 = alloc_stack(rbflength+2);
690 btmp2 = alloc_stack(rbflength+2);
691 bexpand = alloc_stack(rbflength+2);
692
693 limit = 32767.99;
694
695 /* if (bitshift >= 24) limit = 31.99;
696 if (bitshift >= 29) limit = 3.99; */
697 floattobf(blimit,limit);
698 floattobf(bexpand,expand);
699
700 add_bf(bcenterx,bfxmin,bfxmax);
701 half_a_bf(bcenterx);
702
703 /* centery = (yymin+yymax)/2; */
704 add_bf(bcentery,bfymin,bfymax);
705 half_a_bf(bcentery);
706
707 /* if (xxmin == centerx) { */
708 if (cmp_bf(bfxmin,bcenterx)==0) { /* ohoh, infinitely thin, fix it */
709 smallest_add_bf(bfxmax);
710 /* bfxmin -= bfxmax-centerx; */
711 sub_a_bf(bfxmin,sub_bf(btmp1,bfxmax,bcenterx));
712 }
713
714 /* if (bfymin == centery) */
715 if (cmp_bf(bfymin,bcentery)==0) {
716 smallest_add_bf(bfymax);
717 /* bfymin -= bfymax-centery; */
718 sub_a_bf(bfymin,sub_bf(btmp1,bfymax,bcentery));
719 }
720
721 /* if (bfx3rd == centerx) */
722 if (cmp_bf(bfx3rd,bcenterx)==0)
723 smallest_add_bf(bfx3rd);
724
725 /* if (bfy3rd == centery) */
726 if (cmp_bf(bfy3rd,bcentery)==0)
727 smallest_add_bf(bfy3rd);
728
729 /* setup array for easier manipulation */
730 /* cornerx[0] = xxmin; */
731 copy_bf(bcornerx[0],bfxmin);
732
733 /* cornerx[1] = xxmax; */
734 copy_bf(bcornerx[1],bfxmax);
735
736 /* cornerx[2] = xx3rd; */
737 copy_bf(bcornerx[2],bfx3rd);
738
739 /* cornerx[3] = xxmin+(xxmax-xx3rd); */
740 sub_bf(bcornerx[3],bfxmax,bfx3rd);
741 add_a_bf(bcornerx[3],bfxmin);
742
743 /* cornery[0] = yymax; */
744 copy_bf(bcornery[0],bfymax);
745
746 /* cornery[1] = yymin; */
747 copy_bf(bcornery[1],bfymin);
748
749 /* cornery[2] = yy3rd; */
750 copy_bf(bcornery[2],bfy3rd);
751
752 /* cornery[3] = yymin+(yymax-yy3rd); */
753 sub_bf(bcornery[3],bfymax,bfy3rd);
754 add_a_bf(bcornery[3],bfymin);
755
756 /* if caller wants image size adjusted, do that first */
757 if (expand != 1.0)
758 {
759 for (i=0; i<4; ++i) {
760 /* cornerx[i] = centerx + (cornerx[i]-centerx)*expand; */
761 sub_bf(btmp1,bcornerx[i],bcenterx);
762 mult_bf(bcornerx[i],btmp1,bexpand);
763 add_a_bf(bcornerx[i],bcenterx);
764
765 /* cornery[i] = centery + (cornery[i]-centery)*expand; */
766 sub_bf(btmp1,bcornery[i],bcentery);
767 mult_bf(bcornery[i],btmp1,bexpand);
768 add_a_bf(bcornery[i],bcentery);
769 }
770 }
771
772 /* get min/max x/y values */
773 /* lowx = highx = cornerx[0]; */
774 copy_bf(blowx,bcornerx[0]); copy_bf(bhighx,bcornerx[0]);
775
776 /* lowy = highy = cornery[0]; */
777 copy_bf(blowy,bcornery[0]); copy_bf(bhighy,bcornery[0]);
778
779 for (i=1; i<4; ++i) {
780 /* if (cornerx[i] < lowx) lowx = cornerx[i]; */
781 if (cmp_bf(bcornerx[i],blowx) < 0) copy_bf(blowx,bcornerx[i]);
782
783 /* if (cornerx[i] > highx) highx = cornerx[i]; */
784 if (cmp_bf(bcornerx[i],bhighx) > 0) copy_bf(bhighx,bcornerx[i]);
785
786 /* if (cornery[i] < lowy) lowy = cornery[i]; */
787 if (cmp_bf(bcornery[i],blowy) < 0) copy_bf(blowy,bcornery[i]);
788
789 /* if (cornery[i] > highy) highy = cornery[i]; */
790 if (cmp_bf(bcornery[i],bhighy) > 0) copy_bf(bhighy,bcornery[i]);
791 }
792
793 /* if image is too large, downsize it maintaining center */
794 /* ftemp = highx-lowx; */
795 sub_bf(bftemp,bhighx,blowx);
796
797 /* if (highy-lowy > ftemp) ftemp = highy-lowy; */
798 if (cmp_bf(sub_bf(btmp1,bhighy,blowy),bftemp) > 0) copy_bf(bftemp,btmp1);
799
800 /* if image is too large, downsize it maintaining center */
801
802 floattobf(btmp1,limit*2.0);
803 copy_bf(btmp2,bftemp);
804 div_bf(bftemp,btmp1,btmp2);
805 floattobf(btmp1,1.0);
806 if (cmp_bf(bftemp,btmp1) < 0)
807 for (i=0; i<4; ++i) {
808 /* cornerx[i] = centerx + (cornerx[i]-centerx)*ftemp; */
809 sub_bf(btmp1,bcornerx[i],bcenterx);
810 mult_bf(bcornerx[i],btmp1,bftemp);
811 add_a_bf(bcornerx[i],bcenterx);
812
813 /* cornery[i] = centery + (cornery[i]-centery)*ftemp; */
814 sub_bf(btmp1,bcornery[i],bcentery);
815 mult_bf(bcornery[i],btmp1,bftemp);
816 add_a_bf(bcornery[i],bcentery);
817 }
818
819 /* if any corner has x or y past limit, move the image */
820 /* adjx = adjy = 0; */
821 clear_bf(badjx); clear_bf(badjy);
822
823 for (i=0; i<4; ++i) {
824 /* if (cornerx[i] > limit && (ftemp = cornerx[i] - limit) > adjx)
825 adjx = ftemp; */
826 if (cmp_bf(bcornerx[i],blimit) > 0 &&
827 cmp_bf(sub_bf(bftemp,bcornerx[i],blimit),badjx) > 0)
828 copy_bf(badjx,bftemp);
829
830 /* if (cornerx[i] < 0.0-limit && (ftemp = cornerx[i] + limit) < adjx)
831 adjx = ftemp; */
832 if (cmp_bf(bcornerx[i],neg_bf(btmp1,blimit)) < 0 &&
833 cmp_bf(add_bf(bftemp,bcornerx[i],blimit),badjx) < 0)
834 copy_bf(badjx,bftemp);
835
836 /* if (cornery[i] > limit && (ftemp = cornery[i] - limit) > adjy)
837 adjy = ftemp; */
838 if (cmp_bf(bcornery[i],blimit) > 0 &&
839 cmp_bf(sub_bf(bftemp,bcornery[i],blimit),badjy) > 0)
840 copy_bf(badjy,bftemp);
841
842 /* if (cornery[i] < 0.0-limit && (ftemp = cornery[i] + limit) < adjy)
843 adjy = ftemp; */
844 if (cmp_bf(bcornery[i],neg_bf(btmp1,blimit)) < 0 &&
845 cmp_bf(add_bf(bftemp,bcornery[i],blimit),badjy) < 0)
846 copy_bf(badjy,bftemp);
847 }
848
849 /* if (calc_status == 2 && (adjx != 0 || adjy != 0) && (zwidth == 1.0))
850 calc_status = 0; */
851 if (calc_status == 2 && (is_bf_not_zero(badjx)|| is_bf_not_zero(badjy)) && (zwidth == 1.0))
852 calc_status = 0;
853
854 /* xxmin = cornerx[0] - adjx; */
855 sub_bf(bfxmin,bcornerx[0],badjx);
856 /* xxmax = cornerx[1] - adjx; */
857 sub_bf(bfxmax,bcornerx[1],badjx);
858 /* xx3rd = cornerx[2] - adjx; */
859 sub_bf(bfx3rd,bcornerx[2],badjx);
860 /* yymax = cornery[0] - adjy; */
861 sub_bf(bfymax,bcornery[0],badjy);
862 /* yymin = cornery[1] - adjy; */
863 sub_bf(bfymin,bcornery[1],badjy);
864 /* yy3rd = cornery[2] - adjy; */
865 sub_bf(bfy3rd,bcornery[2],badjy);
866
867 adjust_cornerbf(); /* make 3rd corner exact if very near other co-ords */
868 restore_stack(saved);
869 }
870
871 static void _fastcall adjust_to_limits(double expand)
872 {
873 double cornerx[4],cornery[4];
874 double lowx,highx,lowy,highy,limit,ftemp;
875 double centerx,centery,adjx,adjy;
876 int i;
877
878 limit = 32767.99;
879
880 if (integerfractal) {
881 if (save_release > 1940) /* let user reproduce old GIF's and PAR's */
882 limit = 1023.99;
883 if (bitshift >= 24) limit = 31.99;
884 if (bitshift >= 29) limit = 3.99;
885 }
886
887 centerx = (xxmin+xxmax)/2;
888 centery = (yymin+yymax)/2;
889
890 if (xxmin == centerx) { /* ohoh, infinitely thin, fix it */
891 smallest_add(&xxmax);
892 xxmin -= xxmax-centerx;
893 }
894
895 if (yymin == centery) {
896 smallest_add(&yymax);
897 yymin -= yymax-centery;
898 }
899
900 if (xx3rd == centerx)
901 smallest_add(&xx3rd);
902
903 if (yy3rd == centery)
904 smallest_add(&yy3rd);
905
906 /* setup array for easier manipulation */
907 cornerx[0] = xxmin;
908 cornerx[1] = xxmax;
909 cornerx[2] = xx3rd;
910 cornerx[3] = xxmin+(xxmax-xx3rd);
911
912 cornery[0] = yymax;
913 cornery[1] = yymin;
914 cornery[2] = yy3rd;
915 cornery[3] = yymin+(yymax-yy3rd);
916
917 /* if caller wants image size adjusted, do that first */
918 if (expand != 1.0)
919 {
920 for (i=0; i<4; ++i) {
921 cornerx[i] = centerx + (cornerx[i]-centerx)*expand;
922 cornery[i] = centery + (cornery[i]-centery)*expand;
923 }
924 }
925 /* get min/max x/y values */
926 lowx = highx = cornerx[0];
927 lowy = highy = cornery[0];
928
929 for (i=1; i<4; ++i) {
930 if (cornerx[i] < lowx) lowx = cornerx[i];
931 if (cornerx[i] > highx) highx = cornerx[i];
932 if (cornery[i] < lowy) lowy = cornery[i];
933 if (cornery[i] > highy) highy = cornery[i];
934 }
935
936 /* if image is too large, downsize it maintaining center */
937 ftemp = highx-lowx;
938
939 if (highy-lowy > ftemp) ftemp = highy-lowy;
940
941 /* if image is too large, downsize it maintaining center */
942 if ((ftemp = limit*2/ftemp) < 1.0) {
943 for (i=0; i<4; ++i) {
944 cornerx[i] = centerx + (cornerx[i]-centerx)*ftemp;
945 cornery[i] = centery + (cornery[i]-centery)*ftemp;
946 }
947 }
948
949 /* if any corner has x or y past limit, move the image */
950 adjx = adjy = 0;
951
952 for (i=0; i<4; ++i) {
953 if (cornerx[i] > limit && (ftemp = cornerx[i] - limit) > adjx)
954 adjx = ftemp;
955 if (cornerx[i] < 0.0-limit && (ftemp = cornerx[i] + limit) < adjx)
956 adjx = ftemp;
957 if (cornery[i] > limit && (ftemp = cornery[i] - limit) > adjy)
958 adjy = ftemp;
959 if (cornery[i] < 0.0-limit && (ftemp = cornery[i] + limit) < adjy)
960 adjy = ftemp;
961 }
962 if (calc_status == 2 && (adjx != 0 || adjy != 0) && (zwidth == 1.0))
963 calc_status = 0;
964 xxmin = cornerx[0] - adjx;
965 xxmax = cornerx[1] - adjx;
966 xx3rd = cornerx[2] - adjx;
967 yymax = cornery[0] - adjy;
968 yymin = cornery[1] - adjy;
969 yy3rd = cornery[2] - adjy;
970
971 adjust_corner(); /* make 3rd corner exact if very near other co-ords */
972 }
973
974 static void _fastcall smallest_add(double *num)
975 {
976 *num += *num * 5.0e-16;
977 }
978
979 static void _fastcall smallest_add_bf(bf_t num)
980 {
981 bf_t btmp1;
982 int saved; saved = save_stack();
983 btmp1 = alloc_stack(bflength+2);
984 mult_bf(btmp1,floattobf(btmp1, 5.0e-16),num);
985 add_a_bf(num,btmp1);
986 restore_stack(saved);
987 }
988
989 static int _fastcall ratio_bad(double actual, double desired)
990 {
991 double ftemp, tol;
992 if(integerfractal)
993 tol = math_tol[0];
994 else
995 tol = math_tol[1];
996 if(tol <= 0.0)
997 return(1);
998 else if(tol >= 1.0)
999 return(0);
1000 ftemp = 0;
1001 if (desired != 0 && debugflag != 3400)
1002 ftemp = actual / desired;
1003 if (desired != 0 && debugflag != 3400)
1004 if ((ftemp = actual / desired) < (1.0-tol) || ftemp > (1.0+tol))
1005 return(1);
1006 return(0);
1007 }
1008
1009
1010 /* Save/resume stuff:
1011
1012 Engines which aren't resumable can simply ignore all this.
1013
1014 Before calling the (per_image,calctype) routines (engine), calcfract sets:
1015 "resuming" to 0 if new image, nonzero if resuming a partially done image
1016 "calc_status" to 1
1017 If an engine is interrupted and wants to be able to resume it must:
1018 store whatever status info it needs to be able to resume later
1019 set calc_status to 2 and return
1020 If subsequently called with resuming!=0, the engine must restore status
1021 info and continue from where it left off.
1022
1023 Since the info required for resume can get rather large for some types,
1024 it is not stored directly in save_info. Instead, memory is dynamically
1025 allocated as required, and stored in .fra files as a separate block.
1026 To save info for later resume, an engine routine can use:
1027 alloc_resume(maxsize,version)
1028 Maxsize must be >= max bytes subsequently saved + 2; over-allocation
1029 is harmless except for possibility of insufficient mem available;
1030 undersize is not checked and probably causes serious misbehaviour.
1031 Version is an arbitrary number so that subsequent revisions of the
1032 engine can be made backward compatible.
1033 Alloc_resume sets calc_status to 2 (resumable) if it succeeds; to 3
1034 if it cannot allocate memory (and issues warning to user).
1035 put_resume({bytes,&argument,} ... 0)
1036 Can be called as often as required to store the info.
1037 Arguments must not be far addresses.
1038 Is not protected against calls which use more space than allocated.
1039 To reload info when resuming, use:
1040 start_resume()
1041 initializes and returns version number
1042 get_resume({bytes,&argument,} ... 0)
1043 inverse of store_resume
1044 end_resume()
1045 optional, frees the memory area sooner than would happen otherwise
1046
1047 Example, save info:
1048 alloc_resume(sizeof(parmarray)+100,2);
1049 put_resume(sizeof(row),&row, sizeof(col),&col,
1050 sizeof(parmarray),parmarray, 0);
1051 restore info:
1052 vsn = start_resume();
1053 get_resume(sizeof(row),&row, sizeof(col),&col, 0);
1054 if (vsn >= 2)
1055 get_resume(sizeof(parmarray),parmarray,0);
1056 end_resume();
1057
1058 Engines which allocate a large far memory chunk of their own might
1059 directly set resume_info, resume_len, calc_status to avoid doubling
1060 transient memory needs by using these routines.
1061
1062 StandardFractal, calcmand, solidguess, and bound_trace_main are a related
1063 set of engines for escape-time fractals. They use a common worklist
1064 structure for save/resume. Fractals using these must specify calcmand
1065 or StandardFractal as the engine in fractalspecificinfo.
1066 Other engines don't get btm nor ssg, don't get off-axis symmetry nor
1067 panning (the worklist stuff), and are on their own for save/resume.
1068
1069 */
1070
1071 #ifndef USE_VARARGS
1072 int put_resume(int len, ...)
1076 #endif
1077 {
1078 va_list arg_marker; /* variable arg list */
1079 BYTE *source_ptr;
1080 #ifdef USE_VARARGS
1082 #endif
1083
1084 if (resume_info == 0)
1085 return(-1);
1086 #ifndef USE_VARARGS
1087 va_start(arg_marker,len);
1091 #endif
1092 while (len)
1093 {
1094 source_ptr = (BYTE *)va_arg(arg_marker,char *);
1095 /* far_memcpy(resume_info+resume_len,source_ptr,len); */
1096 MoveToMemory(source_ptr,(U16)1,(long)len,resume_len,resume_info);
1097 resume_len += len;
1098 len = va_arg(arg_marker,int);
1099 }
1100 va_end(arg_marker);
1101 return(0);
1102 }
1103
1104 int alloc_resume(int alloclen, int version)
1105 { /* WARNING! if alloclen > 4096B, problems may occur with GIF save/restore */
1106 if (resume_info != 0) /* free the prior area if there is one */
1107 MemoryRelease(resume_info);
1108 if ((resume_info = MemoryAlloc((U16)sizeof(alloclen), (long)alloclen, FARMEM)) == 0)
1109 {
1110 static FCODE msg[] = {"\
1111 Warning - insufficient free memory to save status.\n\
1112 You will not be able to resume calculating this image."};
1113 stopmsg(0,msg);
1114 calc_status = 3;
1115 return(-1);
1116 }
1117 resume_len = 0;
1118 put_resume(sizeof(version),&version,0);
1119 calc_status = 2;
1120 return(0);
1121 }
1122
1123 #ifndef USE_VARARGS
1124 int get_resume(int len, ...)
1128 #endif
1129 {
1130 va_list arg_marker; /* variable arg list */
1131 BYTE *dest_ptr;
1132 #ifdef USE_VARARGS
1134 #endif
1135
1136 if (resume_info == 0)
1137 return(-1);
1138 #ifndef USE_VARARGS
1139 va_start(arg_marker,len);
1143 #endif
1144 while (len)
1145 {
1146 dest_ptr = (BYTE *)va_arg(arg_marker,char *);
1147 /* far_memcpy(dest_ptr,resume_info+resume_offset,len); */
1148 MoveFromMemory(dest_ptr,(U16)1,(long)len,resume_offset,resume_info);
1149 resume_offset += len;
1150 len = va_arg(arg_marker,int);
1151 }
1152 va_end(arg_marker);
1153 return(0);
1154 }
1155
1156 int start_resume(void)
1157 {
1158 int version;
1159 if (resume_info == 0)
1160 return(-1);
1161 resume_offset = 0;
1162 get_resume(sizeof(version),&version,0);
1163 return(version);
1164 }
1165
1166 void end_resume(void)
1167 {
1168 if (resume_info != 0) /* free the prior area if there is one */
1169 {
1170 MemoryRelease(resume_info);
1171 resume_info = 0;
1172 }
1173 }
1174
1175
1176 /* Showing orbit requires converting real co-ords to screen co-ords.
1177 Define:
1178 Xs == xxmax-xx3rd Ys == yy3rd-yymax
1179 W == xdots-1 D == ydots-1
1180 We know that:
1181 realx == lx0[col] + lx1[row]
1182 realy == ly0[row] + ly1[col]
1183 lx0[col] == (col/width) * Xs + xxmin
1184 lx1[row] == row * delxx
1185 ly0[row] == (row/D) * Ys + yymax
1186 ly1[col] == col * (0-delyy)
1187 so:
1188 realx == (col/W) * Xs + xxmin + row * delxx
1189 realy == (row/D) * Ys + yymax + col * (0-delyy)
1190 and therefore:
1191 row == (realx-xxmin - (col/W)*Xs) / Xv (1)
1192 col == (realy-yymax - (row/D)*Ys) / Yv (2)
1193 substitute (2) into (1) and solve for row:
1194 row == ((realx-xxmin)*(0-delyy2)*W*D - (realy-yymax)*Xs*D)
1195 / ((0-delyy2)*W*delxx2*D-Ys*Xs)
1196 */
1197
1198 /* sleep N * a tenth of a millisecond */
1199
1200 void sleepms_old(long ms)
1201 {
1202 static long scalems = 0L;
1203 int savehelpmode,savetabmode;
1204 struct timebx t1,t2;
1205 #define SLEEPINIT 250 /* milliseconds for calibration */
1206 savetabmode = tabmode;
1207 savehelpmode = helpmode;
1208 tabmode = 0;
1209 helpmode = -1;
1210 if(scalems==0L) /* calibrate */
1211 {
1212 /* selects a value of scalems that makes the units
1213 10000 per sec independent of CPU speed */
1214 int i,elapsed;
1215 scalems = 1L;
1216 if(keypressed()) /* check at start, hope to get start of timeslice */
1217 goto sleepexit;
1218 /* calibrate, assume slow computer first */
1219 showtempmsg("Calibrating timer");
1220 do
1221 {
1222 scalems *= 2;
1223 ftimex(&t2);
1224 do { /* wait for the start of a new tick */
1225 ftimex(&t1);
1226 }
1227 while (t2.time == t1.time && t2.millitm == t1.millitm);
1228 sleepms_old(10L * SLEEPINIT); /* about 1/4 sec */
1229 ftimex(&t2);
1230 if(keypressed()) {
1231 scalems = 0L;
1232 cleartempmsg();
1233 goto sleepexit;
1234 }
1235 }
1236 while ((elapsed = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm)
1237 < SLEEPINIT);
1238 /* once more to see if faster (eg multi-tasking) */
1239 do { /* wait for the start of a new tick */
1240 ftimex(&t1);
1241 }
1242 while (t2.time == t1.time && t2.millitm == t1.millitm);
1243 sleepms_old(10L * SLEEPINIT);
1244 ftimex(&t2);
1245 if ((i = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm) < elapsed)
1246 elapsed = (i == 0) ? 1 : i;
1247 scalems = (long)((float)SLEEPINIT/(float)(elapsed) * scalems);
1248 cleartempmsg();
1249 }
1250 if(ms > 10L * SLEEPINIT) { /* using ftime is probably more accurate */
1251 ms /= 10;
1252 ftimex(&t1);
1253 for(;;) {
1254 if(keypressed()) break;
1255 ftimex(&t2);
1256 if ((long)((t2.time-t1.time)*1000 + t2.millitm-t1.millitm) >= ms) break;
1257 }
1258 }
1259 else
1260 if(!keypressed()) {
1261 ms *= scalems;
1262 while(ms-- >= 0);
1263 }
1264 sleepexit:
1265 tabmode = savetabmode;
1266 helpmode = savehelpmode;
1267 }
1268
1269 static void sleepms_new(long ms)
1270 {
1271 uclock_t next_time;
1272 uclock_t now = usec_clock();
1273 next_time = now + ms*100;
1274 while ((now = usec_clock()) < next_time)
1275 if(keypressed()) break;
1276 }
1277
1278 void sleepms(long ms)
1279 {
1280 if(debugflag == 4020)
1281 sleepms_old(ms);
1282 else
1283 sleepms_new(ms);
1284 }
1285
1286 /*
1287 * wait until wait_time microseconds from the
1288 * last call has elapsed.
1289 */
1290 #define MAX_INDEX 2
1291 static uclock_t next_time[MAX_INDEX];
1292 void wait_until(int index, uclock_t wait_time)
1293 {
1294 if(debugflag == 4020)
1295 sleepms_old(wait_time);
1296 else
1297 {
1298 uclock_t now;
1299 while ( (now = usec_clock()) < next_time[index])
1300 if(keypressed()) break;
1301 next_time[index] = now + wait_time*100; /* wait until this time next call */
1302 }
1303 }
1304
1305 void reset_clock(void)
1306 {
1307 int i;
1308 restart_uclock();
1309 for(i=0;i<MAX_INDEX;i++)
1310 next_time[i] = 0;
1311 }
1312
1313 #define LOG2 (float)0.693147180
1314 #define LOG32 (float)3.465735902
1315
1316 static FILE *snd_fp = NULL;
1317
1318 /* open sound file */
1319 int snd_open(void)
1320 {
1321 static char soundname[] = {"sound001.txt"};
1322 if((orbitsave&2) != 0 && snd_fp == NULL)
1323 {
1324 if((snd_fp = fopen(soundname,"w"))==NULL)
1325 {
1326 static FCODE msg[] = {"Can't open SOUND*.TXT"};
1327 stopmsg(0,msg);
1328 }
1329 else
1330 {
1331 updatesavename(soundname);
1332 }
1333 }
1334 return(snd_fp != NULL);
1335 }
1336
1337 /* This routine plays a tone in the speaker and optionally writes a file
1338 if the orbitsave variable is turned on */
1339 void w_snd(int tone)
1340 {
1341 if((orbitsave&2) != 0)
1342 {
1343 if(snd_open())
1344 fprintf(snd_fp,"%-d\n",tone);
1345 }
1346 taborhelp = 0;
1347 if(!keypressed()) { /* keypressed calls mute() if TAB or F1 pressed */
1348 /* must not then call soundoff(), else indexes out of synch */
1349 /* if(20 < tone && tone < 15000) better limits? */
1350 /* if(10 < tone && tone < 5000) better limits? */
1351 if(soundon(tone)) {
1352 wait_until(0,orbit_delay);
1353 if(!taborhelp) /* kludge because wait_until() calls keypressed */
1354 soundoff();
1355 }
1356 }
1357 }
1358
1359 void snd_time_write(void)
1360 {
1361 if(snd_open())
1362 {
1363 fprintf(snd_fp,"time=%-ld\n",(long)clock()*1000/CLK_TCK);
1364 }
1365 }
1366
1367 void close_snd(void)
1368 {
1369 if(snd_fp)
1370 fclose(snd_fp);
1371 snd_fp = NULL;
1372 }
1373
1374 static void _fastcall plotdorbit(double dx, double dy, int color)
1375 {
1376 int i, j, c;
1377 int save_sxoffs,save_syoffs;
1378 if (orbit_ptr >= 1500) return;
1379 i = (int)(dy * plotmx1 - dx * plotmx2); i += sxoffs;
1380 if (i < 0 || i >= sxdots) return;
1381 j = (int)(dx * plotmy1 - dy * plotmy2); j += syoffs;
1382 if (j < 0 || j >= sydots) return;
1383 save_sxoffs = sxoffs;
1384 save_syoffs = syoffs;
1385 sxoffs = syoffs = 0;
1386 /* save orbit value */
1387 if(color == -1)
1388 {
1389 *(save_orbit + orbit_ptr++) = i;
1390 *(save_orbit + orbit_ptr++) = j;
1391 *(save_orbit + orbit_ptr++) = c = getcolor(i,j);
1392 putcolor(i,j,c^orbit_color);
1393 }
1394 else
1395 putcolor(i,j,color);
1396 sxoffs = save_sxoffs;
1397 syoffs = save_syoffs;
1398 if(debugflag == 4030) {
1399 if((soundflag&7) == 2) /* sound = x */
1400 w_snd((int)(i*1000/xdots+basehertz));
1401 else if((soundflag&7) > 2) /* sound = y or z */
1402 w_snd((int)(j*1000/ydots+basehertz));
1403 else if(orbit_delay > 0)
1404 {
1405 wait_until(0,orbit_delay);
1406 }
1407 }
1408 else {
1409 if((soundflag&7) == 2) /* sound = x */
1410 w_snd((int)(i+basehertz));
1411 else if((soundflag&7) == 3) /* sound = y */
1412 w_snd((int)(j+basehertz));
1413 else if((soundflag&7) == 4) /* sound = z */
1414 w_snd((int)(i+j+basehertz));
1415 else if(orbit_delay > 0)
1416 {
1417 wait_until(0,orbit_delay);
1418 }
1419 }
1420
1421 /* placing sleepms here delays each dot */
1422 }
1423
1424 void iplot_orbit(long ix, long iy, int color)
1425 {
1426 plotdorbit((double)ix/fudge-xxmin,(double)iy/fudge-yymax,color);
1427 }
1428
1429 void plot_orbit(double real,double imag,int color)
1430 {
1431 plotdorbit(real-xxmin,imag-yymax,color);
1432 }
1433
1434 void scrub_orbit(void)
1435 {
1436 int i,j,c;
1437 int save_sxoffs,save_syoffs;
1438 mute();
1439 save_sxoffs = sxoffs;
1440 save_syoffs = syoffs;
1441 sxoffs = syoffs = 0;
1442 while(orbit_ptr > 0)
1443 {
1444 c = *(save_orbit + --orbit_ptr);
1445 j = *(save_orbit + --orbit_ptr);
1446 i = *(save_orbit + --orbit_ptr);
1447 putcolor(i,j,c);
1448 }
1449 sxoffs = save_sxoffs;
1450 syoffs = save_syoffs;
1451 }
1452
1453
1454 int add_worklist(int xfrom, int xto, int xbegin,
1455 int yfrom, int yto, int ybegin,
1456 int pass, int sym)
1457 {
1458 if (num_worklist >= MAXCALCWORK)
1459 return(-1);
1460 worklist[num_worklist].xxstart = xfrom;
1461 worklist[num_worklist].xxstop = xto;
1462 worklist[num_worklist].xxbegin = xbegin;
1463 worklist[num_worklist].yystart = yfrom;
1464 worklist[num_worklist].yystop = yto;
1465 worklist[num_worklist].yybegin = ybegin;
1466 worklist[num_worklist].pass = pass;
1467 worklist[num_worklist].sym = sym;
1468 ++num_worklist;
1469 tidy_worklist();
1470 return(0);
1471 }
1472
1473 static int _fastcall combine_worklist(void) /* look for 2 entries which can freely merge */
1474 {
1475 int i,j;
1476 for (i=0; i<num_worklist; ++i)
1477 if (worklist[i].yystart == worklist[i].yybegin)
1478 for (j=i+1; j<num_worklist; ++j)
1479 if (worklist[j].sym == worklist[i].sym
1480 && worklist[j].yystart == worklist[j].yybegin
1481 && worklist[j].xxstart == worklist[j].xxbegin
1482 && worklist[i].pass == worklist[j].pass)
1483 {
1484 if ( worklist[i].xxstart == worklist[j].xxstart
1485 && worklist[i].xxbegin == worklist[j].xxbegin
1486 && worklist[i].xxstop == worklist[j].xxstop)
1487 {
1488 if (worklist[i].yystop+1 == worklist[j].yystart)
1489 {
1490 worklist[i].yystop = worklist[j].yystop;
1491 return(j);
1492 }
1493 if (worklist[j].yystop+1 == worklist[i].yystart)
1494 {
1495 worklist[i].yystart = worklist[j].yystart;
1496 worklist[i].yybegin = worklist[j].yybegin;
1497 return(j);
1498 }
1499 }
1500 if ( worklist[i].yystart == worklist[j].yystart
1501 && worklist[i].yybegin == worklist[j].yybegin
1502 && worklist[i].yystop == worklist[j].yystop)
1503 {
1504 if (worklist[i].xxstop+1 == worklist[j].xxstart)
1505 {
1506 worklist[i].xxstop = worklist[j].xxstop;
1507 return(j);
1508 }
1509 if (worklist[j].xxstop+1 == worklist[i].xxstart)
1510 {
1511 worklist[i].xxstart = worklist[j].xxstart;
1512 worklist[i].xxbegin = worklist[j].xxbegin;
1513 return(j);
1514 }
1515 }
1516 }
1517 return(0); /* nothing combined */
1518 }
1519
1520 void tidy_worklist(void) /* combine mergeable entries, resort */
1521 {
1522 int i,j;
1523 WORKLIST tempwork;
1524 while ((i=combine_worklist()) != 0)
1525 { /* merged two, delete the gone one */
1526 while (++i < num_worklist)
1527 worklist[i-1] = worklist[i];
1528 --num_worklist;
1529 }
1530 for (i=0; i<num_worklist; ++i)
1531 for (j=i+1; j<num_worklist; ++j)
1532 if (worklist[j].pass < worklist[i].pass
1533 || (worklist[j].pass == worklist[i].pass
1534 && (worklist[j].yystart < worklist[i].yystart
1535 || ( worklist[j].yystart == worklist[i].yystart
1536 && worklist[j].xxstart < worklist[i].xxstart))))
1537 { /* dumb sort, swap 2 entries to correct order */
1538 tempwork = worklist[i];
1539 worklist[i] = worklist[j];
1540 worklist[j] = tempwork;
1541 }
1542 }
1543
1544
1545 void get_julia_attractor (double real, double imag)
1546 {
1547 _LCMPLX lresult;
1548 _CMPLX result;
1549 int savper;
1550 long savmaxit;
1551 int i;
1552
1553 if (attractors == 0 && finattract == 0) /* not magnet & not requested */
1554 return;
1555
1556 if (attractors >= N_ATTR) /* space for more attractors ? */
1557 return; /* Bad luck - no room left ! */
1558
1559 savper = periodicitycheck;
1560 savmaxit = maxit;
1561 periodicitycheck = 0;
1562 old.x = real; /* prepare for f.p orbit calc */
1563 old.y = imag;
1564 tempsqrx = sqr(old.x);
1565 tempsqry = sqr(old.y);
1566
1567 lold.x = (long)real; /* prepare for int orbit calc */
1568 lold.y = (long)imag;
1569 ltempsqrx = (long)tempsqrx;
1570 ltempsqry = (long)tempsqry;
1571
1572 lold.x = lold.x << bitshift;
1573 lold.y = lold.y << bitshift;
1574 ltempsqrx = ltempsqrx << bitshift;
1575 ltempsqry = ltempsqry << bitshift;
1576
1577 if (maxit < 500) /* we're going to try at least this hard */
1578 maxit = 500;
1579 coloriter = 0;
1580 overflow = 0;
1581 while (++coloriter < maxit)
1582 if(curfractalspecific->orbitcalc() || overflow)
1583 break;
1584 if (coloriter >= maxit) /* if orbit stays in the lake */
1585 {
1586 if (integerfractal) /* remember where it went to */
1587 lresult = lnew;
1588 else
1589 result = new;
1590 for (i=0;i<10;i++) {
1591 overflow = 0;
1592 if(!curfractalspecific->orbitcalc() && !overflow) /* if it stays in the lake */
1593 { /* and doesn't move far, probably */
1594 if (integerfractal) /* found a finite attractor */
1595 {
1596 if(labs(lresult.x-lnew.x) < lclosenuff
1597 && labs(lresult.y-lnew.y) < lclosenuff)
1598 {
1599 lattr[attractors] = lnew;
1600 attrperiod[attractors] = i+1;
1601 attractors++; /* another attractor - coloured lakes ! */
1602 break;
1603 }
1604 }
1605 else
1606 {
1607 if(fabs(result.x-new.x) < closenuff
1608 && fabs(result.y-new.y) < closenuff)
1609 {
1610 attr[attractors] = new;
1611 attrperiod[attractors] = i+1;
1612 attractors++; /* another attractor - coloured lakes ! */
1613 break;
1614 }
1615 }
1616 } else {
1617 break;
1618 }
1619 }
1620 }
1621 if(attractors==0)
1622 periodicitycheck = savper;
1623 maxit = savmaxit;
1624 }
1625
1626
1627 #define maxyblk 7 /* must match calcfrac.c */
1628 #define maxxblk 202 /* must match calcfrac.c */
1629 int ssg_blocksize(void) /* used by solidguessing and by zoom panning */
1630 {
1631 int blocksize,i;
1632 /* blocksize 4 if <300 rows, 8 if 300-599, 16 if 600-1199, 32 if >=1200 */
1633 blocksize=4;
1634 i=300;
1635 while(i<=ydots)
1636 {
1637 blocksize+=blocksize;
1638 i+=i;
1639 }
1640 /* increase blocksize if prefix array not big enough */
1641 while(blocksize*(maxxblk-2)<xdots || blocksize*(maxyblk-2)*16<ydots)
1642 blocksize+=blocksize;
1643 return(blocksize);
1644 }
1645