File: common\fractals.c
1 /*
2 FRACTALS.C, FRACTALP.C and CALCFRAC.C actually calculate the fractal
3 images (well, SOMEBODY had to do it!). The modules are set up so that
4 all logic that is independent of any fractal-specific code is in
5 CALCFRAC.C, the code that IS fractal-specific is in FRACTALS.C, and the
6 structure that ties (we hope!) everything together is in FRACTALP.C.
7 Original author Tim Wegner, but just about ALL the authors have
8 contributed SOME code to this routine at one time or another, or
9 contributed to one of the many massive restructurings.
10
11 The Fractal-specific routines are divided into three categories:
12
13 1. Routines that are called once-per-orbit to calculate the orbit
14 value. These have names like "XxxxFractal", and their function
15 pointers are stored in fractalspecific[fractype].orbitcalc. EVERY
16 new fractal type needs one of these. Return 0 to continue iterations,
17 1 if we're done. Results for integer fractals are left in 'lnew.x' and
18 'lnew.y', for floating point fractals in 'new.x' and 'new.y'.
19
20 2. Routines that are called once per pixel to set various variables
21 prior to the orbit calculation. These have names like xxx_per_pixel
22 and are fairly generic - chances are one is right for your new type.
23 They are stored in fractalspecific[fractype].per_pixel.
24
25 3. Routines that are called once per screen to set various variables.
26 These have names like XxxxSetup, and are stored in
27 fractalspecific[fractype].per_image.
28
29 4. The main fractal routine. Usually this will be StandardFractal(),
30 but if you have written a stand-alone fractal routine independent
31 of the StandardFractal mechanisms, your routine name goes here,
32 stored in fractalspecific[fractype].calctype.per_image.
33
34 Adding a new fractal type should be simply a matter of adding an item
35 to the 'fractalspecific' structure, writing (or re-using one of the existing)
36 an appropriate setup, per_image, per_pixel, and orbit routines.
37
38 -------------------------------------------------------------------- */
39
40 #include <limits.h>
41 #include <string.h>
42 #ifdef __TURBOC__
44 #elif !defined(__386BSD__)
45 #include <malloc.h>
46 #endif
47 /* see Fractint.c for a description of the "include" hierarchy */
48 #include "port.h"
49 #include "prototyp.h"
50 #include "helpdefs.h"
51 #include "fractype.h"
52 #include "externs.h"
53
54
55 #define NEWTONDEGREELIMIT 100
56
57 _LCMPLX lcoefficient,lold,lnew,lparm, linit,ltmp,ltmp2,lparm2;
58 long ltempsqrx,ltempsqry;
59 int maxcolor;
60 int root, degree,basin;
61 double floatmin,floatmax;
62 double roverd, d1overd, threshold;
63 _CMPLX tmp2;
64 _CMPLX coefficient;
65 _CMPLX staticroots[16]; /* roots array for degree 16 or less */
66 _CMPLX *roots = staticroots;
67 struct MPC *MPCroots;
68 long FgHalf;
69 _CMPLX pwr;
70 int bitshiftless1; /* bit shift less 1 */
71
72 #ifndef sqr
73 #define sqr(x) ((x)*(x)) 74 #endif
75
76 #ifndef lsqr
78 #endif
79
80 #define modulus(z) (sqr((z).x)+sqr((z).y))
81 #define conjugate(pz) ((pz)->y = 0.0 - (pz)->y)
82 #define distance(z1,z2) (sqr((z1).x-(z2).x)+sqr((z1).y-(z2).y))
83 #define pMPsqr(z) (*pMPmul((z),(z)))
84 #define MPdistance(z1,z2) (*pMPadd(pMPsqr(*pMPsub((z1).x,(z2).x)),pMPsqr(*pMPsub((z1).y,(z2).y))))
85
86 double twopi = PI*2.0;
87 int c_exp;
88
89
90 /* These are local but I don't want to pass them as parameters */
91 _CMPLX parm,parm2;
92 _CMPLX *floatparm;
93 _LCMPLX *longparm; /* used here and in jb.c */
94
95 /* -------------------------------------------------------------------- */
96 /* These variables are external for speed's sake only */
97 /* -------------------------------------------------------------------- */
98
99 double sinx,cosx;
100 double siny,cosy;
101 double tmpexp;
102 double tempsqrx,tempsqry;
103
104 double foldxinitx,foldyinity,foldxinity,foldyinitx;
105 long oldxinitx,oldyinity,oldxinity,oldyinitx;
106 long longtmp;
107
108 /* These are for quaternions */
109 double qc,qci,qcj,qck;
110
111 /* temporary variables for trig use */
112 long lcosx, lsinx;
113 long lcosy, lsiny;
114
115 /*
116 ** details of finite attractors (required for Magnet Fractals)
117 ** (can also be used in "coloring in" the lakes of Julia types)
118 */
119
120 /*
121 ** pre-calculated values for fractal types Magnet2M & Magnet2J
122 */
123 _CMPLX T_Cm1; /* 3 * (floatparm - 1) */
124 _CMPLX T_Cm2; /* 3 * (floatparm - 2) */
125 _CMPLX T_Cm1Cm2; /* (floatparm - 1) * (floatparm - 2) */
126
127 void FloatPreCalcMagnet2(void) /* precalculation for Magnet2 (M & J) for speed */
128 {
129 T_Cm1.x = floatparm->x - 1.0; T_Cm1.y = floatparm->y;
130 T_Cm2.x = floatparm->x - 2.0; T_Cm2.y = floatparm->y;
131 T_Cm1Cm2.x = (T_Cm1.x * T_Cm2.x) - (T_Cm1.y * T_Cm2.y);
132 T_Cm1Cm2.y = (T_Cm1.x * T_Cm2.y) + (T_Cm1.y * T_Cm2.x);
133 T_Cm1.x += T_Cm1.x + T_Cm1.x; T_Cm1.y += T_Cm1.y + T_Cm1.y;
134 T_Cm2.x += T_Cm2.x + T_Cm2.x; T_Cm2.y += T_Cm2.y + T_Cm2.y;
135 }
136
137 /* -------------------------------------------------------------------- */
138 /* Bailout Routines Macros */
139 /* -------------------------------------------------------------------- */
140
141 int (near *floatbailout)(void);
142 int (near *longbailout)(void);
143 int (near *bignumbailout)(void);
144 int (near *bigfltbailout)(void);
145
146 #if 0
154 #endif
155 int near fpMODbailout(void)
156 {
157 tempsqrx=sqr(new.x);
158 tempsqry=sqr(new.y);
159 magnitude = tempsqrx + tempsqry;
160 if(magnitude >= rqlim) return(1);
161 old = new;
162 return(0);
163 }
164
165 int near fpREALbailout(void)
166 {
167 tempsqrx=sqr(new.x);
168 tempsqry=sqr(new.y);
169 magnitude = tempsqrx + tempsqry;
170 if(tempsqrx >= rqlim) return(1);
171 old = new;
172 return(0);
173 }
174
175 int near fpIMAGbailout(void)
176 {
177 tempsqrx=sqr(new.x);
178 tempsqry=sqr(new.y);
179 magnitude = tempsqrx + tempsqry;
180 if(tempsqry >= rqlim) return(1);
181 old = new;
182 return(0);
183 }
184
185 int near fpORbailout(void)
186 {
187 tempsqrx=sqr(new.x);
188 tempsqry=sqr(new.y);
189 magnitude = tempsqrx + tempsqry;
190 if(tempsqrx >= rqlim || tempsqry >= rqlim) return(1);
191 old = new;
192 return(0);
193 }
194
195 int near fpANDbailout(void)
196 {
197 tempsqrx=sqr(new.x);
198 tempsqry=sqr(new.y);
199 magnitude = tempsqrx + tempsqry;
200 if(tempsqrx >= rqlim && tempsqry >= rqlim) return(1);
201 old = new;
202 return(0);
203 }
204
205 int near fpMANHbailout(void)
206 {
207 double manhmag;
208 tempsqrx=sqr(new.x);
209 tempsqry=sqr(new.y);
210 magnitude = tempsqrx + tempsqry;
211 manhmag = fabs(new.x) + fabs(new.y);
212 if((manhmag * manhmag) >= rqlim) return(1);
213 old = new;
214 return(0);
215 }
216
217 int near fpMANRbailout(void)
218 {
219 double manrmag;
220 tempsqrx=sqr(new.x);
221 tempsqry=sqr(new.y);
222 magnitude = tempsqrx + tempsqry;
223 manrmag = new.x + new.y; /* don't need abs() since we square it next */
224 if((manrmag * manrmag) >= rqlim) return(1);
225 old = new;
226 return(0);
227 }
228
229 #define FLOATTRIGBAILOUT() \
230 if (fabs(old.y) >= rqlim2) return(1);
231
232 #define LONGTRIGBAILOUT() \
233 if(labs(lold.y) >= llimit2) { return(1);}
234
235 #define LONGXYTRIGBAILOUT() \
236 if(labs(lold.x) >= llimit2 || labs(lold.y) >= llimit2)\
237 { return(1);}
238
239 #define FLOATXYTRIGBAILOUT() \
240 if (fabs(old.x) >= rqlim2 || fabs(old.y) >= rqlim2) return(1);
241
242 #define FLOATHTRIGBAILOUT() \
243 if (fabs(old.x) >= rqlim2) return(1);
244
245 #define LONGHTRIGBAILOUT() \
246 if(labs(lold.x) >= llimit2) { return(1);}
247
248 #define TRIG16CHECK(X) \
249 if(labs((X)) > l16triglim) { return(1);}
250
251 #define OLD_FLOATEXPBAILOUT() \
252 if (fabs(old.y) >= 1.0e8) return(1);\
253 if (fabs(old.x) >= 6.4e2) return(1);
254
255 #define FLOATEXPBAILOUT() \
256 if (fabs(old.y) >= 1.0e3) return(1);\
257 if (fabs(old.x) >= 8) return(1);
258
259 #define LONGEXPBAILOUT() \
260 if (labs(lold.y) >= (1000L<<bitshift)) return(1);\
261 if (labs(lold.x) >= (8L<<bitshift)) return(1);
262
263 #if 0
272 #endif
273
274 #define LTRIGARG(X) \
275 if(labs((X)) > l16triglim)\
276 {\
277 double tmp;\
278 tmp = (X);\
279 tmp /= fudge;\
280 tmp = fmod(tmp,twopi);\
281 tmp *= fudge;\
282 (X) = (long)tmp;\
283 }\
284
285 static int near Halleybailout(void)
286 {
287 if ( fabs(modulus(new)-modulus(old)) < parm2.x)
288 return(1);
289 old = new;
290 return(0);
291 }
292
293 #ifndef XFRACT
294 #define MPCmod(m) (*pMPadd(*pMPmul((m).x, (m).x), *pMPmul((m).y, (m).y)))
295 struct MPC mpcold, mpcnew, mpctmp, mpctmp1;
296 struct MP mptmpparm2x;
297
298 #if (_MSC_VER >= 700)
299 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
300 #endif
301
302 static int near MPCHalleybailout(void)
303 {
304 static struct MP mptmpbailout;
305 mptmpbailout = *MPabs(*pMPsub(MPCmod(mpcnew), MPCmod(mpcold)));
306 if (pMPcmp(mptmpbailout, mptmpparm2x) < 0)
307 return(1);
308 mpcold = mpcnew;
309 return(0);
310 }
311 #if (_MSC_VER >= 700)
312 #pragma code_seg () /* back to normal segment */
313 #endif
314 #endif
315
316 #ifdef XFRACT
338 #endif
339
340 /* -------------------------------------------------------------------- */
341 /* Fractal (once per iteration) routines */
342 /* -------------------------------------------------------------------- */
343 static double xt, yt, t2;
344
345 /* Raise complex number (base) to the (exp) power, storing the result
346 ** in complex (result).
347 */
348 void cpower(_CMPLX *base, int exp, _CMPLX *result)
349 {
350 if (exp<0) {
351 cpower(base,-exp,result);
352 CMPLXrecip(*result,*result);
353 return;
354 }
355
356 xt = base->x; yt = base->y;
357
358 if (exp & 1)
359 {
360 result->x = xt;
361 result->y = yt;
362 }
363 else
364 {
365 result->x = 1.0;
366 result->y = 0.0;
367 }
368
369 exp >>= 1;
370 while (exp)
371 {
372 t2 = xt * xt - yt * yt;
373 yt = 2 * xt * yt;
374 xt = t2;
375
376 if (exp & 1)
377 {
378 t2 = xt * result->x - yt * result->y;
379 result->y = result->y * xt + yt * result->x;
380 result->x = t2;
381 }
382 exp >>= 1;
383 }
384 }
385
386 #ifndef XFRACT
387 /* long version */
388 static long lxt, lyt, lt2;
389 int
390 lcpower(_LCMPLX *base, int exp, _LCMPLX *result, int bitshift)
391 {
392 static long maxarg;
393 maxarg = 64L<<bitshift;
394
395 if (exp<0) {
396 overflow = lcpower(base,-exp,result,bitshift);
397 LCMPLXrecip(*result,*result);
398 return(overflow);
399 }
400
401 overflow = 0;
402 lxt = base->x; lyt = base->y;
403
404 if (exp & 1)
405 {
406 result->x = lxt;
407 result->y = lyt;
408 }
409 else
410 {
411 result->x = 1L<<bitshift;
412 result->y = 0L;
413 }
414
415 exp >>= 1;
416 while (exp)
417 {
418 /*
419 if(labs(lxt) >= maxarg || labs(lyt) >= maxarg)
420 return(-1);
421 */
422 lt2 = multiply(lxt, lxt, bitshift) - multiply(lyt,lyt,bitshift);
423 lyt = multiply(lxt,lyt,bitshiftless1);
424 if(overflow)
425 return(overflow);
426 lxt = lt2;
427
428 if (exp & 1)
429 {
430 lt2 = multiply(lxt,result->x, bitshift) - multiply(lyt,result->y,bitshift);
431 result->y = multiply(result->y,lxt,bitshift) + multiply(lyt,result->x,bitshift);
432 result->x = lt2;
433 }
434 exp >>= 1;
435 }
436 if(result->x == 0 && result->y == 0)
437 overflow = 1;
438 return(overflow);
439 }
440 #if 0
471 #endif
472 #endif
473
474 #ifdef XFRACT /* fractint uses the NewtonFractal2 code in newton.asm */
475
476 int complex_div(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
477 int complex_mult(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
478
479 /* Distance of complex z from unit circle */
480 #define DIST1(z) (((z).x-1.0)*((z).x-1.0)+((z).y)*((z).y))
481 #define LDIST1(z) (lsqr((((z).x)-fudge)) + lsqr(((z).y)))
482
483
484 int NewtonFractal2(void)
485 {
486 static char start=1;
487 if(start)
488 {
489 start = 0;
490 }
491 cpower(&old, degree-1, &tmp);
492 complex_mult(tmp, old, &new);
493
494 if (DIST1(new) < threshold)
495 {
496 if(fractype==NEWTBASIN || fractype==MPNEWTBASIN)
497 {
498 long tmpcolor;
499 int i;
500 tmpcolor = -1;
501 /* this code determines which degree-th root of root the
502 Newton formula converges to. The roots of a 1 are
503 distributed on a circle of radius 1 about the origin. */
504 for(i=0;i<degree;i++)
505 /* color in alternating shades with iteration according to
506 which root of 1 it converged to */
507 if(distance(roots[i],old) < threshold)
508 {
509 if (basin==2) {
510 tmpcolor = 1+(i&7)+((coloriter&1)<<3);
511 } else {
512 tmpcolor = 1+i;
513 }
514 break;
515 }
516 if(tmpcolor == -1)
517 coloriter = maxcolor;
518 else
519 coloriter = tmpcolor;
520 }
521 return(1);
522 }
523 new.x = d1overd * new.x + roverd;
524 new.y *= d1overd;
525
526 /* Watch for divide underflow */
527 if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN)
528 return(1);
529 else
530 {
531 t2 = 1.0 / t2;
532 old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
533 old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
534 }
535 return(0);
536 }
537
538 int
539 complex_mult(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz)
540 {
541 pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
542 pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
543 return(0);
544 }
545
546 int
547 complex_div(_CMPLX numerator,_CMPLX denominator,_CMPLX *pout)
548 {
549 double mod;
550 if((mod = modulus(denominator)) < FLT_MIN)
551 return(1);
552 conjugate(&denominator);
553 complex_mult(numerator,denominator,pout);
554 pout->x = pout->x/mod;
555 pout->y = pout->y/mod;
556 return(0);
557 } 558 #endif /* newton code only used by xfractint */
559
560 #ifndef XFRACT
561 struct MP mproverd, mpd1overd, mpthreshold;
562 struct MP mpt2;
563 struct MP mpone;
564 #endif
565
566 #if (_MSC_VER >= 700)
567 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
568 #endif
569 int MPCNewtonFractal(void)
570 {
571 #ifndef XFRACT
572 MPOverflow = 0;
573 mpctmp = MPCpow(mpcold,degree-1);
574
575 mpcnew.x = *pMPsub(*pMPmul(mpctmp.x,mpcold.x),*pMPmul(mpctmp.y,mpcold.y));
576 mpcnew.y = *pMPadd(*pMPmul(mpctmp.x,mpcold.y),*pMPmul(mpctmp.y,mpcold.x));
577 mpctmp1.x = *pMPsub(mpcnew.x, MPCone.x);
578 mpctmp1.y = *pMPsub(mpcnew.y, MPCone.y);
579 if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
580 {
581 if(fractype==MPNEWTBASIN)
582 {
583 long tmpcolor;
584 int i;
585 tmpcolor = -1;
586 for(i=0;i<degree;i++)
587 if(pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
588 {
589 if(basin==2)
590 tmpcolor = 1+(i&7) + ((coloriter&1)<<3);
591 else
592 tmpcolor = 1+i;
593 break;
594 }
595 if(tmpcolor == -1)
596 coloriter = maxcolor;
597 else
598 coloriter = tmpcolor;
599 }
600 return(1);
601 }
602
603 mpcnew.x = *pMPadd(*pMPmul(mpd1overd,mpcnew.x),mproverd);
604 mpcnew.y = *pMPmul(mpcnew.y,mpd1overd);
605 mpt2 = MPCmod(mpctmp);
606 mpt2 = *pMPdiv(mpone,mpt2);
607 mpcold.x = *pMPmul(mpt2,(*pMPadd(*pMPmul(mpcnew.x,mpctmp.x),*pMPmul(mpcnew.y,mpctmp.y))));
608 mpcold.y = *pMPmul(mpt2,(*pMPsub(*pMPmul(mpcnew.y,mpctmp.x),*pMPmul(mpcnew.x,mpctmp.y))));
609 new.x = *pMP2d(mpcold.x);
610 new.y = *pMP2d(mpcold.y);
611 return(MPOverflow);
614 #endif
615 }
616 #if (_MSC_VER >= 700)
617 #pragma code_seg () /* back to normal segment */
618 #endif
619
620 int
621 Barnsley1Fractal(void)
622 {
623 #ifndef XFRACT
624 /* Barnsley's Mandelbrot type M1 from "Fractals
625 Everywhere" by Michael Barnsley, p. 322 */
626
627 /* calculate intermediate products */
628 oldxinitx = multiply(lold.x, longparm->x, bitshift);
629 oldyinity = multiply(lold.y, longparm->y, bitshift);
630 oldxinity = multiply(lold.x, longparm->y, bitshift);
631 oldyinitx = multiply(lold.y, longparm->x, bitshift);
632 /* orbit calculation */
633 if(lold.x >= 0)
634 {
635 lnew.x = (oldxinitx - longparm->x - oldyinity);
636 lnew.y = (oldyinitx - longparm->y + oldxinity);
637 }
638 else
639 {
640 lnew.x = (oldxinitx + longparm->x - oldyinity);
641 lnew.y = (oldyinitx + longparm->y + oldxinity);
642 }
643 return(longbailout());
646 #endif
647 }
648
649 int
650 Barnsley1FPFractal(void)
651 {
652 /* Barnsley's Mandelbrot type M1 from "Fractals
653 Everywhere" by Michael Barnsley, p. 322 */
654 /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
655
656 /* calculate intermediate products */
657 foldxinitx = old.x * floatparm->x;
658 foldyinity = old.y * floatparm->y;
659 foldxinity = old.x * floatparm->y;
660 foldyinitx = old.y * floatparm->x;
661 /* orbit calculation */
662 if(old.x >= 0)
663 {
664 new.x = (foldxinitx - floatparm->x - foldyinity);
665 new.y = (foldyinitx - floatparm->y + foldxinity);
666 }
667 else
668 {
669 new.x = (foldxinitx + floatparm->x - foldyinity);
670 new.y = (foldyinitx + floatparm->y + foldxinity);
671 }
672 return(floatbailout());
673 }
674
675 int
676 Barnsley2Fractal(void)
677 {
678 #ifndef XFRACT
679 /* An unnamed Mandelbrot/Julia function from "Fractals
680 Everywhere" by Michael Barnsley, p. 331, example 4.2 */
681 /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
682
683 /* calculate intermediate products */
684 oldxinitx = multiply(lold.x, longparm->x, bitshift);
685 oldyinity = multiply(lold.y, longparm->y, bitshift);
686 oldxinity = multiply(lold.x, longparm->y, bitshift);
687 oldyinitx = multiply(lold.y, longparm->x, bitshift);
688
689 /* orbit calculation */
690 if(oldxinity + oldyinitx >= 0)
691 {
692 lnew.x = oldxinitx - longparm->x - oldyinity;
693 lnew.y = oldyinitx - longparm->y + oldxinity;
694 }
695 else
696 {
697 lnew.x = oldxinitx + longparm->x - oldyinity;
698 lnew.y = oldyinitx + longparm->y + oldxinity;
699 }
700 return(longbailout());
703 #endif
704 }
705
706 int
707 Barnsley2FPFractal(void)
708 {
709 /* An unnamed Mandelbrot/Julia function from "Fractals
710 Everywhere" by Michael Barnsley, p. 331, example 4.2 */
711
712 /* calculate intermediate products */
713 foldxinitx = old.x * floatparm->x;
714 foldyinity = old.y * floatparm->y;
715 foldxinity = old.x * floatparm->y;
716 foldyinitx = old.y * floatparm->x;
717
718 /* orbit calculation */
719 if(foldxinity + foldyinitx >= 0)
720 {
721 new.x = foldxinitx - floatparm->x - foldyinity;
722 new.y = foldyinitx - floatparm->y + foldxinity;
723 }
724 else
725 {
726 new.x = foldxinitx + floatparm->x - foldyinity;
727 new.y = foldyinitx + floatparm->y + foldxinity;
728 }
729 return(floatbailout());
730 }
731
732 int
733 JuliaFractal(void)
734 {
735 #ifndef XFRACT
736 /* used for C prototype of fast integer math routines for classic
737 Mandelbrot and Julia */
738 lnew.x = ltempsqrx - ltempsqry + longparm->x;
739 lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
740 return(longbailout());
744 #endif
745 }
746
747 int
748 JuliafpFractal(void)
749 {
750 /* floating point version of classical Mandelbrot/Julia */
751 /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
752 new.x = tempsqrx - tempsqry + floatparm->x;
753 new.y = 2.0 * old.x * old.y + floatparm->y;
754 return(floatbailout());
755 }
756
757 int
758 LambdaFPFractal(void)
759 {
760 /* variation of classical Mandelbrot/Julia */
761 /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
762
763 tempsqrx = old.x - tempsqrx + tempsqry;
764 tempsqry = -(old.y * old.x);
765 tempsqry += tempsqry + old.y;
766
767 new.x = floatparm->x * tempsqrx - floatparm->y * tempsqry;
768 new.y = floatparm->x * tempsqry + floatparm->y * tempsqrx;
769 return(floatbailout());
770 }
771
772 int
773 LambdaFractal(void)
774 {
775 #ifndef XFRACT
776 /* variation of classical Mandelbrot/Julia */
777
778 /* in complex math) temp = Z * (1-Z) */
779 ltempsqrx = lold.x - ltempsqrx + ltempsqry;
780 ltempsqry = lold.y
781 - multiply(lold.y, lold.x, bitshiftless1);
782 /* (in complex math) Z = Lambda * Z */
783 lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
784 - multiply(longparm->y, ltempsqry, bitshift);
785 lnew.y = multiply(longparm->x, ltempsqry, bitshift)
786 + multiply(longparm->y, ltempsqrx, bitshift);
787 return(longbailout());
790 #endif
791 }
792
793 int
794 SierpinskiFractal(void)
795 {
796 #ifndef XFRACT
797 /* following code translated from basic - see "Fractals
798 Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
799 lnew.x = (lold.x << 1); /* new.x = 2 * old.x */
800 lnew.y = (lold.y << 1); /* new.y = 2 * old.y */
801 if(lold.y > ltmp.y) /* if old.y > .5 */
802 lnew.y = lnew.y - ltmp.x; /* new.y = 2 * old.y - 1 */
803 else if(lold.x > ltmp.y) /* if old.x > .5 */
804 lnew.x = lnew.x - ltmp.x; /* new.x = 2 * old.x - 1 */
805 /* end barnsley code */
806 return(longbailout());
809 #endif
810 }
811
812 int
813 SierpinskiFPFractal(void)
814 {
815 /* following code translated from basic - see "Fractals
816 Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
817
818 new.x = old.x + old.x;
819 new.y = old.y + old.y;
820 if(old.y > .5)
821 new.y = new.y - 1;
822 else if (old.x > .5)
823 new.x = new.x - 1;
824
825 /* end barnsley code */
826 return(floatbailout());
827 }
828
829 int
830 LambdaexponentFractal(void)
831 {
832 /* found this in "Science of Fractal Images" */
833 if (save_release > 2002) { /* need braces since these are macros */
834 FLOATEXPBAILOUT();
835 }
836 else {
837 OLD_FLOATEXPBAILOUT();
838 }
839 FPUsincos (&old.y,&siny,&cosy);
840
841 if (old.x >= rqlim && cosy >= 0.0) return(1);
842 tmpexp = exp(old.x);
843 tmp.x = tmpexp*cosy;
844 tmp.y = tmpexp*siny;
845
846 /*multiply by lamda */
847 new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
848 new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
849 old = new;
850 return(0);
851 }
852
853 int
854 LongLambdaexponentFractal(void)
855 {
856 #ifndef XFRACT
857 /* found this in "Science of Fractal Images" */
858 LONGEXPBAILOUT();
859
860 SinCos086 (lold.y, &lsiny, &lcosy);
861
862 if (lold.x >= llimit && lcosy >= 0L) return(1);
863 longtmp = Exp086(lold.x);
864
865 ltmp.x = multiply(longtmp, lcosy, bitshift);
866 ltmp.y = multiply(longtmp, lsiny, bitshift);
867
868 lnew.x = multiply(longparm->x, ltmp.x, bitshift)
869 - multiply(longparm->y, ltmp.y, bitshift);
870 lnew.y = multiply(longparm->x, ltmp.y, bitshift)
871 + multiply(longparm->y, ltmp.x, bitshift);
872 lold = lnew;
873 return(0);
876 #endif
877 }
878
879 int
880 FloatTrigPlusExponentFractal(void)
881 {
882 /* another Scientific American biomorph type */
883 /* z(n+1) = e**z(n) + trig(z(n)) + C */
884
885 if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
886 tmpexp = exp(old.x);
887 FPUsincos (&old.y,&siny,&cosy);
888 CMPLXtrig0(old,new);
889
890 /*new = trig(old) + e**old + C */
891 new.x += tmpexp*cosy + floatparm->x;
892 new.y += tmpexp*siny + floatparm->y;
893 return(floatbailout());
894 }
895
896 int
897 LongTrigPlusExponentFractal(void)
898 {
899 #ifndef XFRACT
900 /* calculate exp(z) */
901
902 /* domain check for fast transcendental functions */
903 TRIG16CHECK(lold.x);
904 TRIG16CHECK(lold.y);
905
906 longtmp = Exp086(lold.x);
907 SinCos086 (lold.y, &lsiny, &lcosy);
908 LCMPLXtrig0(lold,lnew);
909 lnew.x += multiply(longtmp, lcosy, bitshift) + longparm->x;
910 lnew.y += multiply(longtmp, lsiny, bitshift) + longparm->y;
911 return(longbailout());
914 #endif
915 }
916
917 int
918 MarksLambdaFractal(void)
919 {
920 /* Mark Peterson's variation of "lambda" function */
921
922 /* Z1 = (C^(exp-1) * Z**2) + C */
923 #ifndef XFRACT
924 ltmp.x = ltempsqrx - ltempsqry;
925 ltmp.y = multiply(lold.x ,lold.y ,bitshiftless1);
926
927 lnew.x = multiply(lcoefficient.x, ltmp.x, bitshift)
928 - multiply(lcoefficient.y, ltmp.y, bitshift) + longparm->x;
929 lnew.y = multiply(lcoefficient.x, ltmp.y, bitshift)
930 + multiply(lcoefficient.y, ltmp.x, bitshift) + longparm->y;
931
932 return(longbailout());
935 #endif
936 }
937
938 int
939 MarksLambdafpFractal(void)
940 {
941 /* Mark Peterson's variation of "lambda" function */
942
943 /* Z1 = (C^(exp-1) * Z**2) + C */
944 tmp.x = tempsqrx - tempsqry;
945 tmp.y = old.x * old.y *2;
946
947 new.x = coefficient.x * tmp.x - coefficient.y * tmp.y + floatparm->x;
948 new.y = coefficient.x * tmp.y + coefficient.y * tmp.x + floatparm->y;
949
950 return(floatbailout());
951 }
952
953
954 long XXOne, FgOne, FgTwo;
955
956 int
957 UnityFractal(void)
958 {
959 #ifndef XFRACT
960 /* brought to you by Mark Peterson - you won't find this in any fractal
961 books unless they saw it here first - Mark invented it! */
962 XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
963 if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delmin))
964 return(1);
965 lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
966 lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
967 lnew=lold; /* TW added this line */
968 return(0);
971 #endif
972 }
973
974 int
975 UnityfpFractal(void)
976 {
977 double XXOne;
978 /* brought to you by Mark Peterson - you won't find this in any fractal
979 books unless they saw it here first - Mark invented it! */
980
981 XXOne = sqr(old.x) + sqr(old.y);
982 if((XXOne > 2.0) || (fabs(XXOne - 1.0) < ddelmin))
983 return(1);
984 old.y = (2.0 - XXOne)* old.x;
985 old.x = (2.0 - XXOne)* old.y;
986 new=old; /* TW added this line */
987 return(0);
988 }
989
990 int
991 Mandel4Fractal(void)
992 {
993 /* By writing this code, Bert has left behind the excuse "don't
994 know what a fractal is, just know how to make'em go fast".
995 Bert is hereby declared a bonafide fractal expert! Supposedly
996 this routine calculates the Mandelbrot/Julia set based on the
997 polynomial z**4 + lambda, but I wouldn't know -- can't follow
998 all that integer math speedup stuff - Tim */
999
1000 /* first, compute (x + iy)**2 */
1001 #ifndef XFRACT
1002 lnew.x = ltempsqrx - ltempsqry;
1003 lnew.y = multiply(lold.x, lold.y, bitshiftless1);
1004 if (longbailout()) return(1);
1005
1006 /* then, compute ((x + iy)**2)**2 + lambda */
1007 lnew.x = ltempsqrx - ltempsqry + longparm->x;
1008 lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
1009 return(longbailout());
1012 #endif
1013 }
1014
1015 int
1016 Mandel4fpFractal(void)
1017 {
1018 /* first, compute (x + iy)**2 */
1019 new.x = tempsqrx - tempsqry;
1020 new.y = old.x*old.y*2;
1021 if (floatbailout()) return(1);
1022
1023 /* then, compute ((x + iy)**2)**2 + lambda */
1024 new.x = tempsqrx - tempsqry + floatparm->x;
1025 new.y = old.x*old.y*2 + floatparm->y;
1026 return(floatbailout());
1027 }
1028
1029 int
1030 floatZtozPluszpwrFractal(void)
1031 {
1032 cpower(&old,(int)param[2],&new);
1033 old = ComplexPower(old,old);
1034 new.x = new.x + old.x +floatparm->x;
1035 new.y = new.y + old.y +floatparm->y;
1036 return(floatbailout());
1037 }
1038
1039 int
1040 longZpowerFractal(void)
1041 {
1042 #ifndef XFRACT
1043 if(lcpower(&lold,c_exp,&lnew,bitshift))
1044 lnew.x = lnew.y = 8L<<bitshift;
1045 lnew.x += longparm->x;
1046 lnew.y += longparm->y;
1047 return(longbailout());
1050 #endif
1051 }
1052
1053 int
1054 longCmplxZpowerFractal(void)
1055 {
1056 #ifndef XFRACT
1057 _CMPLX x, y;
1058
1059 x.x = (double)lold.x / fudge;
1060 x.y = (double)lold.y / fudge;
1061 y.x = (double)lparm2.x / fudge;
1062 y.y = (double)lparm2.y / fudge;
1063 x = ComplexPower(x, y);
1064 if(fabs(x.x) < fgLimit && fabs(x.y) < fgLimit) {
1065 lnew.x = (long)(x.x * fudge);
1066 lnew.y = (long)(x.y * fudge);
1067 }
1068 else
1069 overflow = 1;
1070 lnew.x += longparm->x;
1071 lnew.y += longparm->y;
1072 return(longbailout());
1075 #endif
1076 }
1077
1078 int
1079 floatZpowerFractal(void)
1080 {
1081 cpower(&old,c_exp,&new);
1082 new.x += floatparm->x;
1083 new.y += floatparm->y;
1084 return(floatbailout());
1085 }
1086
1087 int
1088 floatCmplxZpowerFractal(void)
1089 {
1090 new = ComplexPower(old, parm2);
1091 new.x += floatparm->x;
1092 new.y += floatparm->y;
1093 return(floatbailout());
1094 }
1095
1096 int
1097 Barnsley3Fractal(void)
1098 {
1099 /* An unnamed Mandelbrot/Julia function from "Fractals
1100 Everywhere" by Michael Barnsley, p. 292, example 4.1 */
1101
1102 /* calculate intermediate products */
1103 #ifndef XFRACT
1104 oldxinitx = multiply(lold.x, lold.x, bitshift);
1105 oldyinity = multiply(lold.y, lold.y, bitshift);
1106 oldxinity = multiply(lold.x, lold.y, bitshift);
1107
1108 /* orbit calculation */
1109 if(lold.x > 0)
1110 {
1111 lnew.x = oldxinitx - oldyinity - fudge;
1112 lnew.y = oldxinity << 1;
1113 }
1114 else
1115 {
1116 lnew.x = oldxinitx - oldyinity - fudge
1117 + multiply(longparm->x,lold.x,bitshift);
1118 lnew.y = oldxinity <<1;
1119
1120 /* This term added by Tim Wegner to make dependent on the
1121 imaginary part of the parameter. (Otherwise Mandelbrot
1122 is uninteresting. */
1123 lnew.y += multiply(longparm->y,lold.x,bitshift);
1124 }
1125 return(longbailout());
1128 #endif
1129 }
1130
1131 int
1132 Barnsley3FPFractal(void)
1133 {
1134 /* An unnamed Mandelbrot/Julia function from "Fractals
1135 Everywhere" by Michael Barnsley, p. 292, example 4.1 */
1136
1137
1138 /* calculate intermediate products */
1139 foldxinitx = old.x * old.x;
1140 foldyinity = old.y * old.y;
1141 foldxinity = old.x * old.y;
1142
1143 /* orbit calculation */
1144 if(old.x > 0)
1145 {
1146 new.x = foldxinitx - foldyinity - 1.0;
1147 new.y = foldxinity * 2;
1148 }
1149 else
1150 {
1151 new.x = foldxinitx - foldyinity -1.0 + floatparm->x * old.x;
1152 new.y = foldxinity * 2;
1153
1154 /* This term added by Tim Wegner to make dependent on the
1155 imaginary part of the parameter. (Otherwise Mandelbrot
1156 is uninteresting. */
1157 new.y += floatparm->y * old.x;
1158 }
1159 return(floatbailout());
1160 }
1161
1162 int
1163 TrigPlusZsquaredFractal(void)
1164 {
1165 #ifndef XFRACT
1166 /* From Scientific American, July 1989 */
1167 /* A Biomorph */
1168 /* z(n+1) = trig(z(n))+z(n)**2+C */
1169 LCMPLXtrig0(lold,lnew);
1170 lnew.x += ltempsqrx - ltempsqry + longparm->x;
1171 lnew.y += multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
1172 return(longbailout());
1175 #endif
1176 }
1177
1178 int
1179 TrigPlusZsquaredfpFractal(void)
1180 {
1181 /* From Scientific American, July 1989 */
1182 /* A Biomorph */
1183 /* z(n+1) = trig(z(n))+z(n)**2+C */
1184
1185 CMPLXtrig0(old,new);
1186 new.x += tempsqrx - tempsqry + floatparm->x;
1187 new.y += 2.0 * old.x * old.y + floatparm->y;
1188 return(floatbailout());
1189 }
1190
1191 int
1192 Richard8fpFractal(void)
1193 {
1194 /* Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
1195 CMPLXtrig0(old,new);
1196 /* CMPLXtrig1(*floatparm,tmp); */
1197 new.x += tmp.x;
1198 new.y += tmp.y;
1199 return(floatbailout());
1200 }
1201
1202 int
1203 Richard8Fractal(void)
1204 {
1205 #ifndef XFRACT
1206 /* Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
1207 LCMPLXtrig0(lold,lnew);
1208 /* LCMPLXtrig1(*longparm,ltmp); */
1209 lnew.x += ltmp.x;
1210 lnew.y += ltmp.y;
1211 return(longbailout());
1214 #endif
1215 }
1216
1217 int
1218 PopcornFractal_Old(void)
1219 {
1220 tmp = old;
1221 tmp.x *= 3.0;
1222 tmp.y *= 3.0;
1223 FPUsincos(&tmp.x,&sinx,&cosx);
1224 FPUsincos(&tmp.y,&siny,&cosy);
1225 tmp.x = sinx/cosx + old.x;
1226 tmp.y = siny/cosy + old.y;
1227 FPUsincos(&tmp.x,&sinx,&cosx);
1228 FPUsincos(&tmp.y,&siny,&cosy);
1229 new.x = old.x - parm.x*siny;
1230 new.y = old.y - parm.x*sinx;
1231 /*
1232 new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
1233 new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
1234 */
1235 if(plot == noplot)
1236 {
1237 plot_orbit(new.x,new.y,1+row%colors);
1238 old = new;
1239 }
1240 else
1241 /* FLOATBAILOUT(); */
1242 /* PB The above line was weird, not what it seems to be! But, bracketing
1243 it or always doing it (either of which seem more likely to be what
1244 was intended) changes the image for the worse, so I'm not touching it.
1245 Same applies to int form in next routine. */
1246 /* PB later: recoded inline, still leaving it weird */
1247 tempsqrx = sqr(new.x);
1248 tempsqry = sqr(new.y);
1249 if((magnitude = tempsqrx + tempsqry) >= rqlim) return(1);
1250 old = new;
1251 return(0);
1252 }
1253
1254 int
1255 PopcornFractal(void)
1256 {
1257 tmp = old;
1258 tmp.x *= 3.0;
1259 tmp.y *= 3.0;
1260 FPUsincos(&tmp.x,&sinx,&cosx);
1261 FPUsincos(&tmp.y,&siny,&cosy);
1262 tmp.x = sinx/cosx + old.x;
1263 tmp.y = siny/cosy + old.y;
1264 FPUsincos(&tmp.x,&sinx,&cosx);
1265 FPUsincos(&tmp.y,&siny,&cosy);
1266 new.x = old.x - parm.x*siny;
1267 new.y = old.y - parm.x*sinx;
1268 /*
1269 new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
1270 new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
1271 */
1272 if(plot == noplot)
1273 {
1274 plot_orbit(new.x,new.y,1+row%colors);
1275 old = new;
1276 }
1277 /* else */
1278 /* FLOATBAILOUT(); */
1279 /* PB The above line was weird, not what it seems to be! But, bracketing
1280 it or always doing it (either of which seem more likely to be what
1281 was intended) changes the image for the worse, so I'm not touching it.
1282 Same applies to int form in next routine. */
1283 /* PB later: recoded inline, still leaving it weird */
1284 /* JCO: sqr's should always be done, else magnitude could be wrong */
1285 tempsqrx = sqr(new.x);
1286 tempsqry = sqr(new.y);
1287 if((magnitude = tempsqrx + tempsqry) >= rqlim
1288 || fabs(new.x) > rqlim2 || fabs(new.y) > rqlim2 )
1289 return(1);
1290 old = new;
1291 return(0);
1292 }
1293
1294 int
1295 LPopcornFractal_Old(void)
1296 {
1297 #ifndef XFRACT
1298 ltmp = lold;
1299 ltmp.x *= 3L;
1300 ltmp.y *= 3L;
1301 LTRIGARG(ltmp.x);
1302 LTRIGARG(ltmp.y);
1303 SinCos086(ltmp.x,&lsinx,&lcosx);
1304 SinCos086(ltmp.y,&lsiny,&lcosy);
1305 ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
1306 ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
1307 LTRIGARG(ltmp.x);
1308 LTRIGARG(ltmp.y);
1309 SinCos086(ltmp.x,&lsinx,&lcosx);
1310 SinCos086(ltmp.y,&lsiny,&lcosy);
1311 lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
1312 lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
1313 if(plot == noplot)
1314 {
1315 iplot_orbit(lnew.x,lnew.y,1+row%colors);
1316 lold = lnew;
1317 }
1318 else
1319 /* LONGBAILOUT(); */
1320 /* PB above still the old way, is weird, see notes in FP popcorn case */
1321 {
1322 ltempsqrx = lsqr(lnew.x);
1323 ltempsqry = lsqr(lnew.y);
1324 }
1325 lmagnitud = ltempsqrx + ltempsqry;
1326 if (lmagnitud >= llimit || lmagnitud < 0 || labs(lnew.x) > llimit2
1327 || labs(lnew.y) > llimit2)
1328 return(1);
1329 lold = lnew;
1330 return(0);
1333 #endif
1334 }
1335
1336 int
1337 LPopcornFractal(void)
1338 {
1339 #ifndef XFRACT
1340 ltmp = lold;
1341 ltmp.x *= 3L;
1342 ltmp.y *= 3L;
1343 LTRIGARG(ltmp.x);
1344 LTRIGARG(ltmp.y);
1345 SinCos086(ltmp.x,&lsinx,&lcosx);
1346 SinCos086(ltmp.y,&lsiny,&lcosy);
1347 ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
1348 ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
1349 LTRIGARG(ltmp.x);
1350 LTRIGARG(ltmp.y);
1351 SinCos086(ltmp.x,&lsinx,&lcosx);
1352 SinCos086(ltmp.y,&lsiny,&lcosy);
1353 lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
1354 lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
1355 if(plot == noplot)
1356 {
1357 iplot_orbit(lnew.x,lnew.y,1+row%colors);
1358 lold = lnew;
1359 }
1360 /* else */
1361 /* JCO: sqr's should always be done, else magnitude could be wrong */
1362 ltempsqrx = lsqr(lnew.x);
1363 ltempsqry = lsqr(lnew.y);
1364 lmagnitud = ltempsqrx + ltempsqry;
1365 if (lmagnitud >= llimit || lmagnitud < 0
1366 || labs(lnew.x) > llimit2
1367 || labs(lnew.y) > llimit2)
1368 return(1);
1369 lold = lnew;
1370 return(0);
1373 #endif
1374 }
1375
1376 /* Popcorn generalization proposed by HB */
1377
1378 int
1379 PopcornFractalFn(void)
1380 {
1381 _CMPLX tmpx;
1382 _CMPLX tmpy;
1383
1384 /* tmpx contains the generalized value of the old real "x" equation */
1385 CMPLXtimesreal(parm2,old.y,tmp); /* tmp = (C * old.y) */
1386 CMPLXtrig1(tmp,tmpx); /* tmpx = trig1(tmp) */
1387 tmpx.x += old.y; /* tmpx = old.y + trig1(tmp) */
1388 CMPLXtrig0(tmpx,tmp); /* tmp = trig0(tmpx) */
1389 CMPLXmult(tmp,parm,tmpx); /* tmpx = tmp * h */
1390
1391 /* tmpy contains the generalized value of the old real "y" equation */
1392 CMPLXtimesreal(parm2,old.x,tmp); /* tmp = (C * old.x) */
1393 CMPLXtrig3(tmp,tmpy); /* tmpy = trig3(tmp) */
1394 tmpy.x += old.x; /* tmpy = old.x + trig1(tmp) */
1395 CMPLXtrig2(tmpy,tmp); /* tmp = trig2(tmpy) */
1396
1397 CMPLXmult(tmp,parm,tmpy); /* tmpy = tmp * h */
1398
1399 new.x = old.x - tmpx.x - tmpy.y;
1400 new.y = old.y - tmpy.x - tmpx.y;
1401
1402 if(plot == noplot)
1403 {
1404 plot_orbit(new.x,new.y,1+row%colors);
1405 old = new;
1406 }
1407
1408 tempsqrx = sqr(new.x);
1409 tempsqry = sqr(new.y);
1410 if((magnitude = tempsqrx + tempsqry) >= rqlim
1411 || fabs(new.x) > rqlim2 || fabs(new.y) > rqlim2 )
1412 return(1);
1413 old = new;
1414 return(0);
1415 }
1416
1417 #define FIX_OVERFLOW(arg) if(overflow) \
1418 { \
1419 (arg).x = fudge;\
1420 (arg).y = 0;\
1421 overflow = 0;\
1422 }
1423
1424 int
1425 LPopcornFractalFn(void)
1426 {
1427 #ifndef XFRACT
1428 _LCMPLX ltmpx, ltmpy;
1429
1430 overflow = 0;
1431
1432 /* ltmpx contains the generalized value of the old real "x" equation */
1433 LCMPLXtimesreal(lparm2,lold.y,ltmp); /* tmp = (C * old.y) */
1434 LCMPLXtrig1(ltmp,ltmpx); /* tmpx = trig1(tmp) */
1435 FIX_OVERFLOW(ltmpx);
1436 ltmpx.x += lold.y; /* tmpx = old.y + trig1(tmp) */
1437 LCMPLXtrig0(ltmpx,ltmp); /* tmp = trig0(tmpx) */
1438 FIX_OVERFLOW(ltmp);
1439 LCMPLXmult(ltmp,lparm,ltmpx); /* tmpx = tmp * h */
1440
1441 /* ltmpy contains the generalized value of the old real "y" equation */
1442 LCMPLXtimesreal(lparm2,lold.x,ltmp); /* tmp = (C * old.x) */
1443 LCMPLXtrig3(ltmp,ltmpy); /* tmpy = trig3(tmp) */
1444 FIX_OVERFLOW(ltmpy);
1445 ltmpy.x += lold.x; /* tmpy = old.x + trig1(tmp) */
1446 LCMPLXtrig2(ltmpy,ltmp); /* tmp = trig2(tmpy) */
1447 FIX_OVERFLOW(ltmp);
1448 LCMPLXmult(ltmp,lparm,ltmpy); /* tmpy = tmp * h */
1449
1450 lnew.x = lold.x - ltmpx.x - ltmpy.y;
1451 lnew.y = lold.y - ltmpy.x - ltmpx.y;
1452
1453 if(plot == noplot)
1454 {
1455 iplot_orbit(lnew.x,lnew.y,1+row%colors);
1456 lold = lnew;
1457 }
1458 ltempsqrx = lsqr(lnew.x);
1459 ltempsqry = lsqr(lnew.y);
1460 lmagnitud = ltempsqrx + ltempsqry;
1461 if (lmagnitud >= llimit || lmagnitud < 0
1462 || labs(lnew.x) > llimit2
1463 || labs(lnew.y) > llimit2)
1464 return(1);
1465 lold = lnew;
1466 return(0);
1469 #endif
1470 }
1471
1472 int MarksCplxMand(void)
1473 {
1474 tmp.x = tempsqrx - tempsqry;
1475 tmp.y = 2*old.x*old.y;
1476 FPUcplxmul(&tmp, &coefficient, &new);
1477 new.x += floatparm->x;
1478 new.y += floatparm->y;
1479 return(floatbailout());
1480 }
1481
1482 int SpiderfpFractal(void)
1483 {
1484 /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
1485 new.x = tempsqrx - tempsqry + tmp.x;
1486 new.y = 2 * old.x * old.y + tmp.y;
1487 tmp.x = tmp.x/2 + new.x;
1488 tmp.y = tmp.y/2 + new.y;
1489 return(floatbailout());
1490 }
1491
1492 int
1493 SpiderFractal(void)
1494 {
1495 #ifndef XFRACT
1496 /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
1497 lnew.x = ltempsqrx - ltempsqry + ltmp.x;
1498 lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y;
1499 ltmp.x = (ltmp.x >> 1) + lnew.x;
1500 ltmp.y = (ltmp.y >> 1) + lnew.y;
1501 return(longbailout());
1504 #endif
1505 }
1506
1507 int
1508 TetratefpFractal(void)
1509 {
1510 /* Tetrate(XAXIS) { c=z=pixel: z=c^z, |z|<=(P1+3) } */
1511 new = ComplexPower(*floatparm,old);
1512 return(floatbailout());
1513 }
1514
1515 int
1516 ZXTrigPlusZFractal(void)
1517 {
1518 #ifndef XFRACT
1519 /* z = (p1*z*trig(z))+p2*z */
1520 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(old) */
1521 LCMPLXmult(lparm,ltmp,ltmp); /* ltmp = p1*trig(old) */
1522 LCMPLXmult(lold,ltmp,ltmp2); /* ltmp2 = p1*old*trig(old) */
1523 LCMPLXmult(lparm2,lold,ltmp); /* ltmp = p2*old */
1524 LCMPLXadd(ltmp2,ltmp,lnew); /* lnew = p1*trig(old) + p2*old */
1525 return(longbailout());
1528 #endif
1529 }
1530
1531 int
1532 ScottZXTrigPlusZFractal(void)
1533 {
1534 #ifndef XFRACT
1535 /* z = (z*trig(z))+z */
1536 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(old) */
1537 LCMPLXmult(lold,ltmp,lnew); /* lnew = old*trig(old) */
1538 LCMPLXadd(lnew,lold,lnew); /* lnew = trig(old) + old */
1539 return(longbailout());
1542 #endif
1543 }
1544
1545 int
1546 SkinnerZXTrigSubZFractal(void)
1547 {
1548 #ifndef XFRACT
1549 /* z = (z*trig(z))-z */
1550 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(old) */
1551 LCMPLXmult(lold,ltmp,lnew); /* lnew = old*trig(old) */
1552 LCMPLXsub(lnew,lold,lnew); /* lnew = trig(old) - old */
1553 return(longbailout());
1556 #endif
1557 }
1558
1559 int
1560 ZXTrigPlusZfpFractal(void)
1561 {
1562 /* z = (p1*z*trig(z))+p2*z */
1563 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
1564 CMPLXmult(parm,tmp,tmp); /* tmp = p1*trig(old) */
1565 CMPLXmult(old,tmp,tmp2); /* tmp2 = p1*old*trig(old) */
1566 CMPLXmult(parm2,old,tmp); /* tmp = p2*old */
1567 CMPLXadd(tmp2,tmp,new); /* new = p1*trig(old) + p2*old */
1568 return(floatbailout());
1569 }
1570
1571 int
1572 ScottZXTrigPlusZfpFractal(void)
1573 {
1574 /* z = (z*trig(z))+z */
1575 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
1576 CMPLXmult(old,tmp,new); /* new = old*trig(old) */
1577 CMPLXadd(new,old,new); /* new = trig(old) + old */
1578 return(floatbailout());
1579 }
1580
1581 int
1582 SkinnerZXTrigSubZfpFractal(void)
1583 {
1584 /* z = (z*trig(z))-z */
1585 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
1586 CMPLXmult(old,tmp,new); /* new = old*trig(old) */
1587 CMPLXsub(new,old,new); /* new = trig(old) - old */
1588 return(floatbailout());
1589 }
1590
1591 int
1592 Sqr1overTrigFractal(void)
1593 {
1594 #ifndef XFRACT
1595 /* z = sqr(1/trig(z)) */
1596 LCMPLXtrig0(lold,lold);
1597 LCMPLXrecip(lold,lold);
1598 LCMPLXsqr(lold,lnew);
1599 return(longbailout());
1602 #endif
1603 }
1604
1605 int
1606 Sqr1overTrigfpFractal(void)
1607 {
1608 /* z = sqr(1/trig(z)) */
1609 CMPLXtrig0(old,old);
1610 CMPLXrecip(old,old);
1611 CMPLXsqr(old,new);
1612 return(floatbailout());
1613 }
1614
1615 int
1616 TrigPlusTrigFractal(void)
1617 {
1618 #ifndef XFRACT
1619 /* z = trig(0,z)*p1+trig1(z)*p2 */
1620 LCMPLXtrig0(lold,ltmp);
1621 LCMPLXmult(lparm,ltmp,ltmp);
1622 LCMPLXtrig1(lold,ltmp2);
1623 LCMPLXmult(lparm2,ltmp2,lold);
1624 LCMPLXadd(ltmp,lold,lnew);
1625 return(longbailout());
1628 #endif
1629 }
1630
1631 int
1632 TrigPlusTrigfpFractal(void)
1633 {
1634 /* z = trig0(z)*p1+trig1(z)*p2 */
1635 CMPLXtrig0(old,tmp);
1636 CMPLXmult(parm,tmp,tmp);
1637 CMPLXtrig1(old,old);
1638 CMPLXmult(parm2,old,old);
1639 CMPLXadd(tmp,old,new);
1640 return(floatbailout());
1641 }
1642
1643 /* The following four fractals are based on the idea of parallel
1644 or alternate calculations. The shift is made when the mod
1645 reaches a given value. JCO 5/6/92 */
1646
1647 int
1648 LambdaTrigOrTrigFractal(void)
1649 {
1650 #ifndef XFRACT
1651 /* z = trig0(z)*p1 if mod(old) < p2.x and
1652 trig1(z)*p1 if mod(old) >= p2.x */
1653 if ((LCMPLXmod(lold)) < lparm2.x){
1654 LCMPLXtrig0(lold,ltmp);
1655 LCMPLXmult(*longparm,ltmp,lnew);}
1656 else{
1657 LCMPLXtrig1(lold,ltmp);
1658 LCMPLXmult(*longparm,ltmp,lnew);}
1659 return(longbailout());
1662 #endif
1663 }
1664
1665 int
1666 LambdaTrigOrTrigfpFractal(void)
1667 {
1668 /* z = trig0(z)*p1 if mod(old) < p2.x and
1669 trig1(z)*p1 if mod(old) >= p2.x */
1670 if (CMPLXmod(old) < parm2.x){
1671 CMPLXtrig0(old,old);
1672 FPUcplxmul(floatparm,&old,&new);}
1673 else{
1674 CMPLXtrig1(old,old);
1675 FPUcplxmul(floatparm,&old,&new);}
1676 return(floatbailout());
1677 }
1678
1679 int
1680 JuliaTrigOrTrigFractal(void)
1681 {
1682 #ifndef XFRACT
1683 /* z = trig0(z)+p1 if mod(old) < p2.x and
1684 trig1(z)+p1 if mod(old) >= p2.x */
1685 if (LCMPLXmod(lold) < lparm2.x){
1686 LCMPLXtrig0(lold,ltmp);
1687 LCMPLXadd(*longparm,ltmp,lnew);}
1688 else{
1689 LCMPLXtrig1(lold,ltmp);
1690 LCMPLXadd(*longparm,ltmp,lnew);}
1691 return(longbailout());
1694 #endif
1695 }
1696
1697 int
1698 JuliaTrigOrTrigfpFractal(void)
1699 {
1700 /* z = trig0(z)+p1 if mod(old) < p2.x and
1701 trig1(z)+p1 if mod(old) >= p2.x */
1702 if (CMPLXmod(old) < parm2.x){
1703 CMPLXtrig0(old,old);
1704 CMPLXadd(*floatparm,old,new);}
1705 else{
1706 CMPLXtrig1(old,old);
1707 CMPLXadd(*floatparm,old,new);}
1708 return(floatbailout());
1709 }
1710
1711 int AplusOne, Ap1deg;
1712 struct MP mpAplusOne, mpAp1deg;
1713 struct MPC mpctmpparm;
1714
1715 #if (_MSC_VER >= 700)
1716 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
1717 #endif
1718
1719 int MPCHalleyFractal(void)
1720 {
1721 #ifndef XFRACT
1722 /* X(X^a - 1) = 0, Halley Map */
1723 /* a = parm.x, relaxation coeff. = parm.y, epsilon = parm2.x */
1724
1725 int ihal;
1726 struct MPC mpcXtoAlessOne, mpcXtoA;
1727 struct MPC mpcXtoAplusOne; /* a-1, a, a+1 */
1728 struct MPC mpcFX, mpcF1prime, mpcF2prime, mpcHalnumer1;
1729 struct MPC mpcHalnumer2, mpcHaldenom, mpctmp;
1730
1731 MPOverflow = 0;
1732 mpcXtoAlessOne.x = mpcold.x;
1733 mpcXtoAlessOne.y = mpcold.y;
1734 for(ihal=2; ihal<degree; ihal++) {
1735 mpctmp.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
1736 mpctmp.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
1737 mpcXtoAlessOne.x = mpctmp.x;
1738 mpcXtoAlessOne.y = mpctmp.y;
1739 }
1740 mpcXtoA.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
1741 mpcXtoA.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
1742 mpcXtoAplusOne.x = *pMPsub(*pMPmul(mpcXtoA.x,mpcold.x),*pMPmul(mpcXtoA.y,mpcold.y));
1743 mpcXtoAplusOne.y = *pMPadd(*pMPmul(mpcXtoA.x,mpcold.y),*pMPmul(mpcXtoA.y,mpcold.x));
1744
1745 mpcFX.x = *pMPsub(mpcXtoAplusOne.x, mpcold.x);
1746 mpcFX.y = *pMPsub(mpcXtoAplusOne.y, mpcold.y); /* FX = X^(a+1) - X = F */
1747
1748 mpcF2prime.x = *pMPmul(mpAp1deg, mpcXtoAlessOne.x); /* mpAp1deg in setup */
1749 mpcF2prime.y = *pMPmul(mpAp1deg, mpcXtoAlessOne.y); /* F" */
1750
1751 mpcF1prime.x = *pMPsub(*pMPmul(mpAplusOne, mpcXtoA.x), mpone);
1752 mpcF1prime.y = *pMPmul(mpAplusOne, mpcXtoA.y); /* F' */
1753
1754 mpctmp.x = *pMPsub(*pMPmul(mpcF2prime.x,mpcFX.x),*pMPmul(mpcF2prime.y,mpcFX.y));
1755 mpctmp.y = *pMPadd(*pMPmul(mpcF2prime.x,mpcFX.y),*pMPmul(mpcF2prime.y,mpcFX.x));
1756 /* F * F" */
1757
1758 mpcHaldenom.x = *pMPadd(mpcF1prime.x, mpcF1prime.x);
1759 mpcHaldenom.y = *pMPadd(mpcF1prime.y, mpcF1prime.y); /* 2 * F' */
1760
1761 mpcHalnumer1 = MPCdiv(mpctmp, mpcHaldenom); /* F"F/2F' */
1762 mpctmp.x = *pMPsub(mpcF1prime.x, mpcHalnumer1.x);
1763 mpctmp.y = *pMPsub(mpcF1prime.y, mpcHalnumer1.y); /* F' - F"F/2F' */
1764 mpcHalnumer2 = MPCdiv(mpcFX, mpctmp);
1765
1766 mpctmp = MPCmul(mpctmpparm, mpcHalnumer2); /* mpctmpparm is */
1767 /* relaxation coef. */
1768 #if 0
1777 #endif
1778 mpcnew = MPCsub(mpcold, mpctmp);
1779 new = MPC2cmplx(mpcnew);
1780 return(MPCHalleybailout()||MPOverflow);
1783 #endif
1784 }
1785 #if (_MSC_VER >= 700)
1786 #pragma code_seg () /* back to normal segment */
1787 #endif
1788
1789 int
1790 HalleyFractal(void)
1791 {
1792 /* X(X^a - 1) = 0, Halley Map */
1793 /* a = parm.x = degree, relaxation coeff. = parm.y, epsilon = parm2.x */
1794
1795 int ihal;
1796 _CMPLX XtoAlessOne, XtoA, XtoAplusOne; /* a-1, a, a+1 */
1797 _CMPLX FX, F1prime, F2prime, Halnumer1, Halnumer2, Haldenom;
1798 _CMPLX relax;
1799
1800 XtoAlessOne = old;
1801 for(ihal=2; ihal<degree; ihal++) {
1802 FPUcplxmul(&old, &XtoAlessOne, &XtoAlessOne);
1803 }
1804 FPUcplxmul(&old, &XtoAlessOne, &XtoA);
1805 FPUcplxmul(&old, &XtoA, &XtoAplusOne);
1806
1807 CMPLXsub(XtoAplusOne, old, FX); /* FX = X^(a+1) - X = F */
1808 F2prime.x = Ap1deg * XtoAlessOne.x; /* Ap1deg in setup */
1809 F2prime.y = Ap1deg * XtoAlessOne.y; /* F" */
1810
1811 F1prime.x = AplusOne * XtoA.x - 1.0;
1812 F1prime.y = AplusOne * XtoA.y; /* F' */
1813
1814 FPUcplxmul(&F2prime, &FX, &Halnumer1); /* F * F" */
1815 Haldenom.x = F1prime.x + F1prime.x;
1816 Haldenom.y = F1prime.y + F1prime.y; /* 2 * F' */
1817
1818 FPUcplxdiv(&Halnumer1, &Haldenom, &Halnumer1); /* F"F/2F' */
1819 CMPLXsub(F1prime, Halnumer1, Halnumer2); /* F' - F"F/2F' */
1820 FPUcplxdiv(&FX, &Halnumer2, &Halnumer2);
1821 /* parm.y is relaxation coef. */
1822 /* new.x = old.x - (parm.y * Halnumer2.x);
1823 new.y = old.y - (parm.y * Halnumer2.y); */
1824 relax.x = parm.y;
1825 relax.y = param[3];
1826 FPUcplxmul(&relax, &Halnumer2, &Halnumer2);
1827 new.x = old.x - Halnumer2.x;
1828 new.y = old.y - Halnumer2.y;
1829 return(Halleybailout());
1830 }
1831
1832 int
1833 LongPhoenixFractal(void)
1834 {
1835 #ifndef XFRACT
1836 /* z(n+1) = z(n)^2 + p + qy(n), y(n+1) = z(n) */
1837 ltmp.x = multiply(lold.x, lold.y, bitshift);
1838 lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(longparm->y,ltmp2.x,bitshift);
1839 lnew.y = (ltmp.x + ltmp.x) + multiply(longparm->y,ltmp2.y,bitshift);
1840 ltmp2 = lold; /* set ltmp2 to Y value */
1841 return(longbailout());
1844 #endif
1845 }
1846
1847 int
1848 PhoenixFractal(void)
1849 {
1850 /* z(n+1) = z(n)^2 + p + qy(n), y(n+1) = z(n) */
1851 tmp.x = old.x * old.y;
1852 new.x = tempsqrx - tempsqry + floatparm->x + (floatparm->y * tmp2.x);
1853 new.y = (tmp.x + tmp.x) + (floatparm->y * tmp2.y);
1854 tmp2 = old; /* set tmp2 to Y value */
1855 return(floatbailout());
1856 }
1857
1858 int
1859 LongPhoenixFractalcplx(void)
1860 {
1861 #ifndef XFRACT
1862 /* z(n+1) = z(n)^2 + p + qy(n), y(n+1) = z(n) */
1863 ltmp.x = multiply(lold.x, lold.y, bitshift);
1864 lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(lparm2.x,ltmp2.x,bitshift)-multiply(lparm2.y,ltmp2.y,bitshift);
1865 lnew.y = (ltmp.x + ltmp.x)+longparm->y+multiply(lparm2.x,ltmp2.y,bitshift)+multiply(lparm2.y,ltmp2.x,bitshift);
1866 ltmp2 = lold; /* set ltmp2 to Y value */
1867 return(longbailout());
1870 #endif
1871 }
1872
1873 int
1874 PhoenixFractalcplx(void)
1875 {
1876 /* z(n+1) = z(n)^2 + p1 + p2*y(n), y(n+1) = z(n) */
1877 tmp.x = old.x * old.y;
1878 new.x = tempsqrx - tempsqry + floatparm->x + (parm2.x * tmp2.x) - (parm2.y * tmp2.y);
1879 new.y = (tmp.x + tmp.x) + floatparm->y + (parm2.x * tmp2.y) + (parm2.y * tmp2.x);
1880 tmp2 = old; /* set tmp2 to Y value */
1881 return(floatbailout());
1882 }
1883
1884 int
1885 LongPhoenixPlusFractal(void)
1886 {
1887 #ifndef XFRACT
1888 /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
1889 int i;
1890 _LCMPLX loldplus, lnewminus;
1891 loldplus = lold;
1892 ltmp = lold;
1893 for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
1894 LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
1895 }
1896 loldplus.x += longparm->x;
1897 LCMPLXmult(ltmp, loldplus, lnewminus);
1898 lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
1899 lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
1900 ltmp2 = lold; /* set ltmp2 to Y value */
1901 return(longbailout());
1904 #endif
1905 }
1906
1907 int
1908 PhoenixPlusFractal(void)
1909 {
1910 /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
1911 int i;
1912 _CMPLX oldplus, newminus;
1913 oldplus = old;
1914 tmp = old;
1915 for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
1916 FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
1917 }
1918 oldplus.x += floatparm->x;
1919 FPUcplxmul(&tmp, &oldplus, &newminus);
1920 new.x = newminus.x + (floatparm->y * tmp2.x);
1921 new.y = newminus.y + (floatparm->y * tmp2.y);
1922 tmp2 = old; /* set tmp2 to Y value */
1923 return(floatbailout());
1924 }
1925
1926 int
1927 LongPhoenixMinusFractal(void)
1928 {
1929 #ifndef XFRACT
1930 /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
1931 int i;
1932 _LCMPLX loldsqr, lnewminus;
1933 LCMPLXmult(lold,lold,loldsqr);
1934 ltmp = lold;
1935 for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
1936 LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
1937 }
1938 loldsqr.x += longparm->x;
1939 LCMPLXmult(ltmp, loldsqr, lnewminus);
1940 lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
1941 lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
1942 ltmp2 = lold; /* set ltmp2 to Y value */
1943 return(longbailout());
1946 #endif
1947 }
1948
1949 int
1950 PhoenixMinusFractal(void)
1951 {
1952 /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
1953 int i;
1954 _CMPLX oldsqr, newminus;
1955 FPUcplxmul(&old, &old, &oldsqr);
1956 tmp = old;
1957 for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
1958 FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
1959 }
1960 oldsqr.x += floatparm->x;
1961 FPUcplxmul(&tmp, &oldsqr, &newminus);
1962 new.x = newminus.x + (floatparm->y * tmp2.x);
1963 new.y = newminus.y + (floatparm->y * tmp2.y);
1964 tmp2 = old; /* set tmp2 to Y value */
1965 return(floatbailout());
1966 }
1967
1968 int
1969 LongPhoenixCplxPlusFractal(void)
1970 {
1971 #ifndef XFRACT
1972 /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
1973 int i;
1974 _LCMPLX loldplus, lnewminus;
1975 loldplus = lold;
1976 ltmp = lold;
1977 for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
1978 LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
1979 }
1980 loldplus.x += longparm->x;
1981 loldplus.y += longparm->y;
1982 LCMPLXmult(ltmp, loldplus, lnewminus);
1983 LCMPLXmult(lparm2, ltmp2, ltmp);
1984 lnew.x = lnewminus.x + ltmp.x;
1985 lnew.y = lnewminus.y + ltmp.y;
1986 ltmp2 = lold; /* set ltmp2 to Y value */
1987 return(longbailout());
1990 #endif
1991 }
1992
1993 int
1994 PhoenixCplxPlusFractal(void)
1995 {
1996 /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
1997 int i;
1998 _CMPLX oldplus, newminus;
1999 oldplus = old;
2000 tmp = old;
2001 for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
2002 FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
2003 }
2004 oldplus.x += floatparm->x;
2005 oldplus.y += floatparm->y;
2006 FPUcplxmul(&tmp, &oldplus, &newminus);
2007 FPUcplxmul(&parm2, &tmp2, &tmp);
2008 new.x = newminus.x + tmp.x;
2009 new.y = newminus.y + tmp.y;
2010 tmp2 = old; /* set tmp2 to Y value */
2011 return(floatbailout());
2012 }
2013
2014 int
2015 LongPhoenixCplxMinusFractal(void)
2016 {
2017 #ifndef XFRACT
2018 /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
2019 int i;
2020 _LCMPLX loldsqr, lnewminus;
2021 LCMPLXmult(lold,lold,loldsqr);
2022 ltmp = lold;
2023 for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
2024 LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
2025 }
2026 loldsqr.x += longparm->x;
2027 loldsqr.y += longparm->y;
2028 LCMPLXmult(ltmp, loldsqr, lnewminus);
2029 LCMPLXmult(lparm2, ltmp2, ltmp);
2030 lnew.x = lnewminus.x + ltmp.x;
2031 lnew.y = lnewminus.y + ltmp.y;
2032 ltmp2 = lold; /* set ltmp2 to Y value */
2033 return(longbailout());
2036 #endif
2037 }
2038
2039 int
2040 PhoenixCplxMinusFractal(void)
2041 {
2042 /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
2043 int i;
2044 _CMPLX oldsqr, newminus;
2045 FPUcplxmul(&old, &old, &oldsqr);
2046 tmp = old;
2047 for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
2048 FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
2049 }
2050 oldsqr.x += floatparm->x;
2051 oldsqr.y += floatparm->y;
2052 FPUcplxmul(&tmp, &oldsqr, &newminus);
2053 FPUcplxmul(&parm2, &tmp2, &tmp);
2054 new.x = newminus.x + tmp.x;
2055 new.y = newminus.y + tmp.y;
2056 tmp2 = old; /* set tmp2 to Y value */
2057 return(floatbailout());
2058 }
2059
2060 int
2061 ScottTrigPlusTrigFractal(void)
2062 {
2063 #ifndef XFRACT
2064 /* z = trig0(z)+trig1(z) */
2065 LCMPLXtrig0(lold,ltmp);
2066 LCMPLXtrig1(lold,lold);
2067 LCMPLXadd(ltmp,lold,lnew);
2068 return(longbailout());
2071 #endif
2072 }
2073
2074 int
2075 ScottTrigPlusTrigfpFractal(void)
2076 {
2077 /* z = trig0(z)+trig1(z) */
2078 CMPLXtrig0(old,tmp);
2079 CMPLXtrig1(old,tmp2);
2080 CMPLXadd(tmp,tmp2,new);
2081 return(floatbailout());
2082 }
2083
2084 int
2085 SkinnerTrigSubTrigFractal(void)
2086 {
2087 #ifndef XFRACT
2088 /* z = trig(0,z)-trig1(z) */
2089 LCMPLXtrig0(lold,ltmp);
2090 LCMPLXtrig1(lold,ltmp2);
2091 LCMPLXsub(ltmp,ltmp2,lnew);
2092 return(longbailout());
2095 #endif
2096 }
2097
2098 int
2099 SkinnerTrigSubTrigfpFractal(void)
2100 {
2101 /* z = trig0(z)-trig1(z) */
2102 CMPLXtrig0(old,tmp);
2103 CMPLXtrig1(old,tmp2);
2104 CMPLXsub(tmp,tmp2,new);
2105 return(floatbailout());
2106 }
2107
2108 int
2109 TrigXTrigfpFractal(void)
2110 {
2111 /* z = trig0(z)*trig1(z) */
2112 CMPLXtrig0(old,tmp);
2113 CMPLXtrig1(old,old);
2114 CMPLXmult(tmp,old,new);
2115 return(floatbailout());
2116 }
2117
2118 #ifndef XFRACT
2119 /* call float version of fractal if integer math overflow */
2120 static int TryFloatFractal(int (*fpFractal)(void))
2121 {
2122 overflow=0;
2123 /* lold had better not be changed! */
2124 old.x = lold.x; old.x /= fudge;
2125 old.y = lold.y; old.y /= fudge;
2126 tempsqrx = sqr(old.x);
2127 tempsqry = sqr(old.y);
2128 fpFractal();
2129 if (save_release < 1900) { /* for backwards compatibility */
2130 lnew.x = (long)(new.x/fudge); /* this error has been here a long time */
2131 lnew.y = (long)(new.y/fudge);
2132 } else {
2133 lnew.x = (long)(new.x*fudge);
2134 lnew.y = (long)(new.y*fudge);
2135 }
2136 return(0);
2137 }
2138 #endif
2139
2140 int
2141 TrigXTrigFractal(void)
2142 {
2143 #ifndef XFRACT
2144 _LCMPLX ltmp2;
2145 /* z = trig0(z)*trig1(z) */
2146 LCMPLXtrig0(lold,ltmp);
2147 LCMPLXtrig1(lold,ltmp2);
2148 LCMPLXmult(ltmp,ltmp2,lnew);
2149 if(overflow)
2150 TryFloatFractal(TrigXTrigfpFractal);
2151 return(longbailout());
2154 #endif
2155 }
2156
2157 /********************************************************************/
2158 /* Next six orbit functions are one type - extra functions are */
2159 /* special cases written for speed. */
2160 /********************************************************************/
2161
2162 int
2163 TrigPlusSqrFractal(void) /* generalization of Scott and Skinner types */
2164 {
2165 #ifndef XFRACT
2166 /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
2167 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
2168 LCMPLXmult(lparm,ltmp,lnew); /* lnew = lparm*trig(lold) */
2169 LCMPLXsqr_old(ltmp); /* ltmp = sqr(lold) */
2170 LCMPLXmult(lparm2,ltmp,ltmp);/* ltmp = lparm2*sqr(lold) */
2171 LCMPLXadd(lnew,ltmp,lnew); /* lnew = lparm*trig(lold)+lparm2*sqr(lold) */
2172 return(longbailout());
2175 #endif
2176 }
2177
2178 int
2179 TrigPlusSqrfpFractal(void) /* generalization of Scott and Skinner types */
2180 {
2181 /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
2182 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
2183 CMPLXmult(parm,tmp,new); /* new = parm*trig(old) */
2184 CMPLXsqr_old(tmp); /* tmp = sqr(old) */
2185 CMPLXmult(parm2,tmp,tmp2); /* tmp = parm2*sqr(old) */
2186 CMPLXadd(new,tmp2,new); /* new = parm*trig(old)+parm2*sqr(old) */
2187 return(floatbailout());
2188 }
2189
2190 int
2191 ScottTrigPlusSqrFractal(void)
2192 {
2193 #ifndef XFRACT
2194 /* { z=pixel: z=trig(z)+sqr(z), |z|<BAILOUT } */
2195 LCMPLXtrig0(lold,lnew); /* lnew = trig(lold) */
2196 LCMPLXsqr_old(ltmp); /* lold = sqr(lold) */
2197 LCMPLXadd(ltmp,lnew,lnew); /* lnew = trig(lold)+sqr(lold) */
2198 return(longbailout());
2201 #endif
2202 }
2203
2204 int
2205 ScottTrigPlusSqrfpFractal(void) /* float version */
2206 {
2207 /* { z=pixel: z=sin(z)+sqr(z), |z|<BAILOUT } */
2208 CMPLXtrig0(old,new); /* new = trig(old) */
2209 CMPLXsqr_old(tmp); /* tmp = sqr(old) */
2210 CMPLXadd(new,tmp,new); /* new = trig(old)+sqr(old) */
2211 return(floatbailout());
2212 }
2213
2214 int
2215 SkinnerTrigSubSqrFractal(void)
2216 {
2217 #ifndef XFRACT
2218 /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
2219 LCMPLXtrig0(lold,lnew); /* lnew = trig(lold) */
2220 LCMPLXsqr_old(ltmp); /* lold = sqr(lold) */
2221 LCMPLXsub(lnew,ltmp,lnew); /* lnew = trig(lold)-sqr(lold) */
2222 return(longbailout());
2225 #endif
2226 }
2227
2228 int
2229 SkinnerTrigSubSqrfpFractal(void)
2230 {
2231 /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
2232 CMPLXtrig0(old,new); /* new = trig(old) */
2233 CMPLXsqr_old(tmp); /* old = sqr(old) */
2234 CMPLXsub(new,tmp,new); /* new = trig(old)-sqr(old) */
2235 return(floatbailout());
2236 }
2237
2238 int
2239 TrigZsqrdfpFractal(void)
2240 {
2241 /* { z=pixel: z=trig(z*z), |z|<TEST } */
2242 CMPLXsqr_old(tmp);
2243 CMPLXtrig0(tmp,new);
2244 return(floatbailout());
2245 }
2246
2247 int
2248 TrigZsqrdFractal(void) /* this doesn't work very well */
2249 {
2250 #ifndef XFRACT
2251 /* { z=pixel: z=trig(z*z), |z|<TEST } */
2252 long l16triglim_2 = 8L << 15;
2253 LCMPLXsqr_old(ltmp);
2254 if((labs(ltmp.x) > l16triglim_2 || labs(ltmp.y) > l16triglim_2) &&
2255 save_release > 1900)
2256 overflow = 1;
2257 else
2258 {
2259 LCMPLXtrig0(ltmp,lnew);
2260 }
2261 if(overflow)
2262 TryFloatFractal(TrigZsqrdfpFractal);
2263 return(longbailout());
2266 #endif
2267 }
2268
2269 int
2270 SqrTrigFractal(void)
2271 {
2272 #ifndef XFRACT
2273 /* { z=pixel: z=sqr(trig(z)), |z|<TEST} */
2274 LCMPLXtrig0(lold,ltmp);
2275 LCMPLXsqr(ltmp,lnew);
2276 return(longbailout());
2279 #endif
2280 }
2281
2282 int
2283 SqrTrigfpFractal(void)
2284 {
2285 /* SZSB(XYAXIS) { z=pixel, TEST=(p1+3): z=sin(z)*sin(z), |z|<TEST} */
2286 CMPLXtrig0(old,tmp);
2287 CMPLXsqr(tmp,new);
2288 return(floatbailout());
2289 }
2290
2291 int
2292 Magnet1Fractal(void) /* Z = ((Z**2 + C - 1)/(2Z + C - 2))**2 */
2293 { /* In "Beauty of Fractals", code by Kev Allen. */
2294 _CMPLX top, bot, tmp;
2295 double div;
2296
2297 top.x = tempsqrx - tempsqry + floatparm->x - 1; /* top = Z**2+C-1 */
2298 top.y = old.x * old.y;
2299 top.y = top.y + top.y + floatparm->y;
2300
2301 bot.x = old.x + old.x + floatparm->x - 2; /* bot = 2*Z+C-2 */
2302 bot.y = old.y + old.y + floatparm->y;
2303
2304 div = bot.x*bot.x + bot.y*bot.y; /* tmp = top/bot */
2305 if (div < FLT_MIN) return(1);
2306 tmp.x = (top.x*bot.x + top.y*bot.y)/div;
2307 tmp.y = (top.y*bot.x - top.x*bot.y)/div;
2308
2309 new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y); /* Z = tmp**2 */
2310 new.y = tmp.x * tmp.y;
2311 new.y += new.y;
2312
2313 return(floatbailout());
2314 }
2315
2316 int
2317 Magnet2Fractal(void) /* Z = ((Z**3 + 3(C-1)Z + (C-1)(C-2) ) / */
2318 /* (3Z**2 + 3(C-2)Z + (C-1)(C-2)+1) )**2 */
2319 { /* In "Beauty of Fractals", code by Kev Allen. */
2320 _CMPLX top, bot, tmp;
2321 double div;
2322
2323 top.x = old.x * (tempsqrx-tempsqry-tempsqry-tempsqry + T_Cm1.x)
2324 - old.y * T_Cm1.y + T_Cm1Cm2.x;
2325 top.y = old.y * (tempsqrx+tempsqrx+tempsqrx-tempsqry + T_Cm1.x)
2326 + old.x * T_Cm1.y + T_Cm1Cm2.y;
2327
2328 bot.x = tempsqrx - tempsqry;
2329 bot.x = bot.x + bot.x + bot.x
2330 + old.x * T_Cm2.x - old.y * T_Cm2.y
2331 + T_Cm1Cm2.x + 1.0;
2332 bot.y = old.x * old.y;
2333 bot.y += bot.y;
2334 bot.y = bot.y + bot.y + bot.y
2335 + old.x * T_Cm2.y + old.y * T_Cm2.x
2336 + T_Cm1Cm2.y;
2337
2338 div = bot.x*bot.x + bot.y*bot.y; /* tmp = top/bot */
2339 if (div < FLT_MIN) return(1);
2340 tmp.x = (top.x*bot.x + top.y*bot.y)/div;
2341 tmp.y = (top.y*bot.x - top.x*bot.y)/div;
2342
2343 new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y); /* Z = tmp**2 */
2344 new.y = tmp.x * tmp.y;
2345 new.y += new.y;
2346
2347 return(floatbailout());
2348 }
2349
2350 int
2351 LambdaTrigFractal(void)
2352 {
2353 #ifndef XFRACT
2354 LONGXYTRIGBAILOUT();
2355 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
2356 LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
2357 lold = lnew;
2358 return(0);
2361 #endif
2362 }
2363
2364 int
2365 LambdaTrigfpFractal(void)
2366 {
2367 FLOATXYTRIGBAILOUT();
2368 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
2369 CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
2370 old = new;
2371 return(0);
2372 }
2373
2374 /* bailouts are different for different trig functions */
2375 int
2376 LambdaTrigFractal1(void)
2377 {
2378 #ifndef XFRACT
2379 LONGTRIGBAILOUT(); /* sin,cos */
2380 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
2381 LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
2382 lold = lnew;
2383 return(0);
2386 #endif
2387 }
2388
2389 int
2390 LambdaTrigfpFractal1(void)
2391 {
2392 FLOATTRIGBAILOUT(); /* sin,cos */
2393 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
2394 CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
2395 old = new;
2396 return(0);
2397 }
2398
2399 int
2400 LambdaTrigFractal2(void)
2401 {
2402 #ifndef XFRACT
2403 LONGHTRIGBAILOUT(); /* sinh,cosh */
2404 LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
2405 LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
2406 lold = lnew;
2407 return(0);
2410 #endif
2411 }
2412
2413 int
2414 LambdaTrigfpFractal2(void)
2415 {
2416 #ifndef XFRACT
2417 FLOATHTRIGBAILOUT(); /* sinh,cosh */
2418 CMPLXtrig0(old,tmp); /* tmp = trig(old) */
2419 CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
2420 old = new;
2421 return(0);
2424 #endif
2425 }
2426
2427 int
2428 ManOWarFractal(void)
2429 {
2430 #ifndef XFRACT
2431 /* From Art Matrix via Lee Skinner */
2432 lnew.x = ltempsqrx - ltempsqry + ltmp.x + longparm->x;
2433 lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y + longparm->y;
2434 ltmp = lold;
2435 return(longbailout());
2438 #endif
2439 }
2440
2441 int
2442 ManOWarfpFractal(void)
2443 {
2444 /* From Art Matrix via Lee Skinner */
2445 /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
2446 new.x = tempsqrx - tempsqry + tmp.x + floatparm->x;
2447 new.y = 2.0 * old.x * old.y + tmp.y + floatparm->y;
2448 tmp = old;
2449 return(floatbailout());
2450 }
2451
2452 /*
2453 MarksMandelPwr (XAXIS) {
2454 z = pixel, c = z ^ (z - 1):
2455 z = c * sqr(z) + pixel,
2456 |z| <= 4
2457 }
2458 */
2459
2460 int
2461 MarksMandelPwrfpFractal(void)
2462 {
2463 CMPLXtrig0(old,new);
2464 CMPLXmult(tmp,new,new);
2465 new.x += floatparm->x;
2466 new.y += floatparm->y;
2467 return(floatbailout());
2468 }
2469
2470 int
2471 MarksMandelPwrFractal(void)
2472 {
2473 #ifndef XFRACT
2474 LCMPLXtrig0(lold,lnew);
2475 LCMPLXmult(ltmp,lnew,lnew);
2476 lnew.x += longparm->x;
2477 lnew.y += longparm->y;
2478 return(longbailout());
2481 #endif
2482 }
2483
2484 /* I was coding Marksmandelpower and failed to use some temporary
2485 variables. The result was nice, and since my name is not on any fractal,
2486 I thought I would immortalize myself with this error!
2487 Tim Wegner */
2488
2489 int
2490 TimsErrorfpFractal(void)
2491 {
2492 CMPLXtrig0(old,new);
2493 new.x = new.x * tmp.x - new.y * tmp.y;
2494 new.y = new.x * tmp.y - new.y * tmp.x;
2495 new.x += floatparm->x;
2496 new.y += floatparm->y;
2497 return(floatbailout());
2498 }
2499
2500 int
2501 TimsErrorFractal(void)
2502 {
2503 #ifndef XFRACT
2504 LCMPLXtrig0(lold,lnew);
2505 lnew.x = multiply(lnew.x,ltmp.x,bitshift)-multiply(lnew.y,ltmp.y,bitshift);
2506 lnew.y = multiply(lnew.x,ltmp.y,bitshift)-multiply(lnew.y,ltmp.x,bitshift);
2507 lnew.x += longparm->x;
2508 lnew.y += longparm->y;
2509 return(longbailout());
2512 #endif
2513 }
2514
2515 int
2516 CirclefpFractal(void)
2517 {
2518 long i;
2519 i = (long)(param[0]*(tempsqrx+tempsqry));
2520 coloriter = i%colors;
2521 return(1);
2522 }
2523 /*
2524 CirclelongFractal()
2525 {
2526 long i;
2527 i = multiply(lparm.x,(ltempsqrx+ltempsqry),bitshift);
2528 i = i >> bitshift;
2529 coloriter = i%colors);
2530 return(1);
2531 }
2532 */
2533
2534 /* -------------------------------------------------------------------- */
2535 /* Initialization (once per pixel) routines */
2536 /* -------------------------------------------------------------------- */
2537
2538 #ifdef XFRACT
2555 #endif
2556
2557 int long_julia_per_pixel(void)
2558 {
2559 #ifndef XFRACT
2560 /* integer julia types */
2561 /* lambda */
2562 /* barnsleyj1 */
2563 /* barnsleyj2 */
2564 /* sierpinski */
2565 if(invert)
2566 {
2567 /* invert */
2568 invertz2(&old);
2569
2570 /* watch out for overflow */
2571 if(sqr(old.x)+sqr(old.y) >= 127)
2572 {
2573 old.x = 8; /* value to bail out in one iteration */
2574 old.y = 8;
2575 }
2576
2577 /* convert to fudged longs */
2578 lold.x = (long)(old.x*fudge);
2579 lold.y = (long)(old.y*fudge);
2580 }
2581 else
2582 {
2583 lold.x = lxpixel();
2584 lold.y = lypixel();
2585 }
2586 return(0);
2590 #endif
2591 }
2592
2593 int long_richard8_per_pixel(void)
2594 {
2595 #ifndef XFRACT
2596 long_mandel_per_pixel();
2597 LCMPLXtrig1(*longparm,ltmp);
2598 LCMPLXmult(ltmp,lparm2,ltmp);
2599 return(1);
2602 #endif
2603 }
2604
2605 int long_mandel_per_pixel(void)
2606 {
2607 #ifndef XFRACT
2608 /* integer mandel types */
2609 /* barnsleym1 */
2610 /* barnsleym2 */
2611 linit.x = lxpixel();
2612 if(save_release >= 2004)
2613 linit.y = lypixel();
2614
2615 if(invert)
2616 {
2617 /* invert */
2618 invertz2(&init);
2619
2620 /* watch out for overflow */
2621 if(sqr(init.x)+sqr(init.y) >= 127)
2622 {
2623 init.x = 8; /* value to bail out in one iteration */
2624 init.y = 8;
2625 }
2626
2627 /* convert to fudged longs */
2628 linit.x = (long)(init.x*fudge);
2629 linit.y = (long)(init.y*fudge);
2630 }
2631
2632 if(useinitorbit == 1)
2633 lold = linitorbit;
2634 else
2635 lold = linit;
2636
2637 lold.x += lparm.x; /* initial pertubation of parameters set */
2638 lold.y += lparm.y;
2639 return(1); /* 1st iteration has been done */
2643 #endif
2644 }
2645
2646 int julia_per_pixel(void)
2647 {
2648 /* julia */
2649
2650 if(invert)
2651 {
2652 /* invert */
2653 invertz2(&old);
2654
2655 /* watch out for overflow */
2656 if(bitshift <= 24)
2657 if (sqr(old.x)+sqr(old.y) >= 127)
2658 {
2659 old.x = 8; /* value to bail out in one iteration */
2660 old.y = 8;
2661 }
2662 if(bitshift > 24)
2663 if (sqr(old.x)+sqr(old.y) >= 4.0)
2664 {
2665 old.x = 2; /* value to bail out in one iteration */
2666 old.y = 2;
2667 }
2668
2669 /* convert to fudged longs */
2670 lold.x = (long)(old.x*fudge);
2671 lold.y = (long)(old.y*fudge);
2672 }
2673 else
2674 {
2675 lold.x = lxpixel();
2676 lold.y = lypixel();
2677 }
2678
2679 ltempsqrx = multiply(lold.x, lold.x, bitshift);
2680 ltempsqry = multiply(lold.y, lold.y, bitshift);
2681 ltmp = lold;
2682 return(0);
2683 }
2684
2685 int
2686 marks_mandelpwr_per_pixel(void)
2687 {
2688 #ifndef XFRACT
2689 mandel_per_pixel();
2690 ltmp = lold;
2691 ltmp.x -= fudge;
2692 LCMPLXpwr(lold,ltmp,ltmp);
2693 return(1);
2696 #endif
2697 }
2698
2699 int mandel_per_pixel(void)
2700 {
2701 /* mandel */
2702
2703 if(invert)
2704 {
2705 invertz2(&init);
2706
2707 /* watch out for overflow */
2708 if(bitshift <= 24)
2709 if (sqr(init.x)+sqr(init.y) >= 127)
2710 {
2711 init.x = 8; /* value to bail out in one iteration */
2712 init.y = 8;
2713 }
2714 if(bitshift > 24)
2715 if (sqr(init.x)+sqr(init.y) >= 4)
2716 {
2717 init.x = 2; /* value to bail out in one iteration */
2718 init.y = 2;
2719 }
2720
2721 /* convert to fudged longs */
2722 linit.x = (long)(init.x*fudge);
2723 linit.y = (long)(init.y*fudge);
2724 }
2725 else {
2726 linit.x = lxpixel();
2727 if(save_release >= 2004)
2728 linit.y = lypixel();
2729 }
2730 switch (fractype)
2731 {
2732 case MANDELLAMBDA: /* Critical Value 0.5 + 0.0i */
2733 lold.x = FgHalf;
2734 lold.y = 0;
2735 break;
2736 default:
2737 lold = linit;
2738 break;
2739 }
2740
2741 /* alter init value */
2742 if(useinitorbit == 1)
2743 lold = linitorbit;
2744 else if(useinitorbit == 2)
2745 lold = linit;
2746
2747 if((inside == BOF60 || inside == BOF61) && !nobof)
2748 {
2749 /* kludge to match "Beauty of Fractals" picture since we start
2750 Mandelbrot iteration with init rather than 0 */
2751 lold.x = lparm.x; /* initial pertubation of parameters set */
2752 lold.y = lparm.y;
2753 coloriter = -1;
2754 }
2755 else
2756 {
2757 lold.x += lparm.x; /* initial pertubation of parameters set */
2758 lold.y += lparm.y;
2759 }
2760 ltmp = linit; /* for spider */
2761 ltempsqrx = multiply(lold.x, lold.x, bitshift);
2762 ltempsqry = multiply(lold.y, lold.y, bitshift);
2763 return(1); /* 1st iteration has been done */
2764 }
2765
2766 int marksmandel_per_pixel()
2767 {
2768 #ifndef XFRACT
2769 /* marksmandel */
2770 if(invert)
2771 {
2772 invertz2(&init);
2773
2774 /* watch out for overflow */
2775 if(sqr(init.x)+sqr(init.y) >= 127)
2776 {
2777 init.x = 8; /* value to bail out in one iteration */
2778 init.y = 8;
2779 }
2780
2781 /* convert to fudged longs */
2782 linit.x = (long)(init.x*fudge);
2783 linit.y = (long)(init.y*fudge);
2784 }
2785 else {
2786 linit.x = lxpixel();
2787 if(save_release >= 2004)
2788 linit.y = lypixel();
2789 }
2790
2791 if(useinitorbit == 1)
2792 lold = linitorbit;
2793 else
2794 lold = linit;
2795
2796 lold.x += lparm.x; /* initial pertubation of parameters set */
2797 lold.y += lparm.y;
2798
2799 if(c_exp > 3)
2800 lcpower(&lold,c_exp-1,&lcoefficient,bitshift);
2801 else if(c_exp == 3) {
2802 lcoefficient.x = multiply(lold.x, lold.x, bitshift)
2803 - multiply(lold.y, lold.y, bitshift);
2804 lcoefficient.y = multiply(lold.x, lold.y, bitshiftless1);
2805 }
2806 else if(c_exp == 2)
2807 lcoefficient = lold;
2808 else if(c_exp < 2) {
2809 lcoefficient.x = 1L << bitshift;
2810 lcoefficient.y = 0L;
2811 }
2812
2813 ltempsqrx = multiply(lold.x, lold.x, bitshift);
2814 ltempsqry = multiply(lold.y, lold.y, bitshift);
2815 #endif
2816 return(1); /* 1st iteration has been done */
2817 }
2818
2819 int marksmandelfp_per_pixel()
2820 {
2821 /* marksmandel */
2822
2823 if(invert)
2824 invertz2(&init);
2825 else {
2826 init.x = dxpixel();
2827 if(save_release >= 2004)
2828 init.y = dypixel();
2829 }
2830
2831 if(useinitorbit == 1)
2832 old = initorbit;
2833 else
2834 old = init;
2835
2836 old.x += parm.x; /* initial pertubation of parameters set */
2837 old.y += parm.y;
2838
2839 tempsqrx = sqr(old.x);
2840 tempsqry = sqr(old.y);
2841
2842 if(c_exp > 3)
2843 cpower(&old,c_exp-1,&coefficient);
2844 else if(c_exp == 3) {
2845 coefficient.x = tempsqrx - tempsqry;
2846 coefficient.y = old.x * old.y * 2;
2847 }
2848 else if(c_exp == 2)
2849 coefficient = old;
2850 else if(c_exp < 2) {
2851 coefficient.x = 1.0;
2852 coefficient.y = 0.0;
2853 }
2854
2855 return(1); /* 1st iteration has been done */
2856 }
2857
2858 int
2859 marks_mandelpwrfp_per_pixel(void)
2860 {
2861 mandelfp_per_pixel();
2862 tmp = old;
2863 tmp.x -= 1;
2864 CMPLXpwr(old,tmp,tmp);
2865 return(1);
2866 }
2867
2868 int mandelfp_per_pixel(void)
2869 {
2870 /* floating point mandelbrot */
2871 /* mandelfp */
2872
2873 if(invert)
2874 invertz2(&init);
2875 else {
2876 init.x = dxpixel();
2877 if(save_release >= 2004)
2878 init.y = dypixel();
2879 }
2880 switch (fractype)
2881 {
2882 case MAGNET2M:
2883 FloatPreCalcMagnet2();
2884 case MAGNET1M: /* Crit Val Zero both, but neither */
2885 old.x = old.y = 0.0; /* is of the form f(Z,C) = Z*g(Z)+C */
2886 break;
2887 case MANDELLAMBDAFP: /* Critical Value 0.5 + 0.0i */
2888 old.x = 0.5;
2889 old.y = 0.0;
2890 break;
2891 default:
2892 old = init;
2893 break;
2894 }
2895
2896 /* alter init value */
2897 if(useinitorbit == 1)
2898 old = initorbit;
2899 else if(useinitorbit == 2)
2900 old = init;
2901
2902 if((inside == BOF60 || inside == BOF61) && !nobof)
2903 {
2904 /* kludge to match "Beauty of Fractals" picture since we start
2905 Mandelbrot iteration with init rather than 0 */
2906 old.x = parm.x; /* initial pertubation of parameters set */
2907 old.y = parm.y;
2908 coloriter = -1;
2909 }
2910 else
2911 {
2912 old.x += parm.x;
2913 old.y += parm.y;
2914 }
2915 tmp = init; /* for spider */
2916 tempsqrx = sqr(old.x); /* precalculated value for regular Mandelbrot */
2917 tempsqry = sqr(old.y);
2918 return(1); /* 1st iteration has been done */
2919 }
2920
2921 int juliafp_per_pixel(void)
2922 {
2923 /* floating point julia */
2924 /* juliafp */
2925 if(invert)
2926 invertz2(&old);
2927 else
2928 {
2929 old.x = dxpixel();
2930 old.y = dypixel();
2931 }
2932 tempsqrx = sqr(old.x); /* precalculated value for regular Julia */
2933 tempsqry = sqr(old.y);
2934 tmp = old;
2935 return(0);
2936 }
2937
2938 #if (_MSC_VER >= 700)
2939 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
2940 #endif
2941 int MPCjulia_per_pixel(void)
2942 {
2943 #ifndef XFRACT
2944 /* floating point julia */
2945 /* juliafp */
2946 if(invert)
2947 invertz2(&old);
2948 else
2949 {
2950 old.x = dxpixel();
2951 old.y = dypixel();
2952 }
2953 mpcold.x = *pd2MP(old.x);
2954 mpcold.y = *pd2MP(old.y);
2955 return(0);
2958 #endif
2959 }
2960 #if (_MSC_VER >= 700)
2961 #pragma code_seg () /* back to normal segment */
2962 #endif
2963
2964 int
2965 otherrichard8fp_per_pixel(void)
2966 {
2967 othermandelfp_per_pixel();
2968 CMPLXtrig1(*floatparm,tmp);
2969 CMPLXmult(tmp,parm2,tmp);
2970 return(1);
2971 }
2972
2973 int othermandelfp_per_pixel(void)
2974 {
2975 if(invert)
2976 invertz2(&init);
2977 else {
2978 init.x = dxpixel();
2979 if(save_release >= 2004)
2980 init.y = dypixel();
2981 }
2982
2983 if(useinitorbit == 1)
2984 old = initorbit;
2985 else
2986 old = init;
2987
2988 old.x += parm.x; /* initial pertubation of parameters set */
2989 old.y += parm.y;
2990
2991 return(1); /* 1st iteration has been done */
2992 }
2993
2994 #if (_MSC_VER >= 700)
2995 #pragma code_seg ("mpmath1_text") /* place following in an overlay */
2996 #endif
2997
2998 int MPCHalley_per_pixel(void)
2999 {
3000 #ifndef XFRACT
3001 /* MPC halley */
3002 if(invert)
3003 invertz2(&init);
3004 else {
3005 init.x = dxpixel();
3006 if(save_release >= 2004)
3007 init.y = dypixel();
3008 }
3009
3010 mpcold.x = *pd2MP(init.x);
3011 mpcold.y = *pd2MP(init.y);
3012
3013 return(0);
3016 #endif
3017 }
3018 #if (_MSC_VER >= 700)
3019 #pragma code_seg () /* back to normal segment */
3020 #endif
3021
3022 int Halley_per_pixel(void)
3023 {
3024 if(invert)
3025 invertz2(&init);
3026 else {
3027 init.x = dxpixel();
3028 if(save_release >= 2004)
3029 init.y = dypixel();
3030 }
3031
3032 old = init;
3033
3034 return(0); /* 1st iteration is not done */
3035 }
3036
3037 int otherjuliafp_per_pixel(void)
3038 {
3039 if(invert)
3040 invertz2(&old);
3041 else
3042 {
3043 old.x = dxpixel();
3044 old.y = dypixel();
3045 }
3046 return(0);
3047 }
3048
3049 #if 0
3052 #else
3053 #define Q0 0
3054 #define Q1 0
3055 #endif
3056
3057 int quaternionjulfp_per_pixel(void)
3058 {
3059 old.x = dxpixel();
3060 old.y = dypixel();
3061 floatparm->x = param[4];
3062 floatparm->y = param[5];
3063 qc = param[0];
3064 qci = param[1];
3065 qcj = param[2];
3066 qck = param[3];
3067 return(0);
3068 }
3069
3070 int quaternionfp_per_pixel(void)
3071 {
3072 old.x = 0;
3073 old.y = 0;
3074 floatparm->x = 0;
3075 floatparm->y = 0;
3076 qc = dxpixel();
3077 qci = dypixel();
3078 qcj = param[2];
3079 qck = param[3];
3080 return(0);
3081 }
3082
3083 int MarksCplxMandperp(void)
3084 {
3085 if(invert)
3086 invertz2(&init);
3087 else {
3088 init.x = dxpixel();
3089 if(save_release >= 2004)
3090 init.y = dypixel();
3091 }
3092 old.x = init.x + parm.x; /* initial pertubation of parameters set */
3093 old.y = init.y + parm.y;
3094 tempsqrx = sqr(old.x); /* precalculated value */
3095 tempsqry = sqr(old.y);
3096 coefficient = ComplexPower(init, pwr);
3097 return(1);
3098 }
3099
3100 int long_phoenix_per_pixel(void)
3101 {
3102 #ifndef XFRACT
3103 if(invert)
3104 {
3105 /* invert */
3106 invertz2(&old);
3107
3108 /* watch out for overflow */
3109 if(sqr(old.x)+sqr(old.y) >= 127)
3110 {
3111 old.x = 8; /* value to bail out in one iteration */
3112 old.y = 8;
3113 }
3114
3115 /* convert to fudged longs */
3116 lold.x = (long)(old.x*fudge);
3117 lold.y = (long)(old.y*fudge);
3118 }
3119 else
3120 {
3121 lold.x = lxpixel();
3122 lold.y = lypixel();
3123 }
3124 ltempsqrx = multiply(lold.x, lold.x, bitshift);
3125 ltempsqry = multiply(lold.y, lold.y, bitshift);
3126 ltmp2.x = 0; /* use ltmp2 as the complex Y value */
3127 ltmp2.y = 0;
3128 return(0);
3132 #endif
3133 }
3134
3135 int phoenix_per_pixel(void)
3136 {
3137 if(invert)
3138 invertz2(&old);
3139 else
3140 {
3141 old.x = dxpixel();
3142 old.y = dypixel();
3143 }
3144 tempsqrx = sqr(old.x); /* precalculated value */
3145 tempsqry = sqr(old.y);
3146 tmp2.x = 0; /* use tmp2 as the complex Y value */
3147 tmp2.y = 0;
3148 return(0);
3149 }
3150 int long_mandphoenix_per_pixel(void)
3151 {
3152 #ifndef XFRACT
3153 linit.x = lxpixel();
3154 if(save_release >= 2004)
3155 linit.y = lypixel();
3156
3157 if(invert)
3158 {
3159 /* invert */
3160 invertz2(&init);
3161
3162 /* watch out for overflow */
3163 if(sqr(init.x)+sqr(init.y) >= 127)
3164 {
3165 init.x = 8; /* value to bail out in one iteration */
3166 init.y = 8;
3167 }
3168
3169 /* convert to fudged longs */
3170 linit.x = (long)(init.x*fudge);
3171 linit.y = (long)(init.y*fudge);
3172 }
3173
3174 if(useinitorbit == 1)
3175 lold = linitorbit;
3176 else
3177 lold = linit;
3178
3179 lold.x += lparm.x; /* initial pertubation of parameters set */
3180 lold.y += lparm.y;
3181 ltempsqrx = multiply(lold.x, lold.x, bitshift);
3182 ltempsqry = multiply(lold.y, lold.y, bitshift);
3183 ltmp2.x = 0;
3184 ltmp2.y = 0;
3185 return(1); /* 1st iteration has been done */
3189 #endif
3190 }
3191 int mandphoenix_per_pixel(void)
3192 {
3193 if(invert)
3194 invertz2(&init);
3195 else {
3196 init.x = dxpixel();
3197 if(save_release >= 2004)
3198 init.y = dypixel();
3199 }
3200
3201 if(useinitorbit == 1)
3202 old = initorbit;
3203 else
3204 old = init;
3205
3206 old.x += parm.x; /* initial pertubation of parameters set */
3207 old.y += parm.y;
3208 tempsqrx = sqr(old.x); /* precalculated value */
3209 tempsqry = sqr(old.y);
3210 tmp2.x = 0;
3211 tmp2.y = 0;
3212 return(1); /* 1st iteration has been done */
3213 }
3214
3215 int
3216 QuaternionFPFractal(void)
3217 {
3218 double a0,a1,a2,a3,n0,n1,n2,n3;
3219 a0 = old.x;
3220 a1 = old.y;
3221 a2 = floatparm->x;
3222 a3 = floatparm->y;
3223
3224 n0 = a0*a0-a1*a1-a2*a2-a3*a3 + qc;
3225 n1 = 2*a0*a1 + qci;
3226 n2 = 2*a0*a2 + qcj;
3227 n3 = 2*a0*a3 + qck;
3228 /* Check bailout */
3229 magnitude = a0*a0+a1*a1+a2*a2+a3*a3;
3230 if (magnitude>rqlim) {
3231 return 1;
3232 }
3233 old.x = new.x = n0;
3234 old.y = new.y = n1;
3235 floatparm->x = n2;
3236 floatparm->y = n3;
3237 return(0);
3238 }
3239
3240 int
3241 HyperComplexFPFractal(void)
3242 {
3243 _HCMPLX hold, hnew;
3244 hold.x = old.x;
3245 hold.y = old.y;
3246 hold.z = floatparm->x;
3247 hold.t = floatparm->y;
3248
3249 /* HComplexSqr(&hold,&hnew); */
3250 HComplexTrig0(&hold,&hnew);
3251
3252 hnew.x += qc;
3253 hnew.y += qci;
3254 hnew.z += qcj;
3255 hnew.t += qck;
3256
3257 old.x = new.x = hnew.x;
3258 old.y = new.y = hnew.y;
3259 floatparm->x = hnew.z;
3260 floatparm->y = hnew.t;
3261
3262 /* Check bailout */
3263 magnitude = sqr(old.x)+sqr(old.y)+sqr(floatparm->x)+sqr(floatparm->y);
3264 if (magnitude>rqlim) {
3265 return 1;
3266 }
3267 return(0);
3268 }
3269
3270 int
3271 VLfpFractal(void) /* Beauty of Fractals pp. 125 - 127 */
3272 {
3273 double a, b, ab, half, u, w, xy;
3274
3275 half = param[0] / 2.0;
3276 xy = old.x * old.y;
3277 u = old.x - xy;
3278 w = -old.y + xy;
3279 a = old.x + param[1] * u;
3280 b = old.y + param[1] * w;
3281 ab = a * b;
3282 new.x = old.x + half * (u + (a - ab));
3283 new.y = old.y + half * (w + (-b + ab));
3284 return(floatbailout());
3285 }
3286
3287 int
3288 EscherfpFractal(void) /* Science of Fractal Images pp. 185, 187 */
3289 {
3290 _CMPLX oldtest, newtest, testsqr;
3291 double testsize = 0.0;
3292 long testiter = 0;
3293
3294 new.x = tempsqrx - tempsqry; /* standard Julia with C == (0.0, 0.0i) */
3295 new.y = 2.0 * old.x * old.y;
3296 oldtest.x = new.x * 15.0; /* scale it */
3297 oldtest.y = new.y * 15.0;
3298 testsqr.x = sqr(oldtest.x); /* set up to test with user-specified ... */
3299 testsqr.y = sqr(oldtest.y); /* ... Julia as the target set */
3300 while (testsize <= rqlim && testiter < maxit) /* nested Julia loop */
3301 {
3302 newtest.x = testsqr.x - testsqr.y + param[0];
3303 newtest.y = 2.0 * oldtest.x * oldtest.y + param[1];
3304 testsize = (testsqr.x = sqr(newtest.x)) + (testsqr.y = sqr(newtest.y));
3305 oldtest = newtest;
3306 testiter++;
3307 }
3308 if (testsize > rqlim) return(floatbailout()); /* point not in target set */
3309 else /* make distinct level sets if point stayed in target set */
3310 {
3311 coloriter = ((3L * coloriter) % 255L) + 1L;
3312 return 1;
3313 }
3314 }
3315
3316 /* re-use static roots variable
3317 memory for mandelmix4 */
3318
3319 #define A staticroots[ 0]
3320 #define B staticroots[ 1]
3321 #define C staticroots[ 2]
3322 #define D staticroots[ 3]
3323 #define F staticroots[ 4]
3324 #define G staticroots[ 5]
3325 #define H staticroots[ 6]
3326 #define J staticroots[ 7]
3327 #define K staticroots[ 8]
3328 #define L staticroots[ 9]
3329 #define Z staticroots[10]
3330
3331 int MandelbrotMix4Setup(void)
3332 {
3333 int sign_array = 0;
3334 A.x=param[0]; A.y=0.0; /* a=real(p1), */
3335 B.x=param[1]; B.y=0.0; /* b=imag(p1), */
3336 D.x=param[2]; D.y=0.0; /* d=real(p2), */
3337 F.x=param[3]; F.y=0.0; /* f=imag(p2), */
3338 K.x=param[4]+1.0; K.y=0.0; /* k=real(p3)+1, */
3339 L.x=param[5]+100.0; L.y=0.0; /* l=imag(p3)+100, */
3340 CMPLXrecip(F,G); /* g=1/f, */
3341 CMPLXrecip(D,H); /* h=1/d, */
3342 CMPLXsub(F,B,tmp); /* tmp = f-b */
3343 CMPLXrecip(tmp,J); /* j = 1/(f-b) */
3344 CMPLXneg(A,tmp);
3345 CMPLXmult(tmp,B,tmp); /* z=(-a*b*g*h)^j, */
3346 CMPLXmult(tmp,G,tmp);
3347 CMPLXmult(tmp,H,tmp);
3348
3349 /*
3350 This code kludge attempts to duplicate the behavior
3351 of the parser in determining the sign of zero of the
3352 imaginary part of the argument of the power function. The
3353 reason this is important is that the complex arctangent
3354 returns PI in one case and -PI in the other, depending
3355 on the sign bit of zero, and we wish the results to be
3356 compatible with Jim Muth's mix4 formula using the parser.
3357
3358 First create a number encoding the signs of a, b, g , h. Our
3359 kludge assumes that those signs determine the behavior.
3360 */
3361 if(A.x < 0.0) sign_array += 8;
3362 if(B.x < 0.0) sign_array += 4;
3363 if(G.x < 0.0) sign_array += 2;
3364 if(H.x < 0.0) sign_array += 1;
3365 if(tmp.y == 0.0) /* we know tmp.y IS zero but ... */
3366 {
3367 switch(sign_array)
3368 {
3369 /*
3370 Add to this list the magic numbers of any cases
3371 in which the fractal does not match the formula version
3372 */
3373 case 15: /* 1111 */
3374 case 10: /* 1010 */
3375 case 6: /* 0110 */
3376 case 5: /* 0101 */
3377 case 3: /* 0011 */
3378 case 0: /* 0000 */
3379 tmp.y = -tmp.y; /* swap sign bit */
3380 default: /* do nothing - remaining cases already OK */
3381 ;
3382 }
3383 /* in case our kludge failed, let the user fix it */
3384 if(debugflag == 1012) tmp.y = -tmp.y;
3385 }
3386
3387 CMPLXpwr(tmp,J,tmp); /* note: z is old */
3388 /* in case our kludge failed, let the user fix it */
3389 if(param[6] < 0.0) tmp.y = -tmp.y;
3390
3391 if(bailout == 0)
3392 {
3393 rqlim = L.x;
3394 rqlim2 = rqlim*rqlim;
3395 }
3396 return(1);
3397 }
3398
3399 int MandelbrotMix4fp_per_pixel(void)
3400 {
3401 if(invert)
3402 invertz2(&init);
3403 else {
3404 init.x = dxpixel();
3405 init.y = dypixel();
3406 }
3407 old = tmp;
3408 CMPLXtrig0(init,C); /* c=fn1(pixel): */
3409 return(0); /* 1st iteration has been NOT been done */
3410 }
3411
3412 int
3413 MandelbrotMix4fpFractal(void) /* from formula by Jim Muth */
3414 {
3415 /* z=k*((a*(z^b))+(d*(z^f)))+c, */
3416 _CMPLX z_b, z_f;
3417 CMPLXpwr(old,B,z_b); /* (z^b) */
3418 CMPLXpwr(old,F,z_f); /* (z^f) */
3419 new.x = K.x*A.x*z_b.x + K.x*D.x*z_f.x + C.x;
3420 new.y = K.x*A.x*z_b.y + K.x*D.x*z_f.y + C.y;
3421 return(floatbailout());
3422 }
3423 #undef A
3424 #undef B
3425 #undef C
3426 #undef D
3427 #undef F
3428 #undef G
3429 #undef H
3430 #undef J
3431 #undef K
3432 #undef L
3433
3434 /*
3435 * The following functions calculate the real and imaginary complex
3436 * coordinates of the point in the complex plane corresponding to
3437 * the screen coordinates (col,row) at the current zoom corners
3438 * settings. The functions come in two flavors. One looks up the pixel
3439 * values using the precalculated grid arrays dx0, dx1, dy0, and dy1,
3440 * which has a speed advantage but is limited to MAXPIXELS image
3441 * dimensions. The other calculates the complex coordinates at a
3442 * cost of two additions and two multiplications for each component,
3443 * but works at any resolution.
3444 *
3445 * With Microsoft C's _fastcall keyword, the function call overhead
3446 * appears to be negligible. It also appears that the speed advantage
3447 * of the lookup vs the calculation is negligible on machines with
3448 * coprocessors. Bert Tyler's original implementation was designed for
3449 * machines with no coprocessor; on those machines the saving was
3450 * significant. For the time being, the table lookup capability will
3451 * be maintained.
3452 */
3453
3454 /* Real component, grid lookup version - requires dx0/dx1 arrays */
3455 static double _fastcall dxpixel_grid(void)
3456 {
3457 return(dx0[col]+dx1[row]);
3458 }
3459
3460 /* Real component, calculation version - does not require arrays */
3461 static double _fastcall dxpixel_calc(void)
3462 {
3463 return((double)(xxmin + col*delxx + row*delxx2));
3464 }
3465
3466 /* Imaginary component, grid lookup version - requires dy0/dy1 arrays */
3467 static double _fastcall dypixel_grid(void)
3468 {
3469 return(dy0[row]+dy1[col]);
3470 }
3471
3472 /* Imaginary component, calculation version - does not require arrays */
3473 static double _fastcall dypixel_calc(void)
3474 {
3475 return((double)(yymax - row*delyy - col*delyy2));
3476 }
3477
3478 /* Real component, grid lookup version - requires lx0/lx1 arrays */
3479 static long _fastcall lxpixel_grid(void)
3480 {
3481 return(lx0[col]+lx1[row]);
3482 }
3483
3484 /* Real component, calculation version - does not require arrays */
3485 static long _fastcall lxpixel_calc(void)
3486 {
3487 return(xmin + col*delx + row*delx2);
3488 }
3489
3490 /* Imaginary component, grid lookup version - requires ly0/ly1 arrays */
3491 static long _fastcall lypixel_grid(void)
3492 {
3493 return(ly0[row]+ly1[col]);
3494 }
3495
3496 /* Imaginary component, calculation version - does not require arrays */
3497 static long _fastcall lypixel_calc(void)
3498 {
3499 return(ymax - row*dely - col*dely2);
3500 }
3501
3502 double (_fastcall *dxpixel)(void) = dxpixel_calc;
3503 double (_fastcall *dypixel)(void) = dypixel_calc;
3504 long (_fastcall *lxpixel)(void) = lxpixel_calc;
3505 long (_fastcall *lypixel)(void) = lypixel_calc;
3506
3507 void set_pixel_calc_functions(void)
3508 {
3509 if(use_grid)
3510 {
3511 dxpixel = dxpixel_grid;
3512 dypixel = dypixel_grid;
3513 lxpixel = lxpixel_grid;
3514 lypixel = lypixel_grid;
3515 }
3516 else
3517 {
3518 dxpixel = dxpixel_calc;
3519 dypixel = dypixel_calc;
3520 lxpixel = lxpixel_calc;
3521 lypixel = lypixel_calc;
3522 }
3523 }
3524