home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / fractals.c < prev    next >
C/C++ Source or Header  |  1995-02-15  |  77KB  |  2,887 lines
  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 <stdio.h>
  41. #include <stdlib.h>
  42. #include <float.h>
  43. #include <limits.h>
  44. #include <string.h>
  45. #ifdef __TURBOC__
  46. #include <alloc.h>
  47. #elif !defined(__386BSD__)
  48. #include <malloc.h>
  49. #endif
  50. #include "fractint.h"
  51. #include "mpmath.h"
  52. #include "helpdefs.h"
  53. #include "fractype.h"
  54. #include "prototyp.h"
  55.  
  56. #define NEWTONDEGREELIMIT  100
  57.  
  58. #define CMPLXsqr_old(out)    \
  59.    (out).y = (old.x+old.x) * old.y;\
  60.    (out).x = tempsqrx - tempsqry
  61.  
  62. #define CMPLXpwr(arg1,arg2,out)   (out)= ComplexPower((arg1), (arg2))
  63. #define CMPLXmult1(arg1,arg2,out)    Arg2->d = (arg1); Arg1->d = (arg2);\
  64.      dStkMul(); Arg1++; Arg2++; (out) = Arg2->d
  65.  
  66. #define CMPLXmult(arg1,arg2,out)  \
  67.     {\
  68.        _CMPLX TmP;\
  69.        TmP.x = (arg1).x*(arg2).x - (arg1).y*(arg2).y;\
  70.        TmP.y = (arg1).x*(arg2).y + (arg1).y*(arg2).x;\
  71.        (out) = TmP;\
  72.      }
  73.  
  74. #define CMPLXadd(arg1,arg2,out)    \
  75.     (out).x = (arg1).x + (arg2).x; (out).y = (arg1).y + (arg2).y
  76. #define CMPLXsub(arg1,arg2,out)    \
  77.     (out).x = (arg1).x - (arg2).x; (out).y = (arg1).y - (arg2).y
  78. #define CMPLXtimesreal(arg,real,out)   \
  79.     (out).x = (arg).x*(real);\
  80.     (out).y = (arg).y*(real)
  81.  
  82. #define CMPLXrecip(arg,out)    \
  83.    { double denom; denom = sqr((arg).x) + sqr((arg).y);\
  84.      if(denom==0.0) {(out).x = 1.0e10;(out).y = 1.0e10;}else\
  85.     { (out).x =  (arg).x/denom;\
  86.      (out).y = -(arg).y/denom;}}
  87.  
  88. _LCMPLX lcoefficient,lold,lnew,lparm, linit,ltmp,ltmp2,lparm2;
  89. long ltempsqrx,ltempsqry;
  90. int maxcolor;
  91. int root, degree,basin;
  92. double floatmin,floatmax;
  93. double roverd, d1overd, threshold;
  94. _CMPLX tmp2;
  95. _CMPLX coefficient;
  96. _CMPLX    staticroots[16]; /* roots array for degree 16 or less */
  97. _CMPLX    *roots = staticroots;
  98. struct MPC    *MPCroots;
  99. long FgHalf;
  100. _CMPLX pwr;
  101. int    bitshiftless1;            /* bit shift less 1 */
  102.  
  103. #ifndef sqr
  104. #define sqr(x) ((x)*(x))
  105. #endif
  106.  
  107. #ifndef lsqr
  108. #define lsqr(x) (multiply((x),(x),bitshift))
  109. #endif
  110.  
  111. #define modulus(z)     (sqr((z).x)+sqr((z).y))
  112. #define conjugate(pz)    ((pz)->y = 0.0 - (pz)->y)
  113. #define distance(z1,z2)  (sqr((z1).x-(z2).x)+sqr((z1).y-(z2).y))
  114. #define pMPsqr(z) (*pMPmul((z),(z)))
  115. #define MPdistance(z1,z2)  (*pMPadd(pMPsqr(*pMPsub((z1).x,(z2).x)),pMPsqr(*pMPsub((z1).y,(z2).y))))
  116.  
  117. double twopi = PI*2.0;
  118. int c_exp;
  119.  
  120.  
  121. /* These are local but I don't want to pass them as parameters */
  122. _CMPLX parm,parm2;
  123. _CMPLX *floatparm;
  124. _LCMPLX *longparm; /* used here and in jb.c */
  125.  
  126. /* -------------------------------------------------------------------- */
  127. /*        These variables are external for speed's sake only      */
  128. /* -------------------------------------------------------------------- */
  129.  
  130. double sinx,cosx;
  131. double siny,cosy;
  132. double tmpexp;
  133. double tempsqrx,tempsqry;
  134.  
  135. double foldxinitx,foldyinity,foldxinity,foldyinitx;
  136. long oldxinitx,oldyinity,oldxinity,oldyinitx;
  137. long longtmp;
  138.  
  139. /* These are for quaternions */
  140. double qc,qci,qcj,qck;
  141.  
  142. /* temporary variables for trig use */
  143. long lcosx, lsinx;
  144. long lcosy, lsiny;
  145.  
  146. /*
  147. **  details of finite attractors (required for Magnet Fractals)
  148. **  (can also be used in "coloring in" the lakes of Julia types)
  149. */
  150.  
  151. /*
  152. **  pre-calculated values for fractal types Magnet2M & Magnet2J
  153. */
  154. _CMPLX    T_Cm1;          /* 3 * (floatparm - 1)            */
  155. _CMPLX    T_Cm2;          /* 3 * (floatparm - 2)            */
  156. _CMPLX    T_Cm1Cm2;     /* (floatparm - 1) * (floatparm - 2) */
  157.  
  158. void FloatPreCalcMagnet2() /* precalculation for Magnet2 (M & J) for speed */
  159.   {
  160.     T_Cm1.x = floatparm->x - 1.0;   T_Cm1.y = floatparm->y;
  161.     T_Cm2.x = floatparm->x - 2.0;   T_Cm2.y = floatparm->y;
  162.     T_Cm1Cm2.x = (T_Cm1.x * T_Cm2.x) - (T_Cm1.y * T_Cm2.y);
  163.     T_Cm1Cm2.y = (T_Cm1.x * T_Cm2.y) + (T_Cm1.y * T_Cm2.x);
  164.     T_Cm1.x += T_Cm1.x + T_Cm1.x;   T_Cm1.y += T_Cm1.y + T_Cm1.y;
  165.     T_Cm2.x += T_Cm2.x + T_Cm2.x;   T_Cm2.y += T_Cm2.y + T_Cm2.y;
  166.   }
  167.  
  168.  
  169. /* -------------------------------------------------------------------- */
  170. /*        Bailout Routines Macros                                                 */
  171. /* -------------------------------------------------------------------- */
  172.  
  173. int (near *floatbailout)(void);
  174. int (near *longbailout)(void);
  175. int (near *bignumbailout)(void);
  176. int (near *bigfltbailout)(void);
  177.  
  178. int near fpMODbailout()
  179. {
  180.    if ( ( magnitude = ( tempsqrx=sqr(new.x) )
  181.             + ( tempsqry=sqr(new.y) ) ) >= rqlim ) return(1);
  182.    old = new;
  183.    return(0);
  184. }
  185.  
  186. int near fpREALbailout()
  187. {
  188.    tempsqrx=sqr(new.x);
  189.    tempsqry=sqr(new.y);
  190.    magnitude = tempsqrx + tempsqry;
  191.    if(tempsqrx >= rqlim) return(1);
  192.    old = new;
  193.    return(0);
  194. }
  195.  
  196. int near fpIMAGbailout()
  197. {
  198.    tempsqrx=sqr(new.x);
  199.    tempsqry=sqr(new.y);
  200.    magnitude = tempsqrx + tempsqry;
  201.    if(tempsqry >= rqlim) return(1);
  202.    old = new;
  203.    return(0);
  204. }
  205.  
  206. int near fpORbailout()
  207. {
  208.    tempsqrx=sqr(new.x);
  209.    tempsqry=sqr(new.y);
  210.    magnitude = tempsqrx + tempsqry;
  211.    if(tempsqrx >= rqlim || tempsqry >= rqlim) return(1);
  212.    old = new;
  213.    return(0);
  214. }
  215.  
  216. int near fpANDbailout()
  217. {
  218.    tempsqrx=sqr(new.x);
  219.    tempsqry=sqr(new.y);
  220.    magnitude = tempsqrx + tempsqry;
  221.    if(tempsqrx >= rqlim && tempsqry >= rqlim) return(1);
  222.    old = new;
  223.    return(0);
  224. }
  225.  
  226. /* longbailout() is equivalent to next */
  227. #define LONGBAILOUT()    \
  228.    ltempsqrx = lsqr(lnew.x); ltempsqry = lsqr(lnew.y);\
  229.    lmagnitud = ltempsqrx + ltempsqry;\
  230.    if (lmagnitud >= llimit || lmagnitud < 0 || labs(lnew.x) > llimit2\
  231.      || labs(lnew.y) > llimit2) \
  232.            { return(1);}\
  233.    lold = lnew;
  234.  
  235. #define FLOATTRIGBAILOUT()  \
  236.    if (fabs(old.y) >= rqlim2) return(1);
  237.  
  238. #define LONGTRIGBAILOUT()  \
  239.    if(labs(lold.y) >= llimit2) { return(1);}
  240.  
  241. #define LONGXYTRIGBAILOUT()  \
  242.    if(labs(lold.x) >= llimit2 || labs(lold.y) >= llimit2)\
  243.     { return(1);}
  244.  
  245. #define FLOATXYTRIGBAILOUT()  \
  246.    if (fabs(old.x) >= rqlim2 || fabs(old.y) >= rqlim2) return(1);
  247.  
  248. #define FLOATHTRIGBAILOUT()  \
  249.    if (fabs(old.x) >= rqlim2) return(1);
  250.  
  251. #define LONGHTRIGBAILOUT()  \
  252.    if(labs(lold.x) >= llimit2) { return(1);}
  253.  
  254. #define TRIG16CHECK(X)    \
  255.       if(labs((X)) > l16triglim) { return(1);}
  256.  
  257. #define FLOATEXPBAILOUT()  \
  258.    if (fabs(old.y) >= 1.0e8) return(1);\
  259.    if (fabs(old.x) >= 6.4e2) return(1);
  260.  
  261. #define LONGEXPBAILOUT()  \
  262.    if (labs(lold.y) >= (1000L<<bitshift)) return(1);\
  263.    if (labs(lold.x) >=      (8L<<bitshift)) return(1);
  264.  
  265. #if 0
  266. /* this define uses usual trig instead of fast trig */
  267. #define FPUsincos(px,psinx,pcosx) \
  268.    *(psinx) = sin(*(px));\
  269.    *(pcosx) = cos(*(px));
  270.  
  271. #define FPUsinhcosh(px,psinhx,pcoshx) \
  272.    *(psinhx) = sinh(*(px));\
  273.    *(pcoshx) = cosh(*(px));
  274. #endif
  275.  
  276. #define LTRIGARG(X)    \
  277.    if(labs((X)) > l16triglim)\
  278.    {\
  279.       double tmp;\
  280.       tmp = (X);\
  281.       tmp /= fudge;\
  282.       tmp = fmod(tmp,twopi);\
  283.       tmp *= fudge;\
  284.       (X) = (long)tmp;\
  285.    }\
  286.  
  287. static int near Halleybailout(void)
  288. {
  289.    if ( fabs(modulus(new)-modulus(old)) < parm2.x)
  290.       return(1);
  291.    old = new;
  292.    return(0);
  293. }
  294.  
  295. #ifndef XFRACT
  296. #define MPCmod(m) (*pMPadd(*pMPmul((m).x, (m).x), *pMPmul((m).y, (m).y)))
  297. struct MPC mpcold, mpcnew, mpctmp, mpctmp1;
  298. struct MP mptmpparm2x;
  299.  
  300. static int near MPCHalleybailout(void)
  301. {
  302. static struct MP mptmpbailout;
  303.    mptmpbailout = *MPabs(*pMPsub(MPCmod(mpcnew), MPCmod(mpcold)));
  304.    if (pMPcmp(mptmpbailout, mptmpparm2x) < 0)
  305.       return(1);
  306.    mpcold = mpcnew;
  307.    return(0);
  308. }
  309. #endif
  310.  
  311. #ifdef XFRACT
  312. int asmlMODbailout() { return 0;}
  313. int asmlREALbailout() { return 0;}
  314. int asmlIMAGbailout() { return 0;}
  315. int asmlORbailout() { return 0;}
  316. int asmlANDbailout() { return 0;}
  317. int asm386lMODbailout() { return 0;}
  318. int asm386lREALbailout() { return 0;}
  319. int asm386lIMAGbailout() { return 0;}
  320. int asm386lORbailout() { return 0;}
  321. int asm386lANDbailout() { return 0;}
  322. int asmfpMODbailout() { return 0;}
  323. int asmfpREALbailout() { return 0;}
  324. int asmfpIMAGbailout() { return 0;}
  325. int asmfpORbailout() { return 0;}
  326. int asmfpANDbailout() { return 0;}
  327. #endif
  328.  
  329. /* -------------------------------------------------------------------- */
  330. /*        Fractal (once per iteration) routines            */
  331. /* -------------------------------------------------------------------- */
  332. static double xt, yt, t2;
  333.  
  334. /* Raise complex number (base) to the (exp) power, storing the result
  335. ** in complex (result).
  336. */
  337. void cpower(_CMPLX *base, int exp, _CMPLX *result)
  338. {
  339.     if (exp<0) {
  340.     cpower(base,-exp,result);
  341.     CMPLXrecip(*result,*result);
  342.     return;
  343.     }
  344.  
  345.     xt = base->x;   yt = base->y;
  346.  
  347.     if (exp & 1)
  348.     {
  349.        result->x = xt;
  350.        result->y = yt;
  351.     }
  352.     else
  353.     {
  354.        result->x = 1.0;
  355.        result->y = 0.0;
  356.     }
  357.  
  358.     exp >>= 1;
  359.     while (exp)
  360.     {
  361.     t2 = xt * xt - yt * yt;
  362.     yt = 2 * xt * yt;
  363.     xt = t2;
  364.  
  365.     if (exp & 1)
  366.     {
  367.         t2 = xt * result->x - yt * result->y;
  368.         result->y = result->y * xt + yt * result->x;
  369.         result->x = t2;
  370.     }
  371.     exp >>= 1;
  372.     }
  373. }
  374. /* long version */
  375. static long lxt, lyt, lt2;
  376.  
  377. lcpower(_LCMPLX *base, int exp, _LCMPLX *result, int bitshift)
  378. {
  379.     static long maxarg;
  380.     maxarg = 64L<<bitshift;
  381.  
  382.     if (exp<0) {
  383.     overflow = lcpower(base,-exp,result,bitshift);
  384.     LCMPLXrecip(*result,*result);
  385.     return(overflow);
  386.     }
  387.  
  388.     overflow = 0;
  389.     lxt = base->x;   lyt = base->y;
  390.  
  391.     if (exp & 1)
  392.     {
  393.        result->x = lxt;
  394.        result->y = lyt;
  395.     }
  396.     else
  397.     {
  398.        result->x = 1L<<bitshift;
  399.        result->y = 0L;
  400.     }
  401.  
  402.     exp >>= 1;
  403.     while (exp)
  404.     {
  405.     /*
  406.     if(labs(lxt) >= maxarg || labs(lyt) >= maxarg)
  407.        return(-1);
  408.     */
  409.     lt2 = multiply(lxt, lxt, bitshift) - multiply(lyt,lyt,bitshift);
  410.     lyt = multiply(lxt,lyt,bitshiftless1);
  411.     if(overflow)
  412.        return(overflow);
  413.     lxt = lt2;
  414.  
  415.     if (exp & 1)
  416.     {
  417.         lt2 = multiply(lxt,result->x, bitshift) - multiply(lyt,result->y,bitshift);
  418.         result->y = multiply(result->y,lxt,bitshift) + multiply(lyt,result->x,bitshift);
  419.         result->x = lt2;
  420.     }
  421.     exp >>= 1;
  422.     }
  423.     if(result->x == 0 && result->y == 0)
  424.        overflow = 1;
  425.     return(overflow);
  426. }
  427. #if 0
  428. z_to_the_z(_CMPLX *z, _CMPLX *out)
  429. {
  430.     static _CMPLX tmp1,tmp2;
  431.     /* raises complex z to the z power */
  432.     int errno_xxx;
  433.     errno_xxx = 0;
  434.  
  435.     if(fabs(z->x) < DBL_EPSILON) return(-1);
  436.  
  437.     /* log(x + iy) = 1/2(log(x*x + y*y) + i(arc_tan(y/x)) */
  438.     tmp1.x = .5*log(sqr(z->x)+sqr(z->y));
  439.  
  440.     /* the fabs in next line added to prevent discontinuity in image */
  441.     tmp1.y = atan(fabs(z->y/z->x));
  442.  
  443.     /* log(z)*z */
  444.     tmp2.x = tmp1.x * z->x - tmp1.y * z->y;
  445.     tmp2.y = tmp1.x * z->y + tmp1.y * z->x;
  446.  
  447.     /* z*z = e**(log(z)*z) */
  448.     /* e**(x + iy) =  e**x * (cos(y) + isin(y)) */
  449.  
  450.     tmpexp = exp(tmp2.x);
  451.  
  452.     FPUsincos(&tmp2.y,&siny,&cosy);
  453.     out->x = tmpexp*cosy;
  454.     out->y = tmpexp*siny;
  455.     return(errno_xxx);
  456. }
  457. #endif
  458.  
  459. #ifdef XFRACT /* fractint uses the NewtonFractal2 code in newton.asm */
  460.  
  461. int complex_div(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
  462. int complex_mult(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
  463.  
  464. /* Distance of complex z from unit circle */
  465. #define DIST1(z) (((z).x-1.0)*((z).x-1.0)+((z).y)*((z).y))
  466. #define LDIST1(z) (lsqr((((z).x)-fudge)) + lsqr(((z).y)))
  467.  
  468.  
  469. int NewtonFractal2()
  470. {
  471.     static char start=1;
  472.     if(start)
  473.     {
  474.        start = 0;
  475.     }
  476.     cpower(&old, degree-1, &tmp);
  477.     complex_mult(tmp, old, &new);
  478.  
  479.     if (DIST1(new) < threshold)
  480.     {
  481.        if(fractype==NEWTBASIN || fractype==MPNEWTBASIN)
  482.        {
  483.       long tmpcolor;
  484.       int i;
  485.       tmpcolor = -1;
  486.       /* this code determines which degree-th root of root the
  487.          Newton formula converges to. The roots of a 1 are
  488.          distributed on a circle of radius 1 about the origin. */
  489.       for(i=0;i<degree;i++)
  490.          /* color in alternating shades with iteration according to
  491.         which root of 1 it converged to */
  492.           if(distance(roots[i],old) < threshold)
  493.           {
  494.           if (basin==2) {
  495.                   tmpcolor = 1+(i&7)+((coloriter&1)<<3);
  496.           } else {
  497.               tmpcolor = 1+i;
  498.           }
  499.           break;
  500.           }
  501.        if(tmpcolor == -1)
  502.           coloriter = maxcolor;
  503.        else
  504.           coloriter = tmpcolor;
  505.        }
  506.        return(1);
  507.     }
  508.     new.x = d1overd * new.x + roverd;
  509.     new.y *= d1overd;
  510.  
  511.     /* Watch for divide underflow */
  512.     if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN)
  513.       return(1);
  514.     else
  515.     {
  516.     t2 = 1.0 / t2;
  517.     old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
  518.     old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
  519.     }
  520.     return(0);
  521. }
  522.  
  523. complex_mult(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz)
  524. {
  525.    pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
  526.    pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
  527.    return(0);
  528. }
  529.  
  530.  
  531. complex_div(_CMPLX numerator,_CMPLX denominator,_CMPLX *pout)
  532. {
  533.    double mod;
  534.    if((mod = modulus(denominator)) < FLT_MIN)
  535.       return(1);
  536.    conjugate(&denominator);
  537.    complex_mult(numerator,denominator,pout);
  538.    pout->x = pout->x/mod;
  539.    pout->y = pout->y/mod;
  540.    return(0);
  541. }
  542. #endif /* newton code only used by xfractint */
  543.  
  544. #ifndef XFRACT
  545. struct MP mproverd, mpd1overd, mpthreshold;
  546. struct MP mpt2;
  547. struct MP mpone;
  548. #endif
  549.  
  550. int MPCNewtonFractal()
  551. {
  552. #ifndef XFRACT
  553.     MPOverflow = 0;
  554.     mpctmp   = MPCpow(mpcold,degree-1);
  555.  
  556.     mpcnew.x = *pMPsub(*pMPmul(mpctmp.x,mpcold.x),*pMPmul(mpctmp.y,mpcold.y));
  557.     mpcnew.y = *pMPadd(*pMPmul(mpctmp.x,mpcold.y),*pMPmul(mpctmp.y,mpcold.x));
  558.     mpctmp1.x = *pMPsub(mpcnew.x, MPCone.x);
  559.     mpctmp1.y = *pMPsub(mpcnew.y, MPCone.y);
  560.     if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
  561.     {
  562.       if(fractype==MPNEWTBASIN)
  563.       {
  564.      long tmpcolor;
  565.      int i;
  566.      tmpcolor = -1;
  567.      for(i=0;i<degree;i++)
  568.          if(pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
  569.          {
  570.         if(basin==2)
  571.            tmpcolor = 1+(i&7) + ((coloriter&1)<<3);
  572.         else
  573.            tmpcolor = 1+i;
  574.             break;
  575.          }
  576.       if(tmpcolor == -1)
  577.          coloriter = maxcolor;
  578.       else
  579.          coloriter = tmpcolor;
  580.       }
  581.        return(1);
  582.     }
  583.  
  584.     mpcnew.x = *pMPadd(*pMPmul(mpd1overd,mpcnew.x),mproverd);
  585.     mpcnew.y = *pMPmul(mpcnew.y,mpd1overd);
  586.     mpt2 = MPCmod(mpctmp);
  587.     mpt2 = *pMPdiv(mpone,mpt2);
  588.     mpcold.x = *pMPmul(mpt2,(*pMPadd(*pMPmul(mpcnew.x,mpctmp.x),*pMPmul(mpcnew.y,mpctmp.y))));
  589.     mpcold.y = *pMPmul(mpt2,(*pMPsub(*pMPmul(mpcnew.y,mpctmp.x),*pMPmul(mpcnew.x,mpctmp.y))));
  590.     new.x = *pMP2d(mpcold.x);
  591.     new.y = *pMP2d(mpcold.y);
  592.     return(MPOverflow);
  593. #endif
  594. }
  595.  
  596.  
  597. Barnsley1Fractal()
  598. {
  599. #ifndef XFRACT
  600.    /* Barnsley's Mandelbrot type M1 from "Fractals
  601.    Everywhere" by Michael Barnsley, p. 322 */
  602.  
  603.    /* calculate intermediate products */
  604.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  605.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  606.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  607.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  608.    /* orbit calculation */
  609.    if(lold.x >= 0)
  610.    {
  611.       lnew.x = (oldxinitx - longparm->x - oldyinity);
  612.       lnew.y = (oldyinitx - longparm->y + oldxinity);
  613.    }
  614.    else
  615.    {
  616.       lnew.x = (oldxinitx + longparm->x - oldyinity);
  617.       lnew.y = (oldyinitx + longparm->y + oldxinity);
  618.    }
  619.    return(longbailout());
  620. #endif
  621. }
  622.  
  623. Barnsley1FPFractal()
  624. {
  625.    /* Barnsley's Mandelbrot type M1 from "Fractals
  626.    Everywhere" by Michael Barnsley, p. 322 */
  627.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  628.  
  629.    /* calculate intermediate products */
  630.    foldxinitx = old.x * floatparm->x;
  631.    foldyinity = old.y * floatparm->y;
  632.    foldxinity = old.x * floatparm->y;
  633.    foldyinitx = old.y * floatparm->x;
  634.    /* orbit calculation */
  635.    if(old.x >= 0)
  636.    {
  637.       new.x = (foldxinitx - floatparm->x - foldyinity);
  638.       new.y = (foldyinitx - floatparm->y + foldxinity);
  639.    }
  640.    else
  641.    {
  642.       new.x = (foldxinitx + floatparm->x - foldyinity);
  643.       new.y = (foldyinitx + floatparm->y + foldxinity);
  644.    }
  645.    return(floatbailout());
  646. }
  647.  
  648. Barnsley2Fractal()
  649. {
  650. #ifndef XFRACT
  651.    /* An unnamed Mandelbrot/Julia function from "Fractals
  652.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  653.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  654.  
  655.    /* calculate intermediate products */
  656.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  657.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  658.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  659.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  660.  
  661.    /* orbit calculation */
  662.    if(oldxinity + oldyinitx >= 0)
  663.    {
  664.       lnew.x = oldxinitx - longparm->x - oldyinity;
  665.       lnew.y = oldyinitx - longparm->y + oldxinity;
  666.    }
  667.    else
  668.    {
  669.       lnew.x = oldxinitx + longparm->x - oldyinity;
  670.       lnew.y = oldyinitx + longparm->y + oldxinity;
  671.    }
  672.    return(longbailout());
  673. #endif
  674. }
  675.  
  676. Barnsley2FPFractal()
  677. {
  678.    /* An unnamed Mandelbrot/Julia function from "Fractals
  679.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  680.  
  681.    /* calculate intermediate products */
  682.    foldxinitx = old.x * floatparm->x;
  683.    foldyinity = old.y * floatparm->y;
  684.    foldxinity = old.x * floatparm->y;
  685.    foldyinitx = old.y * floatparm->x;
  686.  
  687.    /* orbit calculation */
  688.    if(foldxinity + foldyinitx >= 0)
  689.    {
  690.       new.x = foldxinitx - floatparm->x - foldyinity;
  691.       new.y = foldyinitx - floatparm->y + foldxinity;
  692.    }
  693.    else
  694.    {
  695.       new.x = foldxinitx + floatparm->x - foldyinity;
  696.       new.y = foldyinitx + floatparm->y + foldxinity;
  697.    }
  698.    return(floatbailout());
  699. }
  700.  
  701. JuliaFractal()
  702. {
  703. #ifndef XFRACT
  704.    /* used for C prototype of fast integer math routines for classic
  705.       Mandelbrot and Julia */
  706.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  707.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  708.    return(longbailout());
  709. #elif !defined(__386BSD__)
  710.    fprintf(stderr,"JuliaFractal called\n");
  711.    exit(-1);
  712. #endif
  713. }
  714.  
  715. JuliafpFractal()
  716. {
  717.    /* floating point version of classical Mandelbrot/Julia */
  718.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  719.    new.x = tempsqrx - tempsqry + floatparm->x;
  720.    new.y = 2.0 * old.x * old.y + floatparm->y;
  721.    return(floatbailout());
  722. }
  723.  
  724. LambdaFPFractal()
  725. {
  726.    /* variation of classical Mandelbrot/Julia */
  727.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  728.  
  729.    tempsqrx = old.x - tempsqrx + tempsqry;
  730.    tempsqry = -(old.y * old.x);
  731.    tempsqry += tempsqry + old.y;
  732.  
  733.    new.x = floatparm->x * tempsqrx - floatparm->y * tempsqry;
  734.    new.y = floatparm->x * tempsqry + floatparm->y * tempsqrx;
  735.    return(floatbailout());
  736. }
  737.  
  738. LambdaFractal()
  739. {
  740. #ifndef XFRACT
  741.    /* variation of classical Mandelbrot/Julia */
  742.  
  743.    /* in complex math) temp = Z * (1-Z) */
  744.    ltempsqrx = lold.x - ltempsqrx + ltempsqry;
  745.    ltempsqry = lold.y
  746.          - multiply(lold.y, lold.x, bitshiftless1);
  747.    /* (in complex math) Z = Lambda * Z */
  748.    lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
  749.     - multiply(longparm->y, ltempsqry, bitshift);
  750.    lnew.y = multiply(longparm->x, ltempsqry, bitshift)
  751.     + multiply(longparm->y, ltempsqrx, bitshift);
  752.    return(longbailout());
  753. #endif
  754. }
  755.  
  756. SierpinskiFractal()
  757. {
  758. #ifndef XFRACT
  759.    /* following code translated from basic - see "Fractals
  760.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  761.    lnew.x = (lold.x << 1);        /* new.x = 2 * old.x  */
  762.    lnew.y = (lold.y << 1);        /* new.y = 2 * old.y  */
  763.    if(lold.y > ltmp.y)    /* if old.y > .5 */
  764.       lnew.y = lnew.y - ltmp.x; /* new.y = 2 * old.y - 1 */
  765.    else if(lold.x > ltmp.y)    /* if old.x > .5 */
  766.       lnew.x = lnew.x - ltmp.x; /* new.x = 2 * old.x - 1 */
  767.    /* end barnsley code */
  768.    return(longbailout());
  769. #endif
  770. }
  771.  
  772. SierpinskiFPFractal()
  773. {
  774.    /* following code translated from basic - see "Fractals
  775.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  776.  
  777.    new.x = old.x + old.x;
  778.    new.y = old.y + old.y;
  779.    if(old.y > .5)
  780.       new.y = new.y - 1;
  781.    else if (old.x > .5)
  782.       new.x = new.x - 1;
  783.  
  784.    /* end barnsley code */
  785.    return(floatbailout());
  786. }
  787.  
  788. LambdaexponentFractal()
  789. {
  790.    /* found this in  "Science of Fractal Images" */
  791.    FLOATEXPBAILOUT();
  792.    FPUsincos  (&old.y,&siny,&cosy);
  793.  
  794.    if (old.x >= rqlim && cosy >= 0.0) return(1);
  795.    tmpexp = exp(old.x);
  796.    tmp.x = tmpexp*cosy;
  797.    tmp.y = tmpexp*siny;
  798.  
  799.    /*multiply by lamda */
  800.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  801.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  802.    old = new;
  803.    return(0);
  804. }
  805.  
  806. LongLambdaexponentFractal()
  807. {
  808. #ifndef XFRACT
  809.    /* found this in  "Science of Fractal Images" */
  810.    LONGEXPBAILOUT();
  811.  
  812.    SinCos086  (lold.y, &lsiny,    &lcosy);
  813.  
  814.    if (lold.x >= llimit && lcosy >= 0L) return(1);
  815.    longtmp = Exp086(lold.x);
  816.  
  817.    ltmp.x = multiply(longtmp,       lcosy,   bitshift);
  818.    ltmp.y = multiply(longtmp,       lsiny,   bitshift);
  819.  
  820.    lnew.x  = multiply(longparm->x, ltmp.x, bitshift)
  821.        - multiply(longparm->y, ltmp.y, bitshift);
  822.    lnew.y  = multiply(longparm->x, ltmp.y, bitshift)
  823.        + multiply(longparm->y, ltmp.x, bitshift);
  824.    lold = lnew;
  825.    return(0);
  826. #endif
  827. }
  828.  
  829. FloatTrigPlusExponentFractal()
  830. {
  831.    /* another Scientific American biomorph type */
  832.    /* z(n+1) = e**z(n) + trig(z(n)) + C */
  833.  
  834.    if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
  835.    tmpexp = exp(old.x);
  836.    FPUsincos  (&old.y,&siny,&cosy);
  837.    CMPLXtrig0(old,new);
  838.  
  839.    /*new =   trig(old) + e**old + C  */
  840.    new.x += tmpexp*cosy + floatparm->x;
  841.    new.y += tmpexp*siny + floatparm->y;
  842.    return(floatbailout());
  843. }
  844.  
  845.  
  846. LongTrigPlusExponentFractal()
  847. {
  848. #ifndef XFRACT
  849.    /* calculate exp(z) */
  850.  
  851.    /* domain check for fast transcendental functions */
  852.    TRIG16CHECK(lold.x);
  853.    TRIG16CHECK(lold.y);
  854.  
  855.    longtmp = Exp086(lold.x);
  856.    SinCos086  (lold.y, &lsiny,    &lcosy);
  857.    LCMPLXtrig0(lold,lnew);
  858.    lnew.x += multiply(longtmp,      lcosy,   bitshift) + longparm->x;
  859.    lnew.y += multiply(longtmp,      lsiny,   bitshift) + longparm->y;
  860.    return(longbailout());
  861. #endif
  862. }
  863.  
  864. MarksLambdaFractal()
  865. {
  866.    /* Mark Peterson's variation of "lambda" function */
  867.  
  868.    /* Z1 = (C^(exp-1) * Z**2) + C */
  869.    ltmp.x = ltempsqrx - ltempsqry;
  870.    ltmp.y = multiply(lold.x ,lold.y ,bitshiftless1);
  871.  
  872.    lnew.x = multiply(lcoefficient.x, ltmp.x, bitshift)
  873.     - multiply(lcoefficient.y, ltmp.y, bitshift) + longparm->x;
  874.    lnew.y = multiply(lcoefficient.x, ltmp.y, bitshift)
  875.     + multiply(lcoefficient.y, ltmp.x, bitshift) + longparm->y;
  876.  
  877.    return(longbailout());
  878. }
  879.  
  880. MarksLambdafpFractal()
  881. {
  882.    /* Mark Peterson's variation of "lambda" function */
  883.  
  884.    /* Z1 = (C^(exp-1) * Z**2) + C */
  885.    tmp.x = tempsqrx - tempsqry;
  886.    tmp.y = old.x * old.y *2;
  887.  
  888.    new.x = coefficient.x * tmp.x - coefficient.y * tmp.y + floatparm->x;
  889.    new.y = coefficient.x * tmp.y + coefficient.y * tmp.x + floatparm->y;
  890.  
  891.    return(floatbailout());
  892. }
  893.  
  894.  
  895. long XXOne, FgOne, FgTwo;
  896.  
  897. UnityFractal()
  898. {
  899. #ifndef XFRACT
  900.    /* brought to you by Mark Peterson - you won't find this in any fractal
  901.       books unless they saw it here first - Mark invented it! */
  902.    XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
  903.    if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delmin))
  904.       return(1);
  905.    lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
  906.    lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
  907.    lnew=lold;  /* TW added this line */
  908.    return(0);
  909. #endif
  910. }
  911.  
  912. #define XXOne new.x
  913.  
  914. UnityfpFractal()
  915. {
  916.    /* brought to you by Mark Peterson - you won't find this in any fractal
  917.       books unless they saw it here first - Mark invented it! */
  918.  
  919.    XXOne = sqr(old.x) + sqr(old.y);
  920.    if((XXOne > 2.0) || (fabs(XXOne - 1.0) < ddelmin))
  921.       return(1);
  922.    old.y = (2.0 - XXOne)* old.x;
  923.    old.x = (2.0 - XXOne)* old.y;
  924.    new=old;  /* TW added this line */
  925.    return(0);
  926. }
  927.  
  928. #undef XXOne
  929.  
  930. Mandel4Fractal()
  931. {
  932.    /* By writing this code, Bert has left behind the excuse "don't
  933.       know what a fractal is, just know how to make'em go fast".
  934.       Bert is hereby declared a bonafide fractal expert! Supposedly
  935.       this routine calculates the Mandelbrot/Julia set based on the
  936.       polynomial z**4 + lambda, but I wouldn't know -- can't follow
  937.       all that integer math speedup stuff - Tim */
  938.  
  939.    /* first, compute (x + iy)**2 */
  940.    lnew.x  = ltempsqrx - ltempsqry;
  941.    lnew.y = multiply(lold.x, lold.y, bitshiftless1);
  942.    if (longbailout()) return(1);
  943.  
  944.    /* then, compute ((x + iy)**2)**2 + lambda */
  945.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  946.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  947.    return(longbailout());
  948. }
  949.  
  950. Mandel4fpFractal()
  951. {
  952.    /* first, compute (x + iy)**2 */
  953.    new.x  = tempsqrx - tempsqry;
  954.    new.y = old.x*old.y*2;
  955.    if (floatbailout()) return(1);
  956.  
  957.    /* then, compute ((x + iy)**2)**2 + lambda */
  958.    new.x  = tempsqrx - tempsqry + floatparm->x;
  959.    new.y =  old.x*old.y*2 + floatparm->y;
  960.    return(floatbailout());
  961. }
  962.  
  963. floatZtozPluszpwrFractal()
  964. {
  965.    cpower(&old,(int)param[2],&new);
  966.    old = ComplexPower(old,old);
  967.    new.x = new.x + old.x +floatparm->x;
  968.    new.y = new.y + old.y +floatparm->y;
  969.    return(floatbailout());
  970. }
  971.  
  972. longZpowerFractal()
  973. {
  974. #ifndef XFRACT
  975.    if(lcpower(&lold,c_exp,&lnew,bitshift))
  976.       lnew.x = lnew.y = 8L<<bitshift;
  977.    lnew.x += longparm->x;
  978.    lnew.y += longparm->y;
  979.    return(longbailout());
  980. #endif
  981. }
  982.  
  983. longCmplxZpowerFractal()
  984. {
  985. #ifndef XFRACT
  986.    _CMPLX x, y;
  987.  
  988.    x.x = (double)lold.x / fudge;
  989.    x.y = (double)lold.y / fudge;
  990.    y.x = (double)lparm2.x / fudge;
  991.    y.y = (double)lparm2.y / fudge;
  992.    x = ComplexPower(x, y);
  993.    if(fabs(x.x) < fgLimit && fabs(x.y) < fgLimit) {
  994.       lnew.x = (long)(x.x * fudge);
  995.       lnew.y = (long)(x.y * fudge);
  996.    }
  997.    else
  998.       overflow = 1;
  999.    lnew.x += longparm->x;
  1000.    lnew.y += longparm->y;
  1001.    return(longbailout());
  1002. #endif
  1003. }
  1004.  
  1005. floatZpowerFractal()
  1006. {
  1007.    cpower(&old,c_exp,&new);
  1008.    new.x += floatparm->x;
  1009.    new.y += floatparm->y;
  1010.    return(floatbailout());
  1011. }
  1012.  
  1013. floatCmplxZpowerFractal()
  1014. {
  1015.    new = ComplexPower(old, parm2);
  1016.    new.x += floatparm->x;
  1017.    new.y += floatparm->y;
  1018.    return(floatbailout());
  1019. }
  1020.  
  1021. Barnsley3Fractal()
  1022. {
  1023.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1024.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  1025.  
  1026.    /* calculate intermediate products */
  1027.    oldxinitx   = multiply(lold.x, lold.x, bitshift);
  1028.    oldyinity   = multiply(lold.y, lold.y, bitshift);
  1029.    oldxinity   = multiply(lold.x, lold.y, bitshift);
  1030.  
  1031.    /* orbit calculation */
  1032.    if(lold.x > 0)
  1033.    {
  1034.       lnew.x = oldxinitx   - oldyinity - fudge;
  1035.       lnew.y = oldxinity << 1;
  1036.    }
  1037.    else
  1038.    {
  1039.       lnew.x = oldxinitx - oldyinity - fudge
  1040.        + multiply(longparm->x,lold.x,bitshift);
  1041.       lnew.y = oldxinity <<1;
  1042.  
  1043.       /* This term added by Tim Wegner to make dependent on the
  1044.      imaginary part of the parameter. (Otherwise Mandelbrot
  1045.      is uninteresting. */
  1046.       lnew.y += multiply(longparm->y,lold.x,bitshift);
  1047.    }
  1048.    return(longbailout());
  1049. }
  1050.  
  1051. Barnsley3FPFractal()
  1052. {
  1053.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1054.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  1055.  
  1056.  
  1057.    /* calculate intermediate products */
  1058.    foldxinitx  = old.x * old.x;
  1059.    foldyinity  = old.y * old.y;
  1060.    foldxinity  = old.x * old.y;
  1061.  
  1062.    /* orbit calculation */
  1063.    if(old.x > 0)
  1064.    {
  1065.       new.x = foldxinitx - foldyinity - 1.0;
  1066.       new.y = foldxinity * 2;
  1067.    }
  1068.    else
  1069.    {
  1070.       new.x = foldxinitx - foldyinity -1.0 + floatparm->x * old.x;
  1071.       new.y = foldxinity * 2;
  1072.  
  1073.       /* This term added by Tim Wegner to make dependent on the
  1074.      imaginary part of the parameter. (Otherwise Mandelbrot
  1075.      is uninteresting. */
  1076.       new.y += floatparm->y * old.x;
  1077.    }
  1078.    return(floatbailout());
  1079. }
  1080.  
  1081. TrigPlusZsquaredFractal()
  1082. {
  1083. #ifndef XFRACT
  1084.    /* From Scientific American, July 1989 */
  1085.    /* A Biomorph              */
  1086.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1087.    LCMPLXtrig0(lold,lnew);
  1088.    lnew.x += ltempsqrx - ltempsqry + longparm->x;
  1089.    lnew.y += multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  1090.    return(longbailout());
  1091. #endif
  1092. }
  1093.  
  1094. TrigPlusZsquaredfpFractal()
  1095. {
  1096.    /* From Scientific American, July 1989 */
  1097.    /* A Biomorph              */
  1098.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1099.  
  1100.    CMPLXtrig0(old,new);
  1101.    new.x += tempsqrx - tempsqry + floatparm->x;
  1102.    new.y += 2.0 * old.x * old.y + floatparm->y;
  1103.    return(floatbailout());
  1104. }
  1105.  
  1106. Richard8fpFractal()
  1107. {
  1108.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1109.    CMPLXtrig0(old,new);
  1110. /*   CMPLXtrig1(*floatparm,tmp); */
  1111.    new.x += tmp.x;
  1112.    new.y += tmp.y;
  1113.    return(floatbailout());
  1114. }
  1115.  
  1116. Richard8Fractal()
  1117. {
  1118. #ifndef XFRACT
  1119.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1120.    LCMPLXtrig0(lold,lnew);
  1121. /*   LCMPLXtrig1(*longparm,ltmp); */
  1122.    lnew.x += ltmp.x;
  1123.    lnew.y += ltmp.y;
  1124.    return(longbailout());
  1125. #endif
  1126. }
  1127.  
  1128. PopcornFractal()
  1129. {
  1130.    tmp = old;
  1131.    tmp.x *= 3.0;
  1132.    tmp.y *= 3.0;
  1133.    FPUsincos(&tmp.x,&sinx,&cosx);
  1134.    FPUsincos(&tmp.y,&siny,&cosy);
  1135.    tmp.x = sinx/cosx + old.x;
  1136.    tmp.y = siny/cosy + old.y;
  1137.    FPUsincos(&tmp.x,&sinx,&cosx);
  1138.    FPUsincos(&tmp.y,&siny,&cosy);
  1139.    new.x = old.x - parm.x*siny;
  1140.    new.y = old.y - parm.x*sinx;
  1141.    /*
  1142.    new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
  1143.    new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
  1144.    */
  1145.    if(plot == noplot)
  1146.    {
  1147.       plot_orbit(new.x,new.y,1+row%colors);
  1148.       old = new;
  1149.    }
  1150.    else
  1151.    /* FLOATBAILOUT(); */
  1152.    /* PB The above line was weird, not what it seems to be!  But, bracketing
  1153.      it or always doing it (either of which seem more likely to be what
  1154.      was intended) changes the image for the worse, so I'm not touching it.
  1155.      Same applies to int form in next routine. */
  1156.    /* PB later: recoded inline, still leaving it weird */
  1157.       tempsqrx = sqr(new.x);
  1158.    tempsqry = sqr(new.y);
  1159.    if((magnitude = tempsqrx + tempsqry) >= rqlim) return(1);
  1160.    old = new;
  1161.    return(0);
  1162. }
  1163.  
  1164. LPopcornFractal()
  1165. {
  1166. #ifndef XFRACT
  1167.    ltmp = lold;
  1168.    ltmp.x *= 3L;
  1169.    ltmp.y *= 3L;
  1170.    LTRIGARG(ltmp.x);
  1171.    LTRIGARG(ltmp.y);
  1172.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1173.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1174.    ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
  1175.    ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
  1176.    LTRIGARG(ltmp.x);
  1177.    LTRIGARG(ltmp.y);
  1178.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1179.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1180.    lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
  1181.    lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
  1182.    if(plot == noplot)
  1183.    {
  1184.       iplot_orbit(lnew.x,lnew.y,1+row%colors);
  1185.       lold = lnew;
  1186.    }
  1187.    else
  1188.       LONGBAILOUT();
  1189.    /* PB above still the old way, is weird, see notes in FP popcorn case */
  1190.    return(0);
  1191. #endif
  1192. }
  1193.  
  1194. int MarksCplxMand(void)
  1195. {
  1196.    tmp.x = tempsqrx - tempsqry;
  1197.    tmp.y = 2*old.x*old.y;
  1198.    FPUcplxmul(&tmp, &coefficient, &new);
  1199.    new.x += floatparm->x;
  1200.    new.y += floatparm->y;
  1201.    return(floatbailout());
  1202. }
  1203.  
  1204. int SpiderfpFractal(void)
  1205. {
  1206.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1207.    new.x = tempsqrx - tempsqry + tmp.x;
  1208.    new.y = 2 * old.x * old.y + tmp.y;
  1209.    tmp.x = tmp.x/2 + new.x;
  1210.    tmp.y = tmp.y/2 + new.y;
  1211.    return(floatbailout());
  1212. }
  1213.  
  1214. SpiderFractal(void)
  1215. {
  1216. #ifndef XFRACT
  1217.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1218.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x;
  1219.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y;
  1220.    ltmp.x = (ltmp.x >> 1) + lnew.x;
  1221.    ltmp.y = (ltmp.y >> 1) + lnew.y;
  1222.    return(longbailout());
  1223. #endif
  1224. }
  1225.  
  1226. TetratefpFractal()
  1227. {
  1228.    /* Tetrate(XAXIS) { c=z=pixel: z=c^z, |z|<=(P1+3) } */
  1229.    new = ComplexPower(*floatparm,old);
  1230.    return(floatbailout());
  1231. }
  1232.  
  1233. ZXTrigPlusZFractal()
  1234. {
  1235. #ifndef XFRACT
  1236.    /* z = (p1*z*trig(z))+p2*z */
  1237.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)         */
  1238.    LCMPLXmult(lparm,ltmp,ltmp);      /* ltmp  = p1*trig(old)          */
  1239.    LCMPLXmult(lold,ltmp,ltmp2);      /* ltmp2 = p1*old*trig(old)      */
  1240.    LCMPLXmult(lparm2,lold,ltmp);     /* ltmp  = p2*old              */
  1241.    LCMPLXadd(ltmp2,ltmp,lnew);         /* lnew  = p1*trig(old) + p2*old */
  1242.    return(longbailout());
  1243. #endif
  1244. }
  1245.  
  1246. ScottZXTrigPlusZFractal()
  1247. {
  1248. #ifndef XFRACT
  1249.    /* z = (z*trig(z))+z */
  1250.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1251.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1252.    LCMPLXadd(lnew,lold,lnew);         /* lnew  = trig(old) + old */
  1253.    return(longbailout());
  1254. #endif
  1255. }
  1256.  
  1257. SkinnerZXTrigSubZFractal()
  1258. {
  1259. #ifndef XFRACT
  1260.    /* z = (z*trig(z))-z */
  1261.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1262.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1263.    LCMPLXsub(lnew,lold,lnew);         /* lnew  = trig(old) - old */
  1264.    return(longbailout());
  1265. #endif
  1266. }
  1267.  
  1268. ZXTrigPlusZfpFractal()
  1269. {
  1270.    /* z = (p1*z*trig(z))+p2*z */
  1271.    CMPLXtrig0(old,tmp);      /* tmp  = trig(old)         */
  1272.    CMPLXmult(parm,tmp,tmp);     /* tmp  = p1*trig(old)      */
  1273.    CMPLXmult(old,tmp,tmp2);     /* tmp2 = p1*old*trig(old)     */
  1274.    CMPLXmult(parm2,old,tmp);     /* tmp  = p2*old         */
  1275.    CMPLXadd(tmp2,tmp,new);     /* new  = p1*trig(old) + p2*old */
  1276.    return(floatbailout());
  1277. }
  1278.  
  1279. ScottZXTrigPlusZfpFractal()
  1280. {
  1281.    /* z = (z*trig(z))+z */
  1282.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1283.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1284.    CMPLXadd(new,old,new);     /* new  = trig(old) + old */
  1285.    return(floatbailout());
  1286. }
  1287.  
  1288. SkinnerZXTrigSubZfpFractal()
  1289. {
  1290.    /* z = (z*trig(z))-z */
  1291.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1292.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1293.    CMPLXsub(new,old,new);     /* new  = trig(old) - old */
  1294.    return(floatbailout());
  1295. }
  1296.  
  1297. Sqr1overTrigFractal()
  1298. {
  1299. #ifndef XFRACT
  1300.    /* z = sqr(1/trig(z)) */
  1301.    LCMPLXtrig0(lold,lold);
  1302.    LCMPLXrecip(lold,lold);
  1303.    LCMPLXsqr(lold,lnew);
  1304.    return(longbailout());
  1305. #endif
  1306. }
  1307.  
  1308. Sqr1overTrigfpFractal()
  1309. {
  1310.    /* z = sqr(1/trig(z)) */
  1311.    CMPLXtrig0(old,old);
  1312.    CMPLXrecip(old,old);
  1313.    CMPLXsqr(old,new);
  1314.    return(floatbailout());
  1315. }
  1316.  
  1317. TrigPlusTrigFractal()
  1318. {
  1319. #ifndef XFRACT
  1320.    /* z = trig(0,z)*p1+trig1(z)*p2 */
  1321.    LCMPLXtrig0(lold,ltmp);
  1322.    LCMPLXmult(lparm,ltmp,ltmp);
  1323.    LCMPLXtrig1(lold,ltmp2);
  1324.    LCMPLXmult(lparm2,ltmp2,lold);
  1325.    LCMPLXadd(ltmp,lold,lnew);
  1326.    return(longbailout());
  1327. #endif
  1328. }
  1329.  
  1330. TrigPlusTrigfpFractal()
  1331. {
  1332.    /* z = trig0(z)*p1+trig1(z)*p2 */
  1333.    CMPLXtrig0(old,tmp);
  1334.    CMPLXmult(parm,tmp,tmp);
  1335.    CMPLXtrig1(old,old);
  1336.    CMPLXmult(parm2,old,old);
  1337.    CMPLXadd(tmp,old,new);
  1338.    return(floatbailout());
  1339. }
  1340.  
  1341. /* The following four fractals are based on the idea of parallel
  1342.    or alternate calculations.  The shift is made when the mod
  1343.    reaches a given value.  JCO  5/6/92 */
  1344.  
  1345. LambdaTrigOrTrigFractal()
  1346. {
  1347. #ifndef XFRACT
  1348.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1349.           trig1(z)*p1 if mod(old) >= p2.x */
  1350.    if ((LCMPLXmod(lold)) < lparm2.x){
  1351.      LCMPLXtrig0(lold,ltmp);
  1352.      LCMPLXmult(*longparm,ltmp,lnew);}
  1353.    else{
  1354.      LCMPLXtrig1(lold,ltmp);
  1355.      LCMPLXmult(*longparm,ltmp,lnew);}
  1356.    return(longbailout());
  1357. #endif
  1358. }
  1359.  
  1360. LambdaTrigOrTrigfpFractal()
  1361. {
  1362.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1363.           trig1(z)*p1 if mod(old) >= p2.x */
  1364.    if (CMPLXmod(old) < parm2.x){
  1365.      CMPLXtrig0(old,old);
  1366.      FPUcplxmul(floatparm,&old,&new);}
  1367.    else{
  1368.      CMPLXtrig1(old,old);
  1369.      FPUcplxmul(floatparm,&old,&new);}
  1370.    return(floatbailout());
  1371. }
  1372.  
  1373. JuliaTrigOrTrigFractal()
  1374. {
  1375. #ifndef XFRACT
  1376.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1377.           trig1(z)+p1 if mod(old) >= p2.x */
  1378.    if (LCMPLXmod(lold) < lparm2.x){
  1379.      LCMPLXtrig0(lold,ltmp);
  1380.      LCMPLXadd(*longparm,ltmp,lnew);}
  1381.    else{
  1382.      LCMPLXtrig1(lold,ltmp);
  1383.      LCMPLXadd(*longparm,ltmp,lnew);}
  1384.    return(longbailout());
  1385. #endif
  1386. }
  1387.  
  1388. JuliaTrigOrTrigfpFractal()
  1389. {
  1390.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1391.           trig1(z)+p1 if mod(old) >= p2.x */
  1392.    if (CMPLXmod(old) < parm2.x){
  1393.      CMPLXtrig0(old,old);
  1394.      CMPLXadd(*floatparm,old,new);}
  1395.    else{
  1396.      CMPLXtrig1(old,old);
  1397.      CMPLXadd(*floatparm,old,new);}
  1398.    return(floatbailout());
  1399. }
  1400.  
  1401. int AplusOne, Ap1deg;
  1402. struct MP mpAplusOne, mpAp1deg;
  1403. struct MPC mpctmpparm; 
  1404.  
  1405. int MPCHalleyFractal()
  1406. {
  1407. #ifndef XFRACT
  1408.    /*  X(X^a - 1) = 0, Halley Map */
  1409.    /*  a = parm.x,  relaxation coeff. = parm.y,  epsilon = parm2.x  */
  1410.  
  1411. int ihal;
  1412. struct MPC mpcXtoAlessOne, mpcXtoA;
  1413. struct MPC mpcXtoAplusOne; /* a-1, a, a+1 */
  1414. struct MPC mpcFX, mpcF1prime, mpcF2prime, mpcHalnumer1;
  1415. struct MPC mpcHalnumer2, mpcHaldenom, mpctmp;
  1416.  
  1417.    MPOverflow = 0;
  1418.    mpcXtoAlessOne.x = mpcold.x;
  1419.    mpcXtoAlessOne.y = mpcold.y;
  1420.    for(ihal=2; ihal<degree; ihal++) {
  1421.      mpctmp.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1422.      mpctmp.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1423.      mpcXtoAlessOne.x = mpctmp.x;
  1424.      mpcXtoAlessOne.y = mpctmp.y;
  1425.    }
  1426.    mpcXtoA.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1427.    mpcXtoA.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1428.    mpcXtoAplusOne.x = *pMPsub(*pMPmul(mpcXtoA.x,mpcold.x),*pMPmul(mpcXtoA.y,mpcold.y));
  1429.    mpcXtoAplusOne.y = *pMPadd(*pMPmul(mpcXtoA.x,mpcold.y),*pMPmul(mpcXtoA.y,mpcold.x));
  1430.  
  1431.    mpcFX.x = *pMPsub(mpcXtoAplusOne.x, mpcold.x);
  1432.    mpcFX.y = *pMPsub(mpcXtoAplusOne.y, mpcold.y); /* FX = X^(a+1) - X  = F */
  1433.  
  1434.    mpcF2prime.x = *pMPmul(mpAp1deg, mpcXtoAlessOne.x); /* mpAp1deg in setup */
  1435.    mpcF2prime.y = *pMPmul(mpAp1deg, mpcXtoAlessOne.y);        /* F" */
  1436.  
  1437.    mpcF1prime.x = *pMPsub(*pMPmul(mpAplusOne, mpcXtoA.x), mpone);
  1438.    mpcF1prime.y = *pMPmul(mpAplusOne, mpcXtoA.y);                   /*  F'  */
  1439.  
  1440.    mpctmp.x = *pMPsub(*pMPmul(mpcF2prime.x,mpcFX.x),*pMPmul(mpcF2prime.y,mpcFX.y));
  1441.    mpctmp.y = *pMPadd(*pMPmul(mpcF2prime.x,mpcFX.y),*pMPmul(mpcF2prime.y,mpcFX.x));
  1442.    /*  F * F"  */
  1443.  
  1444.    mpcHaldenom.x = *pMPadd(mpcF1prime.x, mpcF1prime.x);
  1445.    mpcHaldenom.y = *pMPadd(mpcF1prime.y, mpcF1prime.y);      /*  2 * F'  */
  1446.  
  1447.    mpcHalnumer1 = MPCdiv(mpctmp, mpcHaldenom);        /*  F"F/2F'  */
  1448.    mpctmp.x = *pMPsub(mpcF1prime.x, mpcHalnumer1.x);
  1449.    mpctmp.y = *pMPsub(mpcF1prime.y, mpcHalnumer1.y); /*  F' - F"F/2F'  */
  1450.    mpcHalnumer2 = MPCdiv(mpcFX, mpctmp);
  1451.  
  1452.    mpctmp   =  MPCmul(mpctmpparm, mpcHalnumer2);  /* mpctmpparm is */   
  1453.                                                   /* relaxation coef. */
  1454. #if 0
  1455.    mpctmp.x = *pMPmul(mptmpparmy,mpcHalnumer2.x); /* mptmpparmy is */
  1456.    mpctmp.y = *pMPmul(mptmpparmy,mpcHalnumer2.y); /* relaxation coef. */
  1457.  
  1458.    mpcnew.x = *pMPsub(mpcold.x, mpctmp.x);
  1459.    mpcnew.y = *pMPsub(mpcold.y, mpctmp.y);
  1460.  
  1461.    new.x = *pMP2d(mpcnew.x);
  1462.    new.y = *pMP2d(mpcnew.y);
  1463. #endif
  1464.    mpcnew = MPCsub(mpcold, mpctmp);
  1465.    new    = MPC2cmplx(mpcnew);
  1466.    return(MPCHalleybailout()||MPOverflow);
  1467. #endif
  1468. }
  1469.  
  1470. HalleyFractal()
  1471. {
  1472.    /*  X(X^a - 1) = 0, Halley Map */
  1473.    /*  a = parm.x = degree, relaxation coeff. = parm.y, epsilon = parm2.x  */
  1474.  
  1475. int ihal;
  1476. _CMPLX XtoAlessOne, XtoA, XtoAplusOne; /* a-1, a, a+1 */
  1477. _CMPLX FX, F1prime, F2prime, Halnumer1, Halnumer2, Haldenom;
  1478. _CMPLX relax;
  1479.  
  1480.    XtoAlessOne = old;
  1481.    for(ihal=2; ihal<degree; ihal++) {
  1482.      FPUcplxmul(&old, &XtoAlessOne, &XtoAlessOne);
  1483.    }
  1484.    FPUcplxmul(&old, &XtoAlessOne, &XtoA);
  1485.    FPUcplxmul(&old, &XtoA, &XtoAplusOne);
  1486.  
  1487.    CMPLXsub(XtoAplusOne, old, FX);        /* FX = X^(a+1) - X  = F */
  1488.    F2prime.x = Ap1deg * XtoAlessOne.x; /* Ap1deg in setup */
  1489.    F2prime.y = Ap1deg * XtoAlessOne.y;        /* F" */
  1490.  
  1491.    F1prime.x = AplusOne * XtoA.x - 1.0;
  1492.    F1prime.y = AplusOne * XtoA.y;                             /*  F'  */
  1493.  
  1494.    FPUcplxmul(&F2prime, &FX, &Halnumer1);                  /*  F * F"  */
  1495.    Haldenom.x = F1prime.x + F1prime.x;
  1496.    Haldenom.y = F1prime.y + F1prime.y;                     /*  2 * F'  */
  1497.  
  1498.    FPUcplxdiv(&Halnumer1, &Haldenom, &Halnumer1);         /*  F"F/2F'  */
  1499.    CMPLXsub(F1prime, Halnumer1, Halnumer2);          /*  F' - F"F/2F'  */
  1500.    FPUcplxdiv(&FX, &Halnumer2, &Halnumer2);
  1501.    /* parm.y is relaxation coef. */
  1502.    /* new.x = old.x - (parm.y * Halnumer2.x); 
  1503.    new.y = old.y - (parm.y * Halnumer2.y); */
  1504.    relax.x = parm.y;
  1505.    relax.y = param[3];
  1506.    FPUcplxmul(&relax, &Halnumer2, &Halnumer2);
  1507.    new.x = old.x - Halnumer2.x;
  1508.    new.y = old.y - Halnumer2.y;
  1509.    return(Halleybailout());
  1510. }
  1511.  
  1512. LongPhoenixFractal()
  1513. {
  1514. #ifndef XFRACT
  1515. /* z(n+1) = z(n)^2 + p + qy(n),  y(n+1) = z(n) */
  1516.    ltmp.x = multiply(lold.x, lold.y, bitshift);
  1517.    lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(longparm->y,ltmp2.x,bitshift);
  1518.    lnew.y = (ltmp.x + ltmp.x) + multiply(longparm->y,ltmp2.y,bitshift);
  1519.    ltmp2 = lold; /* set ltmp2 to Y value */
  1520.    return(longbailout());
  1521. #endif
  1522. }
  1523.  
  1524. PhoenixFractal()
  1525. {
  1526. /* z(n+1) = z(n)^2 + p + qy(n),  y(n+1) = z(n) */
  1527.    tmp.x = old.x * old.y;
  1528.    new.x = tempsqrx - tempsqry + floatparm->x + (floatparm->y * tmp2.x);
  1529.    new.y = (tmp.x + tmp.x) + (floatparm->y * tmp2.y);
  1530.    tmp2 = old; /* set tmp2 to Y value */
  1531.    return(floatbailout());
  1532. }
  1533.  
  1534. LongPhoenixFractalcplx()
  1535. {
  1536. #ifndef XFRACT
  1537. /* z(n+1) = z(n)^2 + p + qy(n),  y(n+1) = z(n) */
  1538.    ltmp.x = multiply(lold.x, lold.y, bitshift);
  1539.    lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(lparm2.x,ltmp2.x,bitshift)-multiply(lparm2.y,ltmp2.y,bitshift);
  1540.    lnew.y = (ltmp.x + ltmp.x)+longparm->y+multiply(lparm2.x,ltmp2.y,bitshift)+multiply(lparm2.y,ltmp2.x,bitshift);
  1541.    ltmp2 = lold; /* set ltmp2 to Y value */
  1542.    return(longbailout());
  1543. #endif
  1544. }
  1545.  
  1546. PhoenixFractalcplx()
  1547. {
  1548. /* z(n+1) = z(n)^2 + p1 + p2*y(n),  y(n+1) = z(n) */
  1549.    tmp.x = old.x * old.y;
  1550.    new.x = tempsqrx - tempsqry + floatparm->x + (parm2.x * tmp2.x) - (parm2.y * tmp2.y);
  1551.    new.y = (tmp.x + tmp.x) + floatparm->y + (parm2.x * tmp2.y) + (parm2.y * tmp2.x);
  1552.    tmp2 = old; /* set tmp2 to Y value */
  1553.    return(floatbailout());
  1554. }
  1555.  
  1556. LongPhoenixPlusFractal()
  1557. {
  1558. #ifndef XFRACT
  1559. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1560. int i;
  1561. _LCMPLX loldplus, lnewminus;
  1562.    loldplus = lold;
  1563.    ltmp = lold;
  1564.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1565.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
  1566.    }
  1567.    loldplus.x += longparm->x;
  1568.    LCMPLXmult(ltmp, loldplus, lnewminus);
  1569.    lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
  1570.    lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
  1571.    ltmp2 = lold; /* set ltmp2 to Y value */
  1572.    return(longbailout());
  1573. #endif
  1574. }
  1575.  
  1576. PhoenixPlusFractal()
  1577. {
  1578. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1579. int i;
  1580. _CMPLX oldplus, newminus;
  1581.    oldplus = old;
  1582.    tmp = old;
  1583.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1584.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
  1585.    }
  1586.    oldplus.x += floatparm->x;
  1587.    FPUcplxmul(&tmp, &oldplus, &newminus);
  1588.    new.x = newminus.x + (floatparm->y * tmp2.x);
  1589.    new.y = newminus.y + (floatparm->y * tmp2.y);
  1590.    tmp2 = old; /* set tmp2 to Y value */
  1591.    return(floatbailout());
  1592. }
  1593.  
  1594. LongPhoenixMinusFractal()
  1595. {
  1596. #ifndef XFRACT
  1597. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1598. int i;
  1599. _LCMPLX loldsqr, lnewminus;
  1600.    LCMPLXmult(lold,lold,loldsqr);
  1601.    ltmp = lold;
  1602.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1603.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
  1604.    }
  1605.    loldsqr.x += longparm->x;
  1606.    LCMPLXmult(ltmp, loldsqr, lnewminus);
  1607.    lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
  1608.    lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
  1609.    ltmp2 = lold; /* set ltmp2 to Y value */
  1610.    return(longbailout());
  1611. #endif
  1612. }
  1613.  
  1614. PhoenixMinusFractal()
  1615. {
  1616. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1617. int i;
  1618. _CMPLX oldsqr, newminus;
  1619.    FPUcplxmul(&old, &old, &oldsqr);
  1620.    tmp = old;
  1621.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1622.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
  1623.    }
  1624.    oldsqr.x += floatparm->x;
  1625.    FPUcplxmul(&tmp, &oldsqr, &newminus);
  1626.    new.x = newminus.x + (floatparm->y * tmp2.x);
  1627.    new.y = newminus.y + (floatparm->y * tmp2.y);
  1628.    tmp2 = old; /* set tmp2 to Y value */
  1629.    return(floatbailout());
  1630. }
  1631.  
  1632. LongPhoenixCplxPlusFractal()
  1633. {
  1634. #ifndef XFRACT
  1635. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1636. int i;
  1637. _LCMPLX loldplus, lnewminus;
  1638.    loldplus = lold;
  1639.    ltmp = lold;
  1640.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1641.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
  1642.    }
  1643.    loldplus.x += longparm->x;
  1644.    loldplus.y += longparm->y;
  1645.    LCMPLXmult(ltmp, loldplus, lnewminus);
  1646.    LCMPLXmult(lparm2, ltmp2, ltmp);
  1647.    lnew.x = lnewminus.x + ltmp.x;
  1648.    lnew.y = lnewminus.y + ltmp.y;
  1649.    ltmp2 = lold; /* set ltmp2 to Y value */
  1650.    return(longbailout());
  1651. #endif
  1652. }
  1653.  
  1654. PhoenixCplxPlusFractal()
  1655. {
  1656. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1657. int i;
  1658. _CMPLX oldplus, newminus;
  1659.    oldplus = old;
  1660.    tmp = old;
  1661.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1662.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
  1663.    }
  1664.    oldplus.x += floatparm->x;
  1665.    oldplus.y += floatparm->y;
  1666.    FPUcplxmul(&tmp, &oldplus, &newminus);
  1667.    FPUcplxmul(&parm2, &tmp2, &tmp);
  1668.    new.x = newminus.x + tmp.x;
  1669.    new.y = newminus.y + tmp.y;
  1670.    tmp2 = old; /* set tmp2 to Y value */
  1671.    return(floatbailout());
  1672. }
  1673.  
  1674. LongPhoenixCplxMinusFractal()
  1675. {
  1676. #ifndef XFRACT
  1677. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1678. int i;
  1679. _LCMPLX loldsqr, lnewminus;
  1680.    LCMPLXmult(lold,lold,loldsqr);
  1681.    ltmp = lold;
  1682.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1683.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
  1684.    }
  1685.    loldsqr.x += longparm->x;
  1686.    loldsqr.y += longparm->y;
  1687.    LCMPLXmult(ltmp, loldsqr, lnewminus);
  1688.    LCMPLXmult(lparm2, ltmp2, ltmp);
  1689.    lnew.x = lnewminus.x + ltmp.x;
  1690.    lnew.y = lnewminus.y + ltmp.y;
  1691.    ltmp2 = lold; /* set ltmp2 to Y value */
  1692.    return(longbailout());
  1693. #endif
  1694. }
  1695.  
  1696. PhoenixCplxMinusFractal()
  1697. {
  1698. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1699. int i;
  1700. _CMPLX oldsqr, newminus;
  1701.    FPUcplxmul(&old, &old, &oldsqr);
  1702.    tmp = old;
  1703.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1704.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
  1705.    }
  1706.    oldsqr.x += floatparm->x;
  1707.    oldsqr.y += floatparm->y;
  1708.    FPUcplxmul(&tmp, &oldsqr, &newminus);
  1709.    FPUcplxmul(&parm2, &tmp2, &tmp);
  1710.    new.x = newminus.x + tmp.x;
  1711.    new.y = newminus.y + tmp.y;
  1712.    tmp2 = old; /* set tmp2 to Y value */
  1713.    return(floatbailout());
  1714. }
  1715.  
  1716. ScottTrigPlusTrigFractal()
  1717. {
  1718. #ifndef XFRACT
  1719.    /* z = trig0(z)+trig1(z) */
  1720.    LCMPLXtrig0(lold,ltmp);
  1721.    LCMPLXtrig1(lold,lold);
  1722.    LCMPLXadd(ltmp,lold,lnew);
  1723.    return(longbailout());
  1724. #endif
  1725. }
  1726.  
  1727. ScottTrigPlusTrigfpFractal()
  1728. {
  1729.    /* z = trig0(z)+trig1(z) */
  1730.    CMPLXtrig0(old,tmp);
  1731.    CMPLXtrig1(old,tmp2);
  1732.    CMPLXadd(tmp,tmp2,new);
  1733.    return(floatbailout());
  1734. }
  1735.  
  1736. SkinnerTrigSubTrigFractal()
  1737. {
  1738. #ifndef XFRACT
  1739.    /* z = trig(0,z)-trig1(z) */
  1740.    LCMPLXtrig0(lold,ltmp);
  1741.    LCMPLXtrig1(lold,ltmp2);
  1742.    LCMPLXsub(ltmp,ltmp2,lnew);
  1743.    return(longbailout());
  1744. #endif
  1745. }
  1746.  
  1747. SkinnerTrigSubTrigfpFractal()
  1748. {
  1749.    /* z = trig0(z)-trig1(z) */
  1750.    CMPLXtrig0(old,tmp);
  1751.    CMPLXtrig1(old,tmp2);
  1752.    CMPLXsub(tmp,tmp2,new);
  1753.    return(floatbailout());
  1754. }
  1755.  
  1756.  
  1757. TrigXTrigfpFractal(void)
  1758. {
  1759.    /* z = trig0(z)*trig1(z) */
  1760.    CMPLXtrig0(old,tmp);
  1761.    CMPLXtrig1(old,old);
  1762.    CMPLXmult(tmp,old,new);
  1763.    return(floatbailout());
  1764. }
  1765.  
  1766.  
  1767.  /* call float version of fractal if integer math overflow */
  1768. static TryFloatFractal(int (*fpFractal)(void))
  1769. {
  1770.    overflow=0;
  1771.    /* lold had better not be changed! */
  1772.    old.x = lold.x; old.x /= fudge;
  1773.    old.y = lold.y; old.y /= fudge;
  1774.    tempsqrx = sqr(old.x);
  1775.    tempsqry = sqr(old.y);
  1776.    fpFractal();
  1777.    lnew.x = (long)(new.x*fudge);
  1778.    lnew.y = (long)(new.y*fudge);
  1779.    return(0);
  1780. }
  1781.  
  1782. TrigXTrigFractal()
  1783. {
  1784. #ifndef XFRACT
  1785.    _LCMPLX ltmp2;
  1786.    /* z = trig0(z)*trig1(z) */
  1787.    LCMPLXtrig0(lold,ltmp);
  1788.    LCMPLXtrig1(lold,ltmp2);
  1789.    LCMPLXmult(ltmp,ltmp2,lnew);
  1790.    if(overflow)
  1791.       TryFloatFractal(TrigXTrigfpFractal);
  1792.    return(longbailout());
  1793. #endif
  1794. }
  1795.  
  1796. /********************************************************************/
  1797. /*  Next six orbit functions are one type - extra functions are     */
  1798. /*    special cases written for speed.                    */
  1799. /********************************************************************/
  1800.  
  1801. TrigPlusSqrFractal() /* generalization of Scott and Skinner types */
  1802. {
  1803. #ifndef XFRACT
  1804.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1805.    LCMPLXtrig0(lold,ltmp);     /* ltmp = trig(lold)               */
  1806.    LCMPLXmult(lparm,ltmp,lnew); /* lnew = lparm*trig(lold)            */
  1807.    LCMPLXsqr_old(ltmp);     /* ltmp = sqr(lold)                */
  1808.    LCMPLXmult(lparm2,ltmp,ltmp);/* ltmp = lparm2*sqr(lold)            */
  1809.    LCMPLXadd(lnew,ltmp,lnew);    /* lnew = lparm*trig(lold)+lparm2*sqr(lold) */
  1810.    return(longbailout());
  1811. #endif
  1812. }
  1813.  
  1814. TrigPlusSqrfpFractal() /* generalization of Scott and Skinner types */
  1815. {
  1816.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1817.    CMPLXtrig0(old,tmp);     /* tmp = trig(old)               */
  1818.    CMPLXmult(parm,tmp,new); /* new = parm*trig(old)           */
  1819.    CMPLXsqr_old(tmp);         /* tmp = sqr(old)                */
  1820.    CMPLXmult(parm2,tmp,tmp2); /* tmp = parm2*sqr(old)             */
  1821.    CMPLXadd(new,tmp2,new);    /* new = parm*trig(old)+parm2*sqr(old) */
  1822.    return(floatbailout());
  1823. }
  1824.  
  1825. ScottTrigPlusSqrFractal()
  1826. {
  1827. #ifndef XFRACT
  1828.    /*  { z=pixel: z=trig(z)+sqr(z), |z|<BAILOUT } */
  1829.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1830.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1831.    LCMPLXadd(ltmp,lnew,lnew);  /* lnew = trig(lold)+sqr(lold) */
  1832.    return(longbailout());
  1833. #endif
  1834. }
  1835.  
  1836. ScottTrigPlusSqrfpFractal() /* float version */
  1837. {
  1838.    /* { z=pixel: z=sin(z)+sqr(z), |z|<BAILOUT } */
  1839.    CMPLXtrig0(old,new);       /* new = trig(old)      */
  1840.    CMPLXsqr_old(tmp);           /* tmp = sqr(old)       */
  1841.    CMPLXadd(new,tmp,new);      /* new = trig(old)+sqr(old) */
  1842.    return(floatbailout());
  1843. }
  1844.  
  1845. SkinnerTrigSubSqrFractal()
  1846. {
  1847. #ifndef XFRACT
  1848.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT }           */
  1849.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1850.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1851.    LCMPLXsub(lnew,ltmp,lnew);  /* lnew = trig(lold)-sqr(lold) */
  1852.    return(longbailout());
  1853. #endif
  1854. }
  1855.  
  1856. SkinnerTrigSubSqrfpFractal()
  1857. {
  1858.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
  1859.    CMPLXtrig0(old,new);       /* new = trig(old) */
  1860.    CMPLXsqr_old(tmp);           /* old = sqr(old)  */
  1861.    CMPLXsub(new,tmp,new);      /* new = trig(old)-sqr(old) */
  1862.    return(floatbailout());
  1863. }
  1864.  
  1865. TrigZsqrdfpFractal(void)
  1866. {
  1867.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1868.    CMPLXsqr_old(tmp);
  1869.    CMPLXtrig0(tmp,new);
  1870.    return(floatbailout());
  1871. }
  1872.  
  1873. TrigZsqrdFractal() /* this doesn't work very well */
  1874. {
  1875. #ifndef XFRACT
  1876.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1877.    LCMPLXsqr_old(ltmp);
  1878.    if(labs(ltmp.x) > l16triglim || labs(ltmp.y) > l16triglim)
  1879.       overflow = 1;
  1880.    LCMPLXtrig0(ltmp,lnew);
  1881.    if(overflow)
  1882.       TryFloatFractal(TrigZsqrdfpFractal);
  1883.    return(longbailout());
  1884. #endif
  1885. }
  1886.  
  1887. SqrTrigFractal()
  1888. {
  1889. #ifndef XFRACT
  1890.    /* { z=pixel: z=sqr(trig(z)), |z|<TEST} */
  1891.    LCMPLXtrig0(lold,ltmp);
  1892.    LCMPLXsqr(ltmp,lnew);
  1893.    return(longbailout());
  1894. #endif
  1895. }
  1896.  
  1897. SqrTrigfpFractal()
  1898. {
  1899.    /* SZSB(XYAXIS) { z=pixel, TEST=(p1+3): z=sin(z)*sin(z), |z|<TEST} */
  1900.    CMPLXtrig0(old,tmp);
  1901.    CMPLXsqr(tmp,new);
  1902.    return(floatbailout());
  1903. }
  1904.  
  1905.  
  1906. Magnet1Fractal()    /*      Z = ((Z**2 + C - 1)/(2Z + C - 2))**2      */
  1907.   {              /*  In "Beauty of Fractals", code by Kev Allen. */
  1908.     _CMPLX top, bot, tmp;
  1909.     double div;
  1910.  
  1911.     top.x = tempsqrx - tempsqry + floatparm->x - 1; /* top = Z**2+C-1 */
  1912.     top.y = old.x * old.y;
  1913.     top.y = top.y + top.y + floatparm->y;
  1914.  
  1915.     bot.x = old.x + old.x + floatparm->x - 2;        /* bot = 2*Z+C-2  */
  1916.     bot.y = old.y + old.y + floatparm->y;
  1917.  
  1918.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1919.     if (div < FLT_MIN) return(1);
  1920.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1921.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1922.  
  1923.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1924.     new.y = tmp.x * tmp.y;
  1925.     new.y += new.y;
  1926.  
  1927.     return(floatbailout());
  1928.   }
  1929.  
  1930. Magnet2Fractal()  /* Z = ((Z**3 + 3(C-1)Z + (C-1)(C-2)    ) /     */
  1931.             /*         (3Z**2 + 3(C-2)Z + (C-1)(C-2)+1) )**2  */
  1932.   {            /*     In "Beauty of Fractals", code by Kev Allen.  */
  1933.     _CMPLX top, bot, tmp;
  1934.     double div;
  1935.  
  1936.     top.x = old.x * (tempsqrx-tempsqry-tempsqry-tempsqry + T_Cm1.x)
  1937.       - old.y * T_Cm1.y + T_Cm1Cm2.x;
  1938.     top.y = old.y * (tempsqrx+tempsqrx+tempsqrx-tempsqry + T_Cm1.x)
  1939.       + old.x * T_Cm1.y + T_Cm1Cm2.y;
  1940.  
  1941.     bot.x = tempsqrx - tempsqry;
  1942.     bot.x = bot.x + bot.x + bot.x
  1943.       + old.x * T_Cm2.x - old.y * T_Cm2.y
  1944.       + T_Cm1Cm2.x + 1.0;
  1945.     bot.y = old.x * old.y;
  1946.     bot.y += bot.y;
  1947.     bot.y = bot.y + bot.y + bot.y
  1948.       + old.x * T_Cm2.y + old.y * T_Cm2.x
  1949.       + T_Cm1Cm2.y;
  1950.  
  1951.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1952.     if (div < FLT_MIN) return(1);
  1953.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1954.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1955.  
  1956.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1957.     new.y = tmp.x * tmp.y;
  1958.     new.y += new.y;
  1959.  
  1960.     return(floatbailout());
  1961.   }
  1962.  
  1963. LambdaTrigFractal()
  1964. {
  1965. #ifndef XFRACT
  1966.    LONGXYTRIGBAILOUT();
  1967.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1968.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1969.    lold = lnew;
  1970.    return(0);
  1971. #endif
  1972. }
  1973.  
  1974. LambdaTrigfpFractal()
  1975. {
  1976.    FLOATXYTRIGBAILOUT();
  1977.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1978.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1979.    old = new;
  1980.    return(0);
  1981. }
  1982.  
  1983. /* bailouts are different for different trig functions */
  1984. LambdaTrigFractal1()
  1985. {
  1986. #ifndef XFRACT
  1987.    LONGTRIGBAILOUT(); /* sin,cos */
  1988.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1989.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1990.    lold = lnew;
  1991.    return(0);
  1992. #endif
  1993. }
  1994.  
  1995. LambdaTrigfpFractal1()
  1996. {
  1997.    FLOATTRIGBAILOUT(); /* sin,cos */
  1998.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1999.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  2000.    old = new;
  2001.    return(0);
  2002. }
  2003.  
  2004. LambdaTrigFractal2()
  2005. {
  2006. #ifndef XFRACT
  2007.    LONGHTRIGBAILOUT(); /* sinh,cosh */
  2008.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  2009.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  2010.    lold = lnew;
  2011.    return(0);
  2012. #endif
  2013. }
  2014.  
  2015. LambdaTrigfpFractal2()
  2016. {
  2017. #ifndef XFRACT
  2018.    FLOATHTRIGBAILOUT(); /* sinh,cosh */
  2019.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  2020.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  2021.    old = new;
  2022.    return(0);
  2023. #endif
  2024. }
  2025.  
  2026. ManOWarFractal()
  2027. {
  2028. #ifndef XFRACT
  2029.    /* From Art Matrix via Lee Skinner */
  2030.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x + longparm->x;
  2031.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y + longparm->y;
  2032.    ltmp = lold;
  2033.    return(longbailout());
  2034. #endif
  2035. }
  2036.  
  2037. ManOWarfpFractal()
  2038. {
  2039.    /* From Art Matrix via Lee Skinner */
  2040.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  2041.    new.x = tempsqrx - tempsqry + tmp.x + floatparm->x;
  2042.    new.y = 2.0 * old.x * old.y + tmp.y + floatparm->y;
  2043.    tmp = old;
  2044.    return(floatbailout());
  2045. }
  2046.  
  2047. /*
  2048.    MarksMandelPwr (XAXIS) {
  2049.       z = pixel, c = z ^ (z - 1):
  2050.      z = c * sqr(z) + pixel,
  2051.       |z| <= 4
  2052.    }
  2053. */
  2054.  
  2055. MarksMandelPwrfpFractal()
  2056. {
  2057.    CMPLXtrig0(old,new);
  2058.    CMPLXmult(tmp,new,new);
  2059.    new.x += floatparm->x;
  2060.    new.y += floatparm->y;
  2061.    return(floatbailout());
  2062. }
  2063.  
  2064. MarksMandelPwrFractal()
  2065. {
  2066. #ifndef XFRACT
  2067.    LCMPLXtrig0(lold,lnew);
  2068.    LCMPLXmult(ltmp,lnew,lnew);
  2069.    lnew.x += longparm->x;
  2070.    lnew.y += longparm->y;
  2071.    return(longbailout());
  2072. #endif
  2073. }
  2074.  
  2075. /* I was coding Marksmandelpower and failed to use some temporary
  2076.    variables. The result was nice, and since my name is not on any fractal,
  2077.    I thought I would immortalize myself with this error!
  2078.         Tim Wegner */
  2079.  
  2080.  
  2081. TimsErrorfpFractal()
  2082. {
  2083.    CMPLXtrig0(old,new);
  2084.    new.x = new.x * tmp.x - new.y * tmp.y;
  2085.    new.y = new.x * tmp.y - new.y * tmp.x;
  2086.    new.x += floatparm->x;
  2087.    new.y += floatparm->y;
  2088.    return(floatbailout());
  2089. }
  2090. TimsErrorFractal()
  2091. {
  2092. #ifndef XFRACT
  2093.    LCMPLXtrig0(lold,lnew);
  2094.    lnew.x = multiply(lnew.x,ltmp.x,bitshift)-multiply(lnew.y,ltmp.y,bitshift);
  2095.    lnew.y = multiply(lnew.x,ltmp.y,bitshift)-multiply(lnew.y,ltmp.x,bitshift);
  2096.    lnew.x += longparm->x;
  2097.    lnew.y += longparm->y;
  2098.    return(longbailout());
  2099. #endif
  2100. }
  2101.  
  2102. CirclefpFractal()
  2103. {
  2104.    long i;
  2105.    i = (long)(param[0]*(tempsqrx+tempsqry));
  2106.    coloriter = i%colors;
  2107.    return(1);
  2108. }
  2109. /*
  2110. CirclelongFractal()
  2111. {
  2112.    long i;
  2113.    i = multiply(lparm.x,(ltempsqrx+ltempsqry),bitshift);
  2114.    i = i >> bitshift;
  2115.    coloriter = i%colors);
  2116.    return(1);
  2117. }
  2118. */
  2119.  
  2120. /* -------------------------------------------------------------------- */
  2121. /*        Initialization (once per pixel) routines                        */
  2122. /* -------------------------------------------------------------------- */
  2123.  
  2124. #ifdef XFRACT
  2125. /* this code translated to asm - lives in newton.asm */
  2126. /* transform points with reciprocal function */
  2127. void invertz2(_CMPLX *z)
  2128. {
  2129.    z->x = dx0[col]+dx1[row];
  2130.    z->y = dy0[row]+dy1[col];
  2131.    z->x -= f_xcenter; z->y -= f_ycenter;  /* Normalize values to center of circle */
  2132.  
  2133.    tempsqrx = sqr(z->x) + sqr(z->y);  /* Get old radius */
  2134.    if(fabs(tempsqrx) > FLT_MIN)
  2135.       tempsqrx = f_radius / tempsqrx;
  2136.    else
  2137.       tempsqrx = FLT_MAX;   /* a big number, but not TOO big */
  2138.    z->x *= tempsqrx;      z->y *= tempsqrx;     /* Perform inversion */
  2139.    z->x += f_xcenter; z->y += f_ycenter; /* Renormalize */
  2140. }
  2141. #endif
  2142.  
  2143. int long_julia_per_pixel()
  2144. {
  2145. #ifndef XFRACT
  2146.    /* integer julia types */
  2147.    /* lambda */
  2148.    /* barnsleyj1 */
  2149.    /* barnsleyj2 */
  2150.    /* sierpinski */
  2151.    if(invert)
  2152.    {
  2153.       /* invert */
  2154.       invertz2(&old);
  2155.  
  2156.       /* watch out for overflow */
  2157.       if(sqr(old.x)+sqr(old.y) >= 127)
  2158.       {
  2159.      old.x = 8;  /* value to bail out in one iteration */
  2160.      old.y = 8;
  2161.       }
  2162.  
  2163.       /* convert to fudged longs */
  2164.       lold.x = (long)(old.x*fudge);
  2165.       lold.y = (long)(old.y*fudge);
  2166.    }
  2167.    else
  2168.    {
  2169.       lold.x = lx0[col]+lx1[row];
  2170.       lold.y = ly0[row]+ly1[col];
  2171.    }
  2172.    return(0);
  2173. #else
  2174.    printf("Called long_julia_per_pixel\n");
  2175.    exit(0);
  2176. #endif
  2177. }
  2178.  
  2179. int long_richard8_per_pixel()
  2180. {
  2181. #ifndef XFRACT
  2182.     long_mandel_per_pixel();
  2183.     LCMPLXtrig1(*longparm,ltmp);
  2184.     LCMPLXmult(ltmp,lparm2,ltmp);
  2185.     return(1);
  2186. #endif
  2187. }
  2188.  
  2189. int long_mandel_per_pixel()
  2190. {
  2191. #ifndef XFRACT
  2192.    /* integer mandel types */
  2193.    /* barnsleym1 */
  2194.    /* barnsleym2 */
  2195.    linit.x = lx0[col]+lx1[row];
  2196.  
  2197.    if(invert)
  2198.    {
  2199.       /* invert */
  2200.       invertz2(&init);
  2201.  
  2202.       /* watch out for overflow */
  2203.       if(sqr(init.x)+sqr(init.y) >= 127)
  2204.       {
  2205.      init.x = 8;  /* value to bail out in one iteration */
  2206.      init.y = 8;
  2207.       }
  2208.  
  2209.       /* convert to fudged longs */
  2210.       linit.x = (long)(init.x*fudge);
  2211.       linit.y = (long)(init.y*fudge);
  2212.    }
  2213.  
  2214.    if(useinitorbit == 1)
  2215.       lold = linitorbit;
  2216.    else
  2217.       lold = linit;
  2218.  
  2219.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2220.    lold.y += lparm.y;
  2221.    return(1); /* 1st iteration has been done */
  2222. #else
  2223.    printf("Called long_mandel_per_pixel\n");
  2224.    exit(0);
  2225. #endif
  2226. }
  2227.  
  2228. int julia_per_pixel()
  2229. {
  2230.    /* julia */
  2231.  
  2232.    if(invert)
  2233.    {
  2234.       /* invert */
  2235.       invertz2(&old);
  2236.  
  2237.       /* watch out for overflow */
  2238.       if(bitshift <= 24)
  2239.      if (sqr(old.x)+sqr(old.y) >= 127)
  2240.      {
  2241.         old.x = 8;    /* value to bail out in one iteration */
  2242.         old.y = 8;
  2243.      }
  2244.       if(bitshift >  24)
  2245.      if (sqr(old.x)+sqr(old.y) >= 4.0)
  2246.      {
  2247.         old.x = 2;    /* value to bail out in one iteration */
  2248.         old.y = 2;
  2249.      }
  2250.  
  2251.       /* convert to fudged longs */
  2252.       lold.x = (long)(old.x*fudge);
  2253.       lold.y = (long)(old.y*fudge);
  2254.    }
  2255.    else
  2256.    {
  2257.       lold.x = lx0[col]+lx1[row];
  2258.       lold.y = ly0[row]+ly1[col];
  2259.    }
  2260.  
  2261.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2262.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2263.    ltmp = lold;
  2264.    return(0);
  2265. }
  2266.  
  2267. marks_mandelpwr_per_pixel()
  2268. {
  2269. #ifndef XFRACT
  2270.    mandel_per_pixel();
  2271.    ltmp = lold;
  2272.    ltmp.x -= fudge;
  2273.    LCMPLXpwr(lold,ltmp,ltmp);
  2274.    return(1);
  2275. #endif
  2276. }
  2277.  
  2278. int mandel_per_pixel()
  2279. {
  2280.    /* mandel */
  2281.  
  2282.    if(invert)
  2283.    {
  2284.       invertz2(&init);
  2285.  
  2286.       /* watch out for overflow */
  2287.       if(bitshift <= 24)
  2288.      if (sqr(init.x)+sqr(init.y) >= 127)
  2289.      {
  2290.         init.x = 8;  /* value to bail out in one iteration */
  2291.         init.y = 8;
  2292.      }
  2293.       if(bitshift >  24)
  2294.      if (sqr(init.x)+sqr(init.y) >= 4)
  2295.      {
  2296.         init.x = 2;  /* value to bail out in one iteration */
  2297.         init.y = 2;
  2298.      }
  2299.  
  2300.       /* convert to fudged longs */
  2301.       linit.x = (long)(init.x*fudge);
  2302.       linit.y = (long)(init.y*fudge);
  2303.    }
  2304.    else
  2305.       linit.x = lx0[col]+lx1[row];
  2306.    switch (fractype)
  2307.      {
  2308.     case MANDELLAMBDA:        /* Critical Value 0.5 + 0.0i  */
  2309.         lold.x = FgHalf;
  2310.         lold.y = 0;
  2311.         break;
  2312.     default:
  2313.         lold = linit;
  2314.         break;
  2315.       }
  2316.  
  2317.    /* alter init value */
  2318.    if(useinitorbit == 1)
  2319.       lold = linitorbit;
  2320.    else if(useinitorbit == 2)
  2321.       lold = linit;
  2322.  
  2323.    if(inside == -60 || inside == -61)
  2324.    {
  2325.       /* kludge to match "Beauty of Fractals" picture since we start
  2326.      Mandelbrot iteration with init rather than 0 */
  2327.       lold.x = lparm.x; /* initial pertubation of parameters set */
  2328.       lold.y = lparm.y;
  2329.       coloriter = -1;
  2330.    }
  2331.    else
  2332.    {
  2333.       lold.x += lparm.x; /* initial pertubation of parameters set */
  2334.       lold.y += lparm.y;
  2335.    }
  2336.    ltmp = linit; /* for spider */
  2337.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2338.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2339.    return(1); /* 1st iteration has been done */
  2340. }
  2341.  
  2342.  
  2343. int marksmandel_per_pixel()
  2344. {
  2345.    /* marksmandel */
  2346.  
  2347.    if(invert)
  2348.    {
  2349.       invertz2(&init);
  2350.  
  2351.       /* watch out for overflow */
  2352.       if(sqr(init.x)+sqr(init.y) >= 127)
  2353.       {
  2354.      init.x = 8;  /* value to bail out in one iteration */
  2355.      init.y = 8;
  2356.       }
  2357.  
  2358.       /* convert to fudged longs */
  2359.       linit.x = (long)(init.x*fudge);
  2360.       linit.y = (long)(init.y*fudge);
  2361.    }
  2362.    else
  2363.       linit.x = lx0[col]+lx1[row];
  2364.  
  2365.    if(useinitorbit == 1)
  2366.       lold = linitorbit;
  2367.    else
  2368.       lold = linit;
  2369.  
  2370.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2371.    lold.y += lparm.y;
  2372.  
  2373.    if(c_exp > 3)
  2374.       lcpower(&lold,c_exp-1,&lcoefficient,bitshift);
  2375.    else if(c_exp == 3) {
  2376.       lcoefficient.x = multiply(lold.x, lold.x, bitshift)
  2377.      - multiply(lold.y, lold.y, bitshift);
  2378.       lcoefficient.y = multiply(lold.x, lold.y, bitshiftless1);
  2379.    }
  2380.    else if(c_exp == 2)
  2381.       lcoefficient = lold;
  2382.    else if(c_exp < 2) {
  2383.       lcoefficient.x = 1L << bitshift;
  2384.       lcoefficient.y = 0L;
  2385.    }
  2386.  
  2387.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2388.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2389.    return(1); /* 1st iteration has been done */
  2390. }
  2391.  
  2392. int marksmandelfp_per_pixel()
  2393. {
  2394.    /* marksmandel */
  2395.  
  2396.    if(invert)
  2397.       invertz2(&init);
  2398.    else
  2399.       init.x = dx0[col]+dx1[row];
  2400.  
  2401.    if(useinitorbit == 1)
  2402.       old = initorbit;
  2403.    else
  2404.       old = init;
  2405.  
  2406.    old.x += parm.x;     /* initial pertubation of parameters set */
  2407.    old.y += parm.y;
  2408.  
  2409.    tempsqrx = sqr(old.x);
  2410.    tempsqry = sqr(old.y);
  2411.  
  2412.    if(c_exp > 3)
  2413.       cpower(&old,c_exp-1,&coefficient);
  2414.    else if(c_exp == 3) {
  2415.       coefficient.x = tempsqrx - tempsqry;
  2416.       coefficient.y = old.x * old.y * 2;
  2417.    }
  2418.    else if(c_exp == 2)
  2419.       coefficient = old;
  2420.    else if(c_exp < 2) {
  2421.       coefficient.x = 1.0;
  2422.       coefficient.y = 0.0;
  2423.    }
  2424.  
  2425.    return(1); /* 1st iteration has been done */
  2426. }
  2427.  
  2428. marks_mandelpwrfp_per_pixel()
  2429. {
  2430.    mandelfp_per_pixel();
  2431.    tmp = old;
  2432.    tmp.x -= 1;
  2433.    CMPLXpwr(old,tmp,tmp);
  2434.    return(1);
  2435. }
  2436.  
  2437. int mandelfp_per_pixel()
  2438. {
  2439.    /* floating point mandelbrot */
  2440.    /* mandelfp */
  2441.  
  2442.    if(invert)
  2443.       invertz2(&init);
  2444.    else
  2445.       init.x = dx0[col]+dx1[row];
  2446.     switch (fractype)
  2447.       {
  2448.     case MAGNET2M:
  2449.         FloatPreCalcMagnet2();
  2450.     case MAGNET1M:         /* Crit Val Zero both, but neither   */
  2451.         old.x = old.y = 0.0; /* is of the form f(Z,C) = Z*g(Z)+C  */
  2452.         break;
  2453.     case MANDELLAMBDAFP:        /* Critical Value 0.5 + 0.0i  */
  2454.         old.x = 0.5;
  2455.         old.y = 0.0;
  2456.         break;
  2457.     default:
  2458.         old = init;
  2459.         break;
  2460.       }
  2461.  
  2462.    /* alter init value */
  2463.    if(useinitorbit == 1)
  2464.       old = initorbit;
  2465.    else if(useinitorbit == 2)
  2466.       old = init;
  2467.  
  2468.    if(inside == -60 || inside == -61)
  2469.    {
  2470.       /* kludge to match "Beauty of Fractals" picture since we start
  2471.      Mandelbrot iteration with init rather than 0 */
  2472.       old.x = parm.x; /* initial pertubation of parameters set */
  2473.       old.y = parm.y;
  2474.       coloriter = -1;
  2475.    }
  2476.    else
  2477.    {
  2478.      old.x += parm.x;
  2479.      old.y += parm.y;
  2480.    }
  2481.    tmp = init; /* for spider */
  2482.    tempsqrx = sqr(old.x);  /* precalculated value for regular Mandelbrot */
  2483.    tempsqry = sqr(old.y);
  2484.    return(1); /* 1st iteration has been done */
  2485. }
  2486.  
  2487. int juliafp_per_pixel()
  2488. {
  2489.    /* floating point julia */
  2490.    /* juliafp */
  2491.    if(invert)
  2492.       invertz2(&old);
  2493.    else
  2494.    {
  2495.      old.x = dx0[col]+dx1[row];
  2496.      old.y = dy0[row]+dy1[col];
  2497.    }
  2498.    tempsqrx = sqr(old.x);  /* precalculated value for regular Julia */
  2499.    tempsqry = sqr(old.y);
  2500.    tmp = old;
  2501.    return(0);
  2502. }
  2503.  
  2504. int MPCjulia_per_pixel()
  2505. {
  2506. #ifndef XFRACT
  2507.    /* floating point julia */
  2508.    /* juliafp */
  2509.    if(invert)
  2510.       invertz2(&old);
  2511.    else
  2512.    {
  2513.      old.x = dx0[col]+dx1[row];
  2514.      old.y = dy0[row]+dy1[col];
  2515.    }
  2516.    mpcold.x = *pd2MP(old.x);
  2517.    mpcold.y = *pd2MP(old.y);
  2518.    return(0);
  2519. #endif
  2520. }
  2521.  
  2522. otherrichard8fp_per_pixel()
  2523. {
  2524.     othermandelfp_per_pixel();
  2525.     CMPLXtrig1(*floatparm,tmp);
  2526.     CMPLXmult(tmp,parm2,tmp);
  2527.     return(1);
  2528. }
  2529.  
  2530. int othermandelfp_per_pixel()
  2531. {
  2532.    if(invert)
  2533.       invertz2(&init);
  2534.    else
  2535.       init.x = dx0[col]+dx1[row];
  2536.  
  2537.    if(useinitorbit == 1)
  2538.       old = initorbit;
  2539.    else
  2540.       old = init;
  2541.  
  2542.    old.x += parm.x;     /* initial pertubation of parameters set */
  2543.    old.y += parm.y;
  2544.  
  2545.    return(1); /* 1st iteration has been done */
  2546. }
  2547.  
  2548. int MPCHalley_per_pixel()
  2549. {
  2550. #ifndef XFRACT
  2551.    /* MPC halley */
  2552.    if(invert)
  2553.       invertz2(&init);
  2554.    else
  2555.      init.x = dx0[col]+dx1[row];
  2556.  
  2557.    mpcold.x = *pd2MP(init.x);
  2558.    mpcold.y = *pd2MP(init.y);
  2559.  
  2560.    return(0);
  2561. #endif
  2562. }
  2563.  
  2564. int Halley_per_pixel()
  2565. {
  2566.    if(invert)
  2567.       invertz2(&init);
  2568.    else
  2569.       init.x = dx0[col]+dx1[row];
  2570.  
  2571.    old = init;
  2572.  
  2573.    return(0); /* 1st iteration is not done */
  2574. }
  2575.  
  2576. int otherjuliafp_per_pixel()
  2577. {
  2578.    if(invert)
  2579.       invertz2(&old);
  2580.    else
  2581.    {
  2582.       old.x = dx0[col]+dx1[row];
  2583.       old.y = dy0[row]+dy1[col];
  2584.    }
  2585.    return(0);
  2586. }
  2587.  
  2588. #if 0
  2589. #define Q0 .113
  2590. #define Q1 .01
  2591. #else
  2592. #define Q0 0
  2593. #define Q1 0
  2594. #endif
  2595.  
  2596. int quaternionjulfp_per_pixel()
  2597. {
  2598.    old.x = dx0[col]+dx1[row];
  2599.    old.y = dy0[row]+dy1[col];
  2600.    floatparm->x = param[4];
  2601.    floatparm->y = param[5];
  2602.    qc  = param[0];
  2603.    qci = param[1];
  2604.    qcj = param[2];
  2605.    qck = param[3];
  2606.    return(0);
  2607. }
  2608.  
  2609. int quaternionfp_per_pixel()
  2610. {
  2611.    old.x = 0;
  2612.    old.y = 0;
  2613.    floatparm->x = 0;
  2614.    floatparm->y = 0;
  2615.    qc  = dx0[col]+dx1[row];
  2616.    qci = dy0[row]+dy1[col];
  2617.    qcj = param[2];
  2618.    qck = param[3];
  2619.    return(0);
  2620. }
  2621.  
  2622. int MarksCplxMandperp(void)
  2623. {
  2624.    if(invert)
  2625.       invertz2(&init);
  2626.    else
  2627.       init.x = dx0[col]+dx1[row];
  2628.    old.x = init.x + parm.x; /* initial pertubation of parameters set */
  2629.    old.y = init.y + parm.y;
  2630.    tempsqrx = sqr(old.x);  /* precalculated value */
  2631.    tempsqry = sqr(old.y);
  2632.    coefficient = ComplexPower(init, pwr);
  2633.    return(1);
  2634. }
  2635.  
  2636. int long_phoenix_per_pixel()
  2637. {
  2638. #ifndef XFRACT
  2639.    if(invert)
  2640.    {
  2641.       /* invert */
  2642.       invertz2(&old);
  2643.  
  2644.       /* watch out for overflow */
  2645.       if(sqr(old.x)+sqr(old.y) >= 127)
  2646.       {
  2647.      old.x = 8;  /* value to bail out in one iteration */
  2648.      old.y = 8;
  2649.       }
  2650.  
  2651.       /* convert to fudged longs */
  2652.       lold.x = (long)(old.x*fudge);
  2653.       lold.y = (long)(old.y*fudge);
  2654.    }
  2655.    else
  2656.    {
  2657.       lold.x = lx0[col]+lx1[row];
  2658.       lold.y = ly0[row]+ly1[col];
  2659.    }
  2660.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2661.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2662.    ltmp2.x = 0; /* use ltmp2 as the complex Y value */
  2663.    ltmp2.y = 0;
  2664.    return(0);
  2665. #else
  2666.    printf("Called long_phoenix_per_pixel\n");
  2667.    exit(0);
  2668. #endif
  2669. }
  2670.  
  2671. int phoenix_per_pixel()
  2672. {
  2673.    if(invert)
  2674.       invertz2(&old);
  2675.    else
  2676.    {
  2677.       old.x = dx0[col]+dx1[row];
  2678.       old.y = dy0[row]+dy1[col];
  2679.    }
  2680.    tempsqrx = sqr(old.x);  /* precalculated value */
  2681.    tempsqry = sqr(old.y);
  2682.    tmp2.x = 0; /* use tmp2 as the complex Y value */
  2683.    tmp2.y = 0;
  2684.    return(0);
  2685. }
  2686. int long_mandphoenix_per_pixel()
  2687. {
  2688. #ifndef XFRACT
  2689.    linit.x = lx0[col]+lx1[row];
  2690.  
  2691.    if(invert)
  2692.    {
  2693.       /* invert */
  2694.       invertz2(&init);
  2695.  
  2696.       /* watch out for overflow */
  2697.       if(sqr(init.x)+sqr(init.y) >= 127)
  2698.       {
  2699.      init.x = 8;  /* value to bail out in one iteration */
  2700.      init.y = 8;
  2701.       }
  2702.  
  2703.       /* convert to fudged longs */
  2704.       linit.x = (long)(init.x*fudge);
  2705.       linit.y = (long)(init.y*fudge);
  2706.    }
  2707.  
  2708.    if(useinitorbit == 1)
  2709.       lold = linitorbit;
  2710.    else
  2711.       lold = linit;
  2712.  
  2713.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2714.    lold.y += lparm.y;
  2715.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2716.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2717.    ltmp2.x = 0;
  2718.    ltmp2.y = 0;
  2719.    return(1); /* 1st iteration has been done */
  2720. #else
  2721.    printf("Called long_mandphoenix_per_pixel\n");
  2722.    exit(0);
  2723. #endif
  2724. }
  2725. int mandphoenix_per_pixel()
  2726. {
  2727.    if(invert)
  2728.       invertz2(&init);
  2729.    else
  2730.       init.x = dx0[col]+dx1[row];
  2731.  
  2732.    if(useinitorbit == 1)
  2733.       old = initorbit;
  2734.    else
  2735.       old = init;
  2736.  
  2737.    old.x += parm.x;     /* initial pertubation of parameters set */
  2738.    old.y += parm.y;
  2739.    tempsqrx = sqr(old.x);  /* precalculated value */
  2740.    tempsqry = sqr(old.y);
  2741.    tmp2.x = 0;
  2742.    tmp2.y = 0;
  2743.    return(1); /* 1st iteration has been done */
  2744. }
  2745.  
  2746. QuaternionFPFractal()
  2747. {
  2748.    double a0,a1,a2,a3,n0,n1,n2,n3;
  2749.    a0 = old.x;
  2750.    a1 = old.y;
  2751.    a2 = floatparm->x;
  2752.    a3 = floatparm->y;
  2753.  
  2754.    n0 = a0*a0-a1*a1-a2*a2-a3*a3 + qc;
  2755.    n1 = 2*a0*a1 + qci;
  2756.    n2 = 2*a0*a2 + qcj;
  2757.    n3 = 2*a0*a3 + qck;
  2758.    /* Check bailout */
  2759.    magnitude = a0*a0+a1*a1+a2*a2+a3*a3;
  2760.    if (magnitude>rqlim) {
  2761.        return 1;
  2762.    }
  2763.    old.x = n0;
  2764.    old.y = n1;
  2765.    floatparm->x = n2;
  2766.    floatparm->y = n3;
  2767.    return(0);
  2768. }
  2769.  
  2770. HyperComplexFPFractal()
  2771. {
  2772.    _HCMPLX hold, hnew;
  2773.    hold.x = old.x;
  2774.    hold.y = old.y;
  2775.    hold.z = floatparm->x;
  2776.    hold.t = floatparm->y;
  2777.  
  2778. /*   HComplexSqr(&hold,&hnew); */
  2779.    HComplexTrig0(&hold,&hnew);
  2780.    
  2781.    hnew.x += qc;
  2782.    hnew.y += qci;
  2783.    hnew.z += qcj;
  2784.    hnew.t += qck;
  2785.    
  2786.    old.x = hnew.x;
  2787.    old.y = hnew.y;
  2788.    floatparm->x = hnew.z;
  2789.    floatparm->y = hnew.t;
  2790.  
  2791.    /* Check bailout */
  2792.    magnitude = sqr(old.x)+sqr(old.y)+sqr(floatparm->x)+sqr(floatparm->y);
  2793.    if (magnitude>rqlim) {
  2794.        return 1;
  2795.    }
  2796.    return(0);
  2797. }
  2798.  
  2799. #if 0
  2800. demowalk()
  2801. {
  2802.     float stepsize;            /* average stepsize */
  2803.     int xwalk, ywalk;            /* current position */
  2804.     int xstep, ystep;            /* current step */
  2805.     long steps;                /* number of steps */
  2806.     int color;                /* color to draw this step */
  2807.     float temp, tempadjust;        /* temporary variables */
  2808.  
  2809. if (param[0] != 999) {            /* if 999, do a Mandelbrot instead */
  2810.  
  2811.     srand(rseed);            /* seed the random number generator */
  2812.     if (!rflag) ++rseed;
  2813.     tempadjust = RAND_MAX >> 2;        /* adjustment factor */
  2814.  
  2815.     xwalk = xdots / 2;            /* start in the center of the image */
  2816.     ywalk = ydots / 2;
  2817.  
  2818.     stepsize = min(xdots, ydots)     /* calculate average stepsize */
  2819.                * (param[0]/100.0);    /* as a percentage of the image */
  2820.  
  2821.     color = (int)max(0, min(colors, param[1]));  /* set the initial color */  
  2822.  
  2823.     for (steps = 0; steps < maxit; steps++) { /* take maxit steps */
  2824.         if (keypressed())        /* abort if told to do so */
  2825.             return(0);
  2826.         temp = rand();            /* calculate the next xstep */
  2827.         xstep = (int)(((temp/tempadjust) - 2.0) * stepsize);
  2828.         xstep = min(xwalk + xstep, xdots - 1);
  2829.         xstep = max(0, xstep);
  2830.         temp = rand();            /* calculate the next ystep */
  2831.         ystep = (int)(((temp/tempadjust) - 2.0) * stepsize);
  2832.         ystep = min(ywalk + ystep, ydots - 1);
  2833.         ystep = max(0, ystep);
  2834.         if (param[1] == 0.0)        /* rotate the colors? */
  2835.             if (++color >= colors)    /* rotate the colors, avoiding */
  2836.                 color = 1;        /* the background color 0 */
  2837.         /* the draw_line function is borrowed from the 3D routines */
  2838.         draw_line(xwalk, ywalk,xstep,ystep,color);
  2839.         /* or, we could be on a pogo stick and just displaying
  2840.            where we landed...
  2841.         putcolor(xstep, ystep, color);
  2842.         */
  2843.  
  2844.         xwalk = xstep;            /* remember where we were */
  2845.         ywalk = ystep;
  2846.         }
  2847.     return(1);                          /* we're done */
  2848.  
  2849. } else {                /* a simple Mandelbrot routine */
  2850.  
  2851.     /* the following routine determines the X and Y values of
  2852.        each pixel coordinate and calculates a simple mandelbrot
  2853.        fractal with them - slowly, but surely */
  2854.     int ix, iy;
  2855.     for (iy = 0; iy < ydots; iy++) {
  2856.         for (ix = 0; ix < xdots; ix++) {
  2857.             long iter;
  2858.             double x, y, newx, newy, tempxx, tempxy, tempyy;
  2859.             /* first, obtain the X and Y coordinate values of this pixel */
  2860.             x = dx0[ix]+dx1[iy];
  2861.             y = dy0[iy]+dy1[ix];
  2862.             /* now initialize the temporary values */
  2863.             tempxx = tempyy = tempxy = 0.0;
  2864.             if (keypressed())        /* abort if told to do so */
  2865.                 return(0);
  2866.             /* the inner iteration loop */
  2867.             for (iter = 1; iter < maxit; iter++) {
  2868.                 /* calculate the X and Y values of Z(iter) */
  2869.                 newx = tempxx - tempyy + x;
  2870.                 newy = tempxy + tempxy + y;
  2871.                 /* calculate the temporary values */
  2872.                 tempxx = newx * newx;
  2873.                 tempyy = newy * newy;
  2874.                 tempxy = newx * newy;
  2875.                 /* are we done yet? */
  2876.                 if (tempxx + tempyy > 4.0) break;
  2877.                 }
  2878.             /* color in the pixel */
  2879.             putcolor(ix, iy, (int)iter % colors);
  2880.             }
  2881.         }
  2882.     return(1);                          /* we're done */
  2883.     }
  2884.  
  2885. }
  2886. #endif
  2887.