home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / FORTH / FRASRC11.ZIP / FRACTALS.C < prev    next >
C/C++ Source or Header  |  1989-12-31  |  93KB  |  3,053 lines
  1. /*
  2. FRACTALS.C and CALCFRAC.C actually calculate the fractal images (well,
  3. SOMEBODY had to do it!).  The modules are set up so that all logic that
  4. is independent of any fractal-specific code is in CALCFRAC.C, and the
  5. code that IS fractal-specific is in FRACTALS.C. Original author Tim Wegner,
  6. but just about ALL the authors have contributed SOME code to this routine
  7. at one time or another, or contributed to one of the many massive
  8. restructurings.
  9.  
  10. The Fractal-specific routines are divided into three categories:
  11.  
  12. 1. Routines that are called once-per-orbit to calculate the orbit
  13.    value. These have names live "XxxxFractal", and their function
  14.    pointers are stored in fractalspecific[fractype].orbit_calc. EVERY 
  15.    new fractal type needs one of these. Return 0 to continue iterations,
  16.    1 if we're done. Results for integer fractals are left in 'lnew.x' and 
  17.    'lnew.y', for floating point fractals in 'new.x' and 'new.y'.
  18.    
  19. 2. Routines that are called once per pixel to set various variables
  20.    prior to the orbit calculation. These have names like xxx_per_pixel
  21.    and are fairly generic - chances are one is right for your new type.
  22.    They are stored in fractalspecific[fractype].per_pixel.
  23.    
  24. 3. Routines that are called once per screen to set various variables.
  25.    These have names like XxxxSetup, and are stored in 
  26.    fractalspecific[fractype].per_image.
  27.    
  28. 4. The main fractal routine. Usually this will be StandardFractal(),
  29.    but if you have written a stand-alone fractal routine independent
  30.    of the StandardFractal mechanisms, your routine name goes here, 
  31.    stored in fractalspecific[fractype].calctype.per_image.
  32.  
  33. Adding a new fractal type should be simply a matter of adding an item
  34. to the 'fractalspecific' structure, writing (or re-using one of the existing) 
  35. an appropriate setup, per_image, per_pixel, and orbit routines.
  36.  
  37. --------------------------------------------------------------------   */
  38.  
  39. /* -------------------------------------------------------------------- */
  40. /*    These variables are not specific to any particular fractal    */
  41. /* -------------------------------------------------------------------- */
  42.  
  43. #include <stdio.h>
  44. #include <stdlib.h>
  45. #include <float.h>
  46. #include <limits.h>
  47. #include "fractint.h"
  48. #include "mpmath.h"
  49.  
  50. #ifndef __TURBOC__
  51. #include <malloc.h>
  52. #endif
  53.  
  54. /* defines should match fractalspecific array indices */
  55.  
  56. /* first values defined in fractype.h for historical reasons  - from there:
  57.   #define NOFRACTAL      -1
  58.   #define MANDEL          0 
  59.   #define JULIA           1 
  60.   #define NEWTBASIN       2 
  61.   #define LAMBDA          3 
  62.   #define MANDELFP        4 
  63.   #define NEWTON          5 
  64.   #define JULIAFP         6 
  65.   #define PLASMA          7 
  66. */
  67.  
  68. #define LAMBDASINE        8 
  69. #define LAMBDACOS         9 
  70. #define LAMBDAEXP         10
  71. #define TEST              11
  72. #define SIERPINSKI        12
  73. #define BARNSLEYM1        13
  74. #define BARNSLEYJ1        14
  75. #define BARNSLEYM2        15
  76. #define BARNSLEYJ2        16
  77. #define MANDELSINE        17
  78. #define MANDELCOS         18
  79. #define MANDELEXP         19
  80. #define MANDELLAMBDA      20
  81. #define MARKSMANDEL       21
  82. #define MARKSJULIA        22
  83. #define UNITY             23
  84. #define MANDEL4           24
  85. #define JULIA4            25
  86. #define IFS               26
  87. #define IFS3D             27
  88. #define BARNSLEYM3        28
  89. #define BARNSLEYJ3        29
  90. #define DEMM              30
  91. #define DEMJ              31
  92. #define BIFURCATION       32
  93. #define MANDELSINH        33
  94. #define LAMBDASINH        34
  95. #define MANDELCOSH        35
  96. #define LAMBDACOSH        36
  97. #define LMANDELSINE       37
  98. #define LLAMBDASINE       38
  99. #define LMANDELCOS        39
  100. #define LLAMBDACOS        40
  101. #define LMANDELSINH       41
  102. #define LLAMBDASINH       42
  103. #define LMANDELCOSH       43
  104. #define LLAMBDACOSH       44
  105. #define LMANSINZSQRD      45
  106. #define LJULSINZSQRD      46
  107. #define FPMANSINZSQRD     47
  108. #define FPJULSINZSQRD     48
  109. #define LMANDELEXP        49
  110. #define LLAMBDAEXP        50
  111. #define LMANDELZPOWER     51
  112. #define LJULIAZPOWER      52
  113. #define FPMANDELZPOWER    53
  114. #define FPJULIAZPOWER     54
  115. #define FPMANZTOZPLUSZPWR 55
  116. #define FPJULZTOZPLUSZPWR 56
  117. #define LMANSINEXP        57
  118. #define LJULSINEXP        58
  119. #define FPMANSINEXP       59
  120. #define FPJULSINEXP       60
  121. #define FPPOPCORN         61
  122. #define LPOPCORN          62
  123. #define FPLORENZ          63
  124. #define LLORENZ           64
  125. #define LORENZ3D          65
  126. #define MPNEWTON          66 
  127. #define MPNEWTBASIN       67 
  128. #define COMPLEXNEWTON     68
  129.  
  130. /* DEFINITIONS PRIOR TO THIS POINT ARE FROZEN BY THE VERSION 11.0 RELEASE! */
  131.  
  132. #define NEWTONDEGREELIMIT  100
  133. extern void draw_line(int,int,int,int,int);
  134. long Exp086(long);
  135. double fmod(double,double);
  136. extern int biomorph;
  137. extern int forcesymmetry;
  138. struct lcomplex lcoefficient,lold,lnew,llambda, linit,ltemp;
  139. long ltempsqrx,ltempsqry;
  140. extern int decomp[];
  141. extern double param[];
  142. extern double inversion[];                     /* inversion parameters */
  143. extern double f_radius,f_xcenter,f_ycenter;    /* inversion radius, center */
  144. extern double xxmin,xxmax,yymin,yymax;         /* corners */
  145. extern VECTOR view;
  146. extern int init3d[];
  147. extern int overflow;
  148.  
  149. int maxcolor;
  150. int root, degree,basin;
  151. double floatmin,floatmax;
  152. double roverd, d1overd, threshold;
  153. long lroverd,ld1overd, lthreshold;
  154. int switchedtofloat;
  155. extern struct complex init,tmp,old,new,saved;
  156. struct complex  staticroots[16]; /* roots array for degree 16 or less */
  157. struct complex  *roots = staticroots;
  158. struct MPC      *MPCroots;
  159. extern int color, oldcolor, oldmax, oldmax10, row, col, passes;
  160. extern int iterations, invert;
  161. extern double far *dx0, far *dy0;
  162. extern long XXOne, FgOne, FgTwo, LowerLimit;
  163. struct complex one;
  164.  
  165. extern void (*plot)();
  166.  
  167. extern int    xdots, ydots;        /* coordinates of dots on the screen  */
  168. extern int    colors;                /* maximum colors available */
  169. extern int    inside;                /* "inside" color to use    */
  170. extern int    maxit;                /* try this many iterations */
  171. extern int    fractype;            /* fractal type */
  172. extern int    numpasses;            /* 0 = 1 pass, 1 = double pass */
  173. extern int    solidguessing;        /* 1 if solid-guessing */
  174. extern int    debugflag;            /* for debugging purposes */
  175. extern int    timerflag;            /* ditto */
  176. extern    int    diskvideo;            /* for disk-video klooges */
  177.  
  178. extern double    param[];        /* parameters */
  179. extern double    potparam[];        /* potential parameters */
  180. extern long    lx0[], ly0[];        /* X, Y points */
  181. extern long    delx,dely;            /* X, Y increments */
  182. extern long    fudge;                /* fudge factor (2**n) */
  183. extern int    bitshift;            /* bit shift for fudge */
  184. extern char potfile[];          /* potential filename */
  185.  
  186. #ifndef sqr
  187. #define sqr(x) ((x)*(x))
  188. #endif
  189.  
  190. #ifndef lsqr
  191. #define lsqr(x) (multiply((x),(x),bitshift))
  192. #endif
  193.  
  194. #define modulus(z)       (sqr((z).x)+sqr((z).y))
  195. #define conjugate(pz)   ((pz)->y = 0.0 - (pz)->y)
  196. #define distance(z1,z2)  (sqr((z1).x-(z2).x)+sqr((z1).y-(z2).y))
  197. #define pMPsqr(z) (pMPmul((z),(z)))
  198. #define MPdistance(z1,z2)  (pMPadd(pMPsqr(pMPsub((z1).r,(z2).r)),pMPsqr(pMPsub((z1).i,(z2).i))))
  199.  
  200. double twopi = PI*2.0;
  201. static int c_exp;
  202.  
  203. /* These are local but I don't want to pass them as parameters */
  204. extern struct complex lambda;
  205. extern double deltaX, deltaY;
  206. extern double magnitude, rqlim, rqlim2;
  207. extern int XXdots, YYdots; /* local dots */
  208. struct complex parm;
  209. struct complex *floatparm;
  210. struct lcomplex *longparm;
  211. extern int (*calctype)();
  212. extern double closenuff;
  213. extern int pixelpi; /* value of pi in pixels */
  214. extern FILE *fp_pot;
  215. extern int potflag; /* tells if continuous potential on  */
  216. extern unsigned long lm;        /* magnitude limit (CALCMAND) */
  217.  
  218. extern int    ixstart, ixstop, iystart, iystop;    /* (for CALCMAND) start, stop here */
  219.  
  220. /* -------------------------------------------------------------------- */
  221. /*        These variables are external for speed's sake only    */
  222. /* -------------------------------------------------------------------- */
  223.  
  224. double sinx,cosx,sinhx,coshx;
  225. double siny,cosy,sinhy,coshy;
  226. double tmpexp;
  227. double tempsqrx,tempsqry;
  228.  
  229. long oldxinitx,oldyinity,oldxinity,oldyinitx;
  230. long modulusinit;
  231. long ltmp1,ltmp2,ltmp3;
  232. struct lcomplex ltmp;
  233. extern long lmagnitud, llimit, llimit2, lclosenuff, l16triglim;
  234. extern periodicitycheck;
  235. extern char floatflag;
  236.  
  237. extern int StandardFractal();
  238. extern int NewtonFractal2(); /* Lee Crocker's Newton code */
  239.  
  240. /* these are in mpmath_c.c */
  241. extern int ComplexNewtonSetup(void);
  242. extern int ComplexNewton(void), ComplexBasin(void);
  243. complex_mult(struct complex arg1,struct complex arg2,struct complex *pz);
  244. complex_div(struct complex arg1,struct complex arg2,struct complex *pz);
  245.  
  246. /* -------------------------------------------------------------------- */
  247. /*        Stand-alone routines                                            */
  248. /* -------------------------------------------------------------------- */
  249.  
  250. extern char far plasmamessage[];
  251.  
  252. /* Thanks to Rob Beyer */
  253. floatlorenz() 
  254. {
  255.    int count;
  256.    double dx,dy,dz,x,y,z,dt,a,b,c;
  257.    double adt,bdt,cdt,xdt,ydt;
  258.    int oldrow, oldcol;
  259.    count = 0;
  260.    if(inside > 0)
  261.       color = inside;
  262.    if(color >= colors)
  263.       color = 1;   
  264.    oldcol = oldrow = -1;
  265.    x = 1;  /* initial conditions */
  266.    y = 1;
  267.    z = 1;
  268.    if(param[0] > 0 && param[0] <= .05)
  269.      dt = param[0];
  270.    else
  271.      dt = .02; /* time step */
  272.    if(param[1])
  273.       a = param[1];
  274.    else
  275.       a = 5;
  276.    if(param[2])
  277.       b = param[2];
  278.    else
  279.       b = 15;
  280.    if(param[3])
  281.       c = param[3];
  282.    else
  283.       c = 1;
  284.          
  285.    /* precalculations for speed */
  286.    adt = a*dt;
  287.    bdt = b*dt;
  288.    cdt = c*dt;
  289.    
  290.    while(1)
  291.    {
  292.       if (++count > 1000) 
  293.       {        /* time to switch colors? */
  294.          count = 0;
  295.          if (++color >= colors)   /* another color to switch to? */
  296.               color = 1;        /* (don't use the background color) */
  297.       }
  298.       if(check_key())
  299.          return(-1);
  300.       if ( x > xxmin && x < xxmax && z > yymin && z < yymax )
  301.       {
  302.          col =          (( x-xxmin) / deltaX);
  303.          row = YYdots - (( z-yymin) / deltaY);
  304.          if(oldcol != -1)
  305.             draw_line(col,row,oldcol,oldrow,color&(colors-1));
  306.          else            
  307.             (*plot)(col,row,color&(colors-1));
  308.          oldcol = col;
  309.          oldrow = row;    
  310.       }
  311.       else
  312.          oldrow = oldcol = -1;
  313.  
  314.       /* Calculate the next point */
  315.       xdt = x*dt;
  316.       ydt = y*dt;
  317.       dx = -(adt * x) + (adt * y);
  318.       dy =  (bdt * x) - (ydt) - (z * xdt);
  319.       dz = -(cdt * z) + (x * ydt);
  320.  
  321.       x += dx;
  322.       y += dy;
  323.       z += dz;
  324.    }
  325.    return(0); 
  326. }
  327. longlorenz() 
  328. {
  329.    int count;
  330.    long delx,dely,xmin,xmax,ymin,ymax;
  331.    long dx,dy,dz,x,y,z,dt,a,b,c;
  332.    long adt,bdt,cdt,xdt,ydt;
  333.    int oldrow, oldcol;
  334.    count = 0;
  335.    if(inside > 0)
  336.       color = inside;
  337.    if(color >= colors)
  338.       color = 1;   
  339.    oldcol = oldrow = -1;
  340.    fudge = 1L<<bitshift;
  341.    delx = deltaX*fudge;
  342.    dely = deltaY*fudge;
  343.    xmin = xxmin*fudge;
  344.    ymin = yymin*fudge;
  345.    xmax = xxmax*fudge;
  346.    ymax = yymax*fudge;
  347.    x = fudge;  /* initial conditions */
  348.    y = fudge;
  349.    z = fudge;
  350.    if(param[0] > 0 && param[0] <= .05)
  351.      dt = param[0]*fudge;
  352.    else
  353.      dt = .02*fudge; /* time step */
  354.    if(param[1])
  355.       a = param[1];
  356.    else
  357.       a = 5;
  358.    if(param[2])
  359.       b = param[2];
  360.    else
  361.       b = 15;
  362.    if(param[3])
  363.       c = param[3];
  364.    else
  365.       c = 1;
  366.    
  367.    /* precalculations for speed */
  368.    adt = a*dt;
  369.    bdt = b*dt;
  370.    cdt = c*dt;
  371.    while(1)
  372.    {
  373.       if (++count > 1000) 
  374.       {        /* time to switch colors? */
  375.          count = 0;
  376.          if (++color >= colors)   /* another color to switch to? */
  377.               color = 1;        /* (don't use the background color) */
  378.       }
  379.       if(check_key())
  380.          return(-1);
  381.       if ( x > xmin && x < xmax && z > ymin && z < ymax )
  382.       {
  383.          col =          (( x-xmin) / delx);
  384.          row = YYdots - (( z-ymin) / dely);
  385.          if(oldcol != -1)
  386.             draw_line(col,row,oldcol,oldrow,color&(colors-1));
  387.          else            
  388.             (*plot)(col,row,color&(colors-1));
  389.          oldcol = col;
  390.          oldrow = row;    
  391.       }
  392.       else
  393.          oldrow = oldcol = -1;
  394.  
  395.       /* Calculate the next point */
  396.       
  397.       xdt = multiply(x,dt,bitshift);
  398.       ydt = multiply(y,dt,bitshift);
  399.       dx  = -multiply(adt,x,bitshift) + multiply(adt,y,bitshift);
  400.       dy  =  multiply(bdt,x,bitshift) -ydt -multiply(z,xdt,bitshift);
  401.       dz  = -multiply(cdt,z,bitshift) + multiply(x,ydt,bitshift);
  402.  
  403.       x += dx;
  404.       y += dy;
  405.       z += dz;
  406.    }
  407.    return(0); 
  408. }
  409.  
  410. /* this ought to be combined with ifs3d */
  411. lorenz3dlong()
  412. {
  413.    unsigned count;
  414.    long dx,dy,dz,x,y,z,dt,a,b,c;
  415.    long adt,bdt,cdt,xdt,ydt;
  416.    int oldcol,oldrow;
  417.    long delx,dely;
  418.       
  419.    double tmpx, tmpy, tmpz;
  420.    extern VECTOR view;
  421.    long iview[3];         /* perspective viewer's coordinates */      
  422.    long tmp;   
  423.    long maxifs[3];
  424.    long minifs[3];
  425.    extern int init3d[];
  426.    unsigned long maxct,ct;
  427.    register int col;
  428.    register int row;
  429.    register int color;
  430.  
  431.    long localifs[NUMIFS][IFS3DPARM];        /* local IFS values */
  432.    long newx,newy,newz,r,sum, xmin, xmax, ymin, ymax;
  433.     
  434.    int i,j,k;
  435.  
  436.    /* 3D viewing stuff */
  437.    MATRIX doublemat;           /* transformation matrix */
  438.    long longmat[4][4];         /* long version of matrix */
  439.    long ifsvect[3];               /* interated function orbit value */
  440.    long viewvect[3];        /* orbit transformed for viewing */
  441.  
  442.    /*
  443.    printf("xrot %d yrot %d zrot %d \nxshift %d yshift %d perspective %d\n",
  444.              XROT,YROT,ZROT,XSHIFT,YSHIFT,ZVIEWER);
  445.    */
  446.    oldcol = oldrow = -1;
  447.    color = 2;
  448.    fudge = 1L << bitshift;
  449.  
  450.    delx = deltaX*fudge;
  451.    dely = deltaY*fudge;
  452.       
  453.    for(i=0;i<3;i++)
  454.    {
  455.       minifs[i] =  1L << 30;
  456.       maxifs[i] = -minifs[i];
  457.    }
  458.    
  459.    /* build transformation matrix */
  460.    identity (doublemat);
  461.  
  462.    /* apply rotations - uses the same rotation variables as line3d.c */
  463.    xrot ((double)XROT / 57.29577,doublemat);
  464.    yrot ((double)YROT / 57.29577,doublemat);
  465.    zrot ((double)ZROT / 57.29577,doublemat);
  466.  
  467.    /* apply scale */
  468.    scale((double)XSCALE/100.0,(double)YSCALE/100.0,(double)ROUGH/100.0,doublemat);
  469.  
  470.    if(!ZVIEWER)
  471.       trans((double)XSHIFT*(xxmax-xxmin)/100.0,(double)YSHIFT*(yymax-yymin)/100.0,0.0,doublemat);
  472.  
  473.    /* copy xform matrix to long for for fixed point math */
  474.    for (i = 0; i < 4; i++)
  475.    for (j = 0; j < 4; j++)
  476.       longmat[i][j] = doublemat[i][j] * fudge;
  477.  
  478.    if(diskvideo)        /* this would KILL a disk drive! */
  479.    {
  480.       setvideomode(3,0,0,0);
  481.       buzzer(2);
  482.       helpmessage(plasmamessage);
  483.       return(-1);
  484.    }
  485.  
  486.    for (i = 0; i < NUMIFS; i++)    /* fill in the local IFS array */
  487.    for (j = 0; j < IFS3DPARM; j++)
  488.          localifs[i][j] = initifs3d[i][j] * fudge;
  489.  
  490.    xmin  = xxmin*fudge;   
  491.    xmax  = xxmax*fudge;
  492.    ymax  = yymax*fudge;
  493.    ymin  = yymin*fudge;
  494.  
  495.    x = fudge;  /* initial conditions */
  496.    y = fudge;
  497.    z = fudge;
  498.    if(param[0] > 0 && param[0] <= .05)
  499.      dt = param[0]*fudge;
  500.    else
  501.      dt = .02*fudge; /* time step */
  502.    if(param[1])
  503.       a = param[1];
  504.    else
  505.       a = 5;
  506.    if(param[2])
  507.       b = param[2];
  508.    else
  509.       b = 15;
  510.    if(param[3])
  511.       c = param[3];
  512.    else
  513.       c = 1;
  514.    
  515.    /* precalculations for speed */
  516.    adt = a*dt;
  517.    bdt = b*dt;
  518.    cdt = c*dt;
  519.  
  520.    ifsvect[0] = x;
  521.    ifsvect[1] = y;
  522.    ifsvect[2] = z;
  523.    
  524.    /* make maxct a function of screen size               */
  525.    /* 1k times maxit at EGA resolution seems about right */
  526.    maxct = (float)maxit*(1024.0*xdots*ydots)/(640.0*350.0); 
  527.    ct = 0L;
  528.    while(ct++ < maxct) /* loop until keypress or maxit */ 
  529.    {
  530.       if (++count > 1000) 
  531.       {        /* time to switch colors? */
  532.          count = 0;
  533.          if (++color >= colors)   /* another color to switch to? */
  534.               color = 1;        /* (don't use the background color) */
  535.       }
  536.       if(check_key())
  537.          return(-1);
  538.  
  539.       x = ifsvect[0];
  540.       y = ifsvect[1];
  541.       z = ifsvect[2];
  542.  
  543.       xdt = multiply(x,dt,bitshift);
  544.       ydt = multiply(y,dt,bitshift);
  545.       dx  = -multiply(adt,x,bitshift) + multiply(adt,y,bitshift);
  546.       dy  =  multiply(bdt,x,bitshift) -ydt -multiply(z,xdt,bitshift);
  547.       dz  = -multiply(cdt,z,bitshift) + multiply(x,ydt,bitshift);
  548.  
  549.       x += dx;
  550.       y += dy;
  551.       z += dz;
  552.  
  553.       ifsvect[0] = x;
  554.       ifsvect[1] = y;
  555.       ifsvect[2] = z;
  556.  
  557.  
  558.       /* 3D VIEWING TRANSFORM */
  559.       longvmult(ifsvect,longmat,viewvect, bitshift);
  560.  
  561.       /* find minz and maxz */
  562.       if( ZVIEWER && ct <= 100) /* waste this many points to find minz and maxz */
  563.       {
  564.          for(i=0;i<3;i++)
  565.             if ((tmp = viewvect[i]) < minifs[i])
  566.                minifs[i] = tmp;
  567.             else if (tmp > maxifs[i])
  568.                maxifs[i] = tmp;
  569.          if(ct == 100)
  570.          {
  571.             /* set of view vector */
  572.             for(i=0;i<2;i++)
  573.                iview[i] = (maxifs[i]+minifs[i])/2;
  574.       
  575.             /* z value of user's eye - should be more negative than extreme 
  576.                negative part of image */
  577.             iview[2] = (long)((minifs[2]-maxifs[2])*(double)ZVIEWER/100.0);
  578.  
  579.             /*
  580.             printf("min %ld max %ld iview[2] %d\n",minifs[2],maxifs[2],iview[2]);
  581.             getch();
  582.             */ 
  583.  
  584.             /* translate back exactly amount so maximum values are non-positive */
  585.             tmpx = ((double)XSHIFT*(xmax-xmin))/(100.0*fudge);
  586.             tmpy = ((double)YSHIFT*(ymax-ymin))/(100.0*fudge);
  587.             tmpz = -((double)maxifs[2]);
  588.             tmpz *= 1.1;
  589.             tmpz /= fudge; 
  590.             trans(tmpx,tmpy,tmpz,doublemat); 
  591.             
  592.             for(i=0;i<3;i++)
  593.             {
  594.                view[i] = iview[i];
  595.                view[i] /= fudge;
  596.             } 
  597.    
  598.             /* copy xform matrix to long for for fixed point math */
  599.             for (i = 0; i < 4; i++)
  600.             for (j = 0; j < 4; j++)
  601.                longmat[i][j] = doublemat[i][j] * fudge;
  602.          }
  603.       }      
  604.       
  605.       /* apply perspective if requested */
  606.       if(ZVIEWER && ct > 100)
  607.       {
  608.          if(debugflag==22 || ZVIEWER < 100) /* use float for small persp */
  609.          {
  610.             /* use float perspective calc */
  611.             VECTOR tmpv;
  612.             for(i=0;i<3;i++)
  613.             {
  614.                tmpv[i] = viewvect[i];
  615.                tmpv[i] /= fudge;
  616.             } 
  617.             perspective(tmpv);
  618.             for(i=0;i<3;i++)
  619.                viewvect[i] = tmpv[i]*fudge;
  620.          }
  621.          else
  622.             longpersp(viewvect,iview,bitshift);
  623.       }
  624.             
  625.       /* plot if inside window */
  626.       if ( viewvect[0] > xmin && viewvect[0] < xmax 
  627.         && viewvect[1] > ymin && viewvect[1] < ymax && ct > 100 )
  628.       {
  629.          col =          (( viewvect[0]-xmin) / delx);
  630.          row = YYdots - (( viewvect[1]-ymin) / dely);
  631.          if(oldcol != -1)
  632.             draw_line(col,row,oldcol,oldrow,color&(colors-1));
  633.          else            
  634.             (*plot)(col,row,color&(colors-1));
  635.          oldcol = col;
  636.          oldrow = row;    
  637.       }
  638.       else
  639.          oldrow = oldcol = -1;
  640.    }
  641.    return(0);
  642. }
  643.  
  644. ifs()            /* IFS logic shamelessly converted to integer math */
  645. {
  646.    long     *lifsptr;
  647.    unsigned long maxct,ct;
  648.    register int col;
  649.    register int row;
  650.    register int color;
  651.  
  652.    long localifs[NUMIFS][7];        /* local IFS values */
  653.    long x,y,newx,newy,r,sum, xmin, xmax, ymin, ymax, tempr;
  654.  
  655.    int i,j,k;
  656.    
  657.    if(diskvideo) {        /* this would KILL a disk drive! */
  658.     setvideomode(3,0,0,0);
  659.     buzzer(2);
  660.     helpmessage(plasmamessage);
  661.     return(-1);
  662.     }
  663.  
  664.    for (i = 0; i < NUMIFS; i++)    /* fill in the local IFS array */
  665.    for (j = 0; j < IFSPARM; j++)
  666.          localifs[i][j] = initifs[i][j] * fudge;
  667.  
  668.    xmin  = lx0[0];              /* find the screen corners */
  669.    xmax  = lx0[xdots-1];
  670.    ymax  = ly0[0];
  671.    ymin  = ly0[ydots-1];
  672.  
  673.    tempr = fudge / 32767;        /* find the proper rand() fudge */
  674.  
  675.    x = 0;
  676.    y = 0;
  677.  
  678.    /* make maxct a function of screen size               */
  679.    /* 1k times maxit at EGA resolution seems about right */
  680.    maxct = (float)maxit*(1024.0*xdots*ydots)/(640.0*350.0); 
  681.    ct = 0L;
  682.    while(ct++ < maxct) /* loop until keypress or maxit */ 
  683.    {
  684.       if (ct & 127)           /* reduce the number of keypress checks */
  685.          if( check_key() )    /* keypress bails out */
  686.             return(-1);
  687.       r = rand();        /* generate fudged random number between 0 and 1 */
  688.       r *= tempr;
  689.  
  690.       /* pick which iterated function to execute, weighted by probability */
  691.       sum = localifs[0][6];
  692.       k = 0;
  693.       while ( sum < r)
  694.       {
  695.          k++;
  696.          sum += localifs[k][6];
  697.          if (localifs[k+1][6] == 0) break;    /* for safety */
  698.       }
  699.  
  700.       /* calculate image of last point under selected iterated function */
  701.       newx = multiply(localifs[k][0], x, bitshift) +
  702.           multiply(localifs[k][1], y, bitshift) +
  703.           localifs[k][4];
  704.       newy = multiply(localifs[k][2], x, bitshift) +
  705.           multiply(localifs[k][3], y, bitshift) +
  706.           localifs[k][5];
  707.       x = newx;
  708.       y = newy;
  709.  
  710.       /* plot if inside window */
  711.       if ( x > xmin && x < xmax && y > ymin && y < ymax )
  712.       {
  713.          col =          (( x-xmin) / delx);
  714.          row = YYdots - (( y-ymin) / dely);
  715.  
  716.          /* color is count of hits on this pixel */
  717.          color = getcolor(col,row)+1;
  718.          if( color < colors ) /* color sticks on last value */
  719.             (*plot)(col,row,color);
  720.       }
  721.    }
  722.    return(0);
  723. }
  724.  
  725. ifs3d()
  726. {
  727.    if(floatflag)
  728.       return(ifs3dfloat()); /* double version of ifs3d */
  729.    else
  730.       return(ifs3dlong());  /* long version of ifs3d   */
  731. }
  732.  
  733. /* double version - mainly for testing */
  734. ifs3dfloat()
  735. {
  736.    double tmp;   
  737.    VECTOR minifs;
  738.    VECTOR maxifs;
  739.    unsigned long maxct,ct;
  740.    register int col;
  741.    register int row;
  742.    register int color;
  743.  
  744.    double newx,newy,newz,r,sum;
  745.     
  746.    int i,k;
  747.  
  748.    /* 3D viewing stuff */
  749.    MATRIX doublemat;           /* transformation matrix */
  750.    VECTOR ifsvect;               /* interated function orbit value */
  751.    VECTOR viewvect;        /* orbit transformed for viewing */
  752.    /*
  753.    printf("xrot %d yrot %d zrot %d \nxshift %d yshift %d perspective %d\n",
  754.              XROT,YROT,ZROT,XSHIFT,YSHIFT,ZVIEWER);
  755.    */
  756.    for(i=0;i<3;i++)
  757.    {
  758.       minifs[i] =  100000.0; /* impossible value */
  759.       maxifs[i] = -100000.0;
  760.    }
  761.    
  762.    /* build transformation matrix */
  763.    identity (doublemat);
  764.  
  765.    /* apply rotations - uses the same rotation variables as line3d.c */
  766.    xrot ((double)XROT / 57.29577,doublemat);
  767.    yrot ((double)YROT / 57.29577,doublemat);
  768.    zrot ((double)ZROT / 57.29577,doublemat);
  769.    /*
  770.    scale((double)XSCALE/100.0,(double)YSCALE/100.0,(double)ROUGH/100.0,doublemat);
  771.    */
  772.    if(!ZVIEWER)
  773.       trans((double)XSHIFT*(xxmax-xxmin)/100.0,(double)YSHIFT*(yymax-yymin)/100.0,0.0,doublemat);
  774.  
  775.    if(diskvideo)        /* this would KILL a disk drive! */
  776.    {
  777.       setvideomode(3,0,0,0);
  778.       buzzer(2);
  779.       helpmessage(plasmamessage);
  780.       return(-1);
  781.    }
  782.  
  783.    ifsvect[0] = 0;
  784.    ifsvect[1] = 0;
  785.    ifsvect[2] = 0;
  786.    
  787.    /* make maxct a function of screen size               */
  788.    /* 1k times maxit at EGA resolution seems about right */
  789.    maxct = (double)maxit*(1024.0*xdots*ydots)/(640.0*350.0); 
  790.    ct = 0L;
  791.    while(ct++ < maxct) /* loop until keypress or maxit */ 
  792.    {
  793.       if (ct & 127)           /* reduce the number of keypress checks */
  794.          if( check_key() )    /* keypress bails out */
  795.             return(-1);
  796.       r = rand();        /* generate fudged random number between 0 and 1 */
  797.       r /= 32767;
  798.  
  799.       /* pick which iterated function to execute, weighted by probability */
  800.       sum = initifs3d[0][12];
  801.       k = 0;
  802.       while ( sum < r)
  803.       {
  804.          k++;
  805.          sum += initifs3d[k][12];
  806.          if (initifs3d[k+1][12] == 0) break;    /* for safety */
  807.       }
  808.  
  809.       /* calculate image of last point under selected iterated function */
  810.       newx = initifs3d[k][0] * ifsvect[0] +
  811.              initifs3d[k][1] * ifsvect[1] +
  812.              initifs3d[k][2] * ifsvect[2] + initifs3d[k][9];
  813.       newy = initifs3d[k][3] * ifsvect[0] +
  814.              initifs3d[k][4] * ifsvect[1] +
  815.              initifs3d[k][5] * ifsvect[2] + initifs3d[k][10];
  816.       newz = initifs3d[k][6] * ifsvect[0] +
  817.              initifs3d[k][7] * ifsvect[1] +
  818.              initifs3d[k][8] * ifsvect[2] + initifs3d[k][11];
  819.             
  820.       ifsvect[0] = newx;
  821.       ifsvect[1] = newy;
  822.       ifsvect[2] = newz;
  823.  
  824.       /* 3D VIEWING TRANSFORM */
  825.       vmult(ifsvect,doublemat,viewvect);
  826.  
  827.       /* find minz and maxz */
  828.       if(ZVIEWER && ct <= 100) /* waste this many points to find minz and maxz */
  829.       {
  830.          for(i=0;i<3;i++)
  831.             if ((tmp = viewvect[i]) < minifs[i])
  832.                minifs[i] = tmp;
  833.             else if (tmp > maxifs[i])
  834.                maxifs[i] = tmp;
  835.          if(ct == 100)
  836.          {
  837.             /* set of view vector */
  838.             for(i=0;i<2;i++)
  839.                view[i] = (maxifs[i]+minifs[i])/2.0;
  840.       
  841.             /* z value of user's eye - should be more negative than extreme 
  842.                negative part of image */
  843.             view[2] = (minifs[2]-maxifs[2])*(double)ZVIEWER/100.0;
  844.  
  845.             /* translate back exactly amount so maximum values are non-positive */
  846.             trans((double)XSHIFT*(xxmax-xxmin)/100.0,(double)YSHIFT*(yymax-yymin)/100.0,-maxifs[2],doublemat); 
  847. /*
  848.             printf("minx %lf maxx %lf\n",minifs[0],maxifs[0]);
  849.             printf("miny %lf maxy %lf\n",minifs[1],maxifs[1]);
  850.             printf("minz %lf maxz %lf\n",minifs[2],maxifs[2]);
  851.             printf("view %lf %lf %lf\n",view[0],view[1],view[2]); */
  852.          }
  853.       }      
  854.       
  855.       /* apply perspective if requested */
  856.       if(ZVIEWER)
  857.          perspective(viewvect);
  858.  
  859.       /* plot if inside window */
  860.       if ( viewvect[0] > xxmin && viewvect[0] < xxmax 
  861.         && viewvect[1] > yymin && viewvect[1] < yymax && ct > 100 )
  862.       {
  863.          col =          (( viewvect[0]-xxmin) / deltaX);
  864.          row = YYdots - (( viewvect[1]-yymin) / deltaY);
  865.          /* color is count of hits on this pixel */
  866.          color = getcolor(col,row)+1;
  867.          if( color < colors ) /* color sticks on last value */
  868.             (*plot)(col,row,color);
  869.       }
  870.    } /* end while */
  871.    return(0);
  872. }
  873. ifs3dlong()
  874. {
  875.    int bitshift;   /* want these local */
  876.    long fudge;
  877.    long delx,dely;
  878.       
  879.    double tmpx, tmpy, tmpz;
  880.    extern VECTOR view;
  881.    long iview[3];         /* perspective viewer's coordinates */      
  882.    long tmp;   
  883.    long maxifs[3];
  884.    long minifs[3];
  885.    extern int init3d[];
  886.    unsigned long maxct,ct;
  887.    register int col;
  888.    register int row;
  889.    register int color;
  890.  
  891.    long localifs[NUMIFS][IFS3DPARM];        /* local IFS values */
  892.    long newx,newy,newz,r,sum, xmin, xmax, ymin, ymax, tempr;
  893.     
  894.    int i,j,k;
  895.  
  896.    /* 3D viewing stuff */
  897.    MATRIX doublemat;           /* transformation matrix */
  898.    long longmat[4][4];         /* long version of matrix */
  899.    long ifsvect[3];               /* interated function orbit value */
  900.    long viewvect[3];        /* orbit transformed for viewing */
  901.  
  902.    /*
  903.    printf("xrot %d yrot %d zrot %d \nxshift %d yshift %d perspective %d\n",
  904.              XROT,YROT,ZROT,XSHIFT,YSHIFT,ZVIEWER);
  905.    */
  906.    bitshift = 16;
  907.    fudge = 1L << bitshift;
  908.  
  909.    delx = deltaX*fudge;
  910.    dely = deltaY*fudge;
  911.       
  912.    for(i=0;i<3;i++)
  913.    {
  914.       minifs[i] =  1L << 30;
  915.       maxifs[i] = -minifs[i];
  916.    }
  917.    
  918.    /* build transformation matrix */
  919.    identity (doublemat);
  920.  
  921.    /* apply rotations - uses the same rotation variables as line3d.c */
  922.    xrot ((double)XROT / 57.29577,doublemat);
  923.    yrot ((double)YROT / 57.29577,doublemat);
  924.    zrot ((double)ZROT / 57.29577,doublemat);
  925.  
  926.    /* apply scale */
  927.    scale((double)XSCALE/100.0,(double)YSCALE/100.0,(double)ROUGH/100.0,doublemat);
  928.  
  929.    if(!ZVIEWER)
  930.       trans((double)XSHIFT*(xxmax-xxmin)/100.0,(double)YSHIFT*(yymax-yymin)/100.0,0.0,doublemat);
  931.  
  932.    /* copy xform matrix to long for for fixed point math */
  933.    for (i = 0; i < 4; i++)
  934.    for (j = 0; j < 4; j++)
  935.       longmat[i][j] = doublemat[i][j] * fudge;
  936.  
  937.    if(diskvideo)        /* this would KILL a disk drive! */
  938.    {
  939.       setvideomode(3,0,0,0);
  940.       buzzer(2);
  941.       helpmessage(plasmamessage);
  942.       return(-1);
  943.    }
  944.  
  945.    for (i = 0; i < NUMIFS; i++)    /* fill in the local IFS array */
  946.    for (j = 0; j < IFS3DPARM; j++)
  947.          localifs[i][j] = initifs3d[i][j] * fudge;
  948.  
  949.    xmin  = xxmin*fudge;   
  950.    xmax  = xxmax*fudge;
  951.    ymax  = yymax*fudge;
  952.    ymin  = yymin*fudge;
  953.  
  954.    tempr = fudge / 32767;        /* find the proper rand() fudge */
  955.  
  956.    ifsvect[0] = 0;
  957.    ifsvect[1] = 0;
  958.    ifsvect[2] = 0;
  959.    
  960.    /* make maxct a function of screen size               */
  961.    /* 1k times maxit at EGA resolution seems about right */
  962.    maxct = (float)maxit*(1024.0*xdots*ydots)/(640.0*350.0); 
  963.    ct = 0L;
  964.    while(ct++ < maxct) /* loop until keypress or maxit */ 
  965.    {
  966.       if (ct & 127)           /* reduce the number of keypress checks */
  967.          if( check_key() )    /* keypress bails out */
  968.             return(-1);
  969.       r = rand();        /* generate fudged random number between 0 and 1 */
  970.       r *= tempr;
  971.  
  972.       /* pick which iterated function to execute, weighted by probability */
  973.       sum = localifs[0][12];
  974.       k = 0;
  975.       while ( sum < r)
  976.       {
  977.          k++;
  978.          sum += localifs[k][12];
  979.          if (localifs[k+1][12] == 0) break;    /* for safety */
  980.       }
  981.  
  982.       /* calculate image of last point under selected iterated function */
  983.       newx = multiply(localifs[k][0], ifsvect[0], bitshift) +
  984.              multiply(localifs[k][1], ifsvect[1], bitshift) +
  985.              multiply(localifs[k][2], ifsvect[2], bitshift) + localifs[k][9];
  986.       newy = multiply(localifs[k][3], ifsvect[0], bitshift) +
  987.              multiply(localifs[k][4], ifsvect[1], bitshift) +
  988.              multiply(localifs[k][5], ifsvect[2], bitshift) + localifs[k][10];
  989.       newz = multiply(localifs[k][6], ifsvect[0], bitshift) +
  990.              multiply(localifs[k][7], ifsvect[1], bitshift) +
  991.              multiply(localifs[k][8], ifsvect[2], bitshift) + localifs[k][11];
  992.             
  993.       ifsvect[0] = newx;
  994.       ifsvect[1] = newy;
  995.       ifsvect[2] = newz;
  996.  
  997.  
  998.       /* 3D VIEWING TRANSFORM */
  999.       longvmult(ifsvect,longmat,viewvect, bitshift);
  1000.  
  1001.       /* find minz and maxz */
  1002.       if( ZVIEWER && ct <= 100) /* waste this many points to find minz and maxz */
  1003.       {
  1004.          for(i=0;i<3;i++)
  1005.             if ((tmp = viewvect[i]) < minifs[i])
  1006.                minifs[i] = tmp;
  1007.             else if (tmp > maxifs[i])
  1008.                maxifs[i] = tmp;
  1009.          if(ct == 100)
  1010.          {
  1011.             /* set of view vector */
  1012.             for(i=0;i<2;i++)
  1013.                iview[i] = (maxifs[i]+minifs[i])/2;
  1014.       
  1015.             /* z value of user's eye - should be more negative than extreme 
  1016.                negative part of image */
  1017.             iview[2] = (long)((minifs[2]-maxifs[2])*(double)ZVIEWER/100.0);
  1018.  
  1019.             /*
  1020.             printf("min %ld max %ld iview[2] %d\n",minifs[2],maxifs[2],iview[2]);
  1021.             getch();
  1022.             */ 
  1023.  
  1024.             /* translate back exactly amount so maximum values are non-positive */
  1025.             tmpx = ((double)XSHIFT*(xmax-xmin))/(100.0*fudge);
  1026.             tmpy = ((double)YSHIFT*(ymax-ymin))/(100.0*fudge);
  1027.             tmpz = -((double)maxifs[2]);
  1028.             tmpz *= 1.1;
  1029.             tmpz /= fudge; 
  1030.             trans(tmpx,tmpy,tmpz,doublemat); 
  1031.             
  1032.             for(i=0;i<3;i++)
  1033.             {
  1034.                view[i] = iview[i];
  1035.                view[i] /= fudge;
  1036.             } 
  1037.    
  1038.             /* copy xform matrix to long for for fixed point math */
  1039.             for (i = 0; i < 4; i++)
  1040.             for (j = 0; j < 4; j++)
  1041.                longmat[i][j] = doublemat[i][j] * fudge;
  1042.          }
  1043.       }      
  1044.       
  1045.       /* apply perspective if requested */
  1046.       if(ZVIEWER && ct > 100)
  1047.       {
  1048.          if(debugflag==22 || ZVIEWER < 100) /* use float for small persp */
  1049.          {
  1050.             /* use float perspective calc */
  1051.             VECTOR tmpv;
  1052.             for(i=0;i<3;i++)
  1053.             {
  1054.                tmpv[i] = viewvect[i];
  1055.                tmpv[i] /= fudge;
  1056.             } 
  1057.             perspective(tmpv);
  1058.             for(i=0;i<3;i++)
  1059.                viewvect[i] = tmpv[i]*fudge;
  1060.          }
  1061.          else
  1062.             longpersp(viewvect,iview,bitshift);
  1063.       }
  1064.             
  1065.       /* plot if inside window */
  1066.       if ( viewvect[0] > xmin && viewvect[0] < xmax 
  1067.         && viewvect[1] > ymin && viewvect[1] < ymax && ct > 100 )
  1068.       {
  1069.          col =          (( viewvect[0]-xmin) / delx);
  1070.          row = YYdots - (( viewvect[1]-ymin) / dely);
  1071.  
  1072.          /* color is count of hits on this pixel */
  1073.          color = getcolor(col,row)+1;
  1074.          if( color < colors ) /* color sticks on last value */
  1075.             (*plot)(col,row,color);
  1076.       }
  1077.    }
  1078.    return(0);
  1079. }
  1080.  
  1081. /* START Phil Wilson's Code */
  1082.  
  1083. static    int        dem_start( void );
  1084. static    void        dem_end( void );
  1085. static    int        dem_pt( void );
  1086. static    double    MSetDist( void );
  1087. static    double    JSetDist( void );
  1088.  
  1089. #ifndef __TURBOC__
  1090.     void far * cdecl _fmalloc(size_t); /* from MSC's <malloc.h> */
  1091. #endif
  1092.  
  1093. #define BIG 100000.0
  1094. #define XOVERFLOW 100000000000000.0
  1095.  
  1096. static    int        mono, outside;
  1097. static    double    delta, pixelwidth;
  1098. static    double    (*DistEstimate)();
  1099. static    struct    complex far *orbit;
  1100. static    struct    complex deriv, squared;
  1101.  
  1102. DistanceEstimatorMethod()
  1103. {
  1104.    if ( ! dem_start()) return(-1);
  1105.    for (passes=0; passes <= numpasses ; passes++)
  1106.    for (row = passes; row < YYdots; row=row+1+numpasses)
  1107.    {
  1108.       register int col;
  1109.       init.y = dy0[row];
  1110.       for (col = 0; col < XXdots; col++)       /* look at each point on screen */
  1111.       {
  1112.          register color;
  1113.          init.x = dx0[col];
  1114.          if(check_key()) 
  1115.          {
  1116.             dem_end();
  1117.             return(-1);
  1118.          }
  1119.          color = dem_pt();
  1120.          (*plot)(col,row,color);
  1121.          if(numpasses && (passes == 0))
  1122.              (*plot)(col,row+1,color);
  1123.       }
  1124.    }
  1125.    dem_end();
  1126.    return(0);
  1127. }
  1128.  
  1129. static int dem_start( void )
  1130. {
  1131.     unsigned long array_size = (maxit+1) * sizeof(*orbit);
  1132.  
  1133. #ifdef __TURBOC__
  1134.     if (( orbit = (struct complex far *) farmemalloc( array_size )) == NULL)
  1135.     {
  1136.         printf("\n\n\n\n\n\n\n\n\n\n\n\n");
  1137.         printf("Not enough Memory for array with %d orbits\n", maxit);
  1138.         printf(" (%ul Bytes Required)", array_size);
  1139.         buzzer(2);
  1140.         dem_end();
  1141.         return(0);
  1142.     }
  1143. #else    /* Assume MSC, whose _fmalloc seems able to find more memory than */
  1144.         /* farmemalloc().  Important for those maxiter=4096 all-nighters  */
  1145.  
  1146.     if (( orbit = (struct complex far *) _fmalloc( array_size )) == NULL)
  1147.     {
  1148.         printf("\n\n\n\n\n\n\n\n\n\n\n\n");
  1149.         printf("Not enough Memory for array with %d orbits\n", maxit);
  1150.         printf(" (%ul Bytes Required)", array_size);
  1151.         buzzer(2);
  1152.         dem_end();
  1153.         return(0);
  1154.     }
  1155. #endif    
  1156.  
  1157.     delta = (deltaX + deltaY) / 4;    /* half a pixel width */
  1158.     pixelwidth = 2*delta;
  1159.  
  1160.     if ( colors == 2 ) mono = 1;
  1161.     if ( mono ) outside = !inside;
  1162.  
  1163.     if ( strncmp( fractalspecific[fractype].name, "demm", 4 ) == 0)
  1164.         DistEstimate = MSetDist;
  1165.     else if ( strncmp( fractalspecific[fractype].name, "demj", 4 ) == 0)
  1166.         DistEstimate = JSetDist;
  1167.     else return(0);
  1168.  
  1169.     return(1); /* everything's OK, proceed */
  1170. }
  1171.  
  1172. static void dem_end( void )
  1173. {
  1174. #ifdef __TURBOC__
  1175.     farmemfree((char far *) orbit);
  1176. #else
  1177.     _ffree((void far *) orbit);
  1178. #endif
  1179. }
  1180.  
  1181. static int dem_pt( void )
  1182. {
  1183.     int        color;
  1184.     double    tempcolor, distance;
  1185.  
  1186.     if ((distance = DistEstimate()) < delta)    color = inside;
  1187.     else {
  1188.         if (mono) color = outside;
  1189.         else {
  1190.             tempcolor = 1 + (( distance/pixelwidth ));
  1191.             color = (int)(tempcolor / sqrt( tempcolor )) % colors;
  1192.         }
  1193.     }
  1194.     return(color);
  1195. }
  1196.  
  1197. /*     Distance estimator for points near Mandelbrot set                            */
  1198. /*     Algorithms from Peitgen & Saupe, Science of Fractal Images, p.198     */
  1199.  
  1200. static double MSetDist()
  1201. {
  1202.     int        iter = 0, i, flag;
  1203.     double    x, y, temp, dist;
  1204.  
  1205.     x = y = squared.x = squared.y = dist = orbit[0].x = orbit[0].y = 0.0;
  1206.  
  1207.     while ((iter <    maxit) && (squared.x + squared.y < BIG)) 
  1208.     {
  1209.         temp = squared.x - squared.y + init.x;
  1210.         iter++;
  1211.         orbit[iter].y = y = 2 * x * y + init.y;
  1212.         orbit[iter].x = x = temp;
  1213.         squared.x = x * x;
  1214.         squared.y = y * y;
  1215.     }
  1216.     if ( squared.x + squared.y > BIG ) 
  1217.     {
  1218.         deriv.x = deriv.y = 0.0;
  1219.         i = flag = 0;
  1220.         while ((i < iter) && (!flag)) 
  1221.         {
  1222.             temp = 2 * (orbit[i].x * deriv.x - orbit[i].y * deriv.y) + 1;
  1223.             deriv.y = 2 * (orbit[i].y * deriv.x + orbit[i].x * deriv.y);
  1224.             deriv.x = temp;
  1225.             if ( max( fabs( deriv.x ), fabs( deriv.y )) > XOVERFLOW ) flag++;
  1226.             i++;
  1227.         }
  1228.         if ( !flag )
  1229.             dist = log(squared.x + squared.y) * sqrt(squared.x + squared.y)
  1230.                                      /    sqrt(deriv.x * deriv.x     + deriv.y * deriv.y);
  1231.     }
  1232.     return(dist);
  1233. }
  1234.  
  1235. static double JSetDist()
  1236. {
  1237.     int        iter = 0, i;
  1238.     int        flag;
  1239.     double    x, y, temp,    dist;
  1240.  
  1241.     x = init.x;
  1242.     y = init.y;
  1243.     squared.x = x * x;
  1244.     squared.y = y * y;
  1245.     dist = orbit[0].x = orbit[0].y = 0.0;
  1246.     
  1247.     while ((iter <    maxit)  && (squared.x + squared.y < BIG) ) {
  1248.         temp = squared.x - squared.y + parm.x;
  1249.         y = 2 * x * y + parm.y;
  1250.         x = temp;
  1251.         squared.x = x * x;
  1252.         squared.y = y * y;
  1253.         iter++;
  1254.         orbit[ iter ].x = x;
  1255.         orbit[ iter ].y = y;
  1256.     }
  1257.     if ( squared.x + squared.y > BIG ) {
  1258.         deriv.x = deriv.y = 0.0;
  1259.         i = 0;
  1260.         flag = 0;
  1261.         while ((i < iter) && (!flag)) {
  1262.             temp = 2 * (orbit[i].x * deriv.x - orbit[i].y * deriv.y) + 1;
  1263.             deriv.y = 2 * (orbit[i].y * deriv.x + orbit[i].x * deriv.y);
  1264.             deriv.x = temp;
  1265.             if ( max( fabs( deriv.x ), fabs( deriv.y )) > XOVERFLOW ) flag++;
  1266.             i++;
  1267.         }
  1268.         if ( !flag )
  1269.             dist = log(squared.x + squared.y) * sqrt(squared.x + squared.y) / sqrt(deriv.x * deriv.x
  1270.                                                                     + deriv.y * deriv.y);
  1271.     }
  1272.     return(dist);
  1273. }
  1274.  
  1275. #define DEFAULTFILTER 1000    /* "Beauty of Fractals" recommends using 5000 
  1276.                                (p.25), but that seems unnecessary. Can 
  1277.                                override this value with a nonzero param1 */
  1278.  
  1279. #define SEED 0.66        /* starting value for population */
  1280.  
  1281. static void verhulst( double rate );
  1282.  
  1283. int far *verhulst_array;    
  1284. unsigned int filter_cycles;
  1285.  
  1286. int Bifurcation( void )
  1287. {
  1288.    unsigned long array_size;
  1289.    int row, column;
  1290.  
  1291.    array_size =  (YYdots + 1) * sizeof(int);
  1292.  
  1293.     if ( (verhulst_array = (int far *) farmemalloc(array_size)) == NULL)
  1294.     {
  1295.         printf("Better bail out, not enough heap for array of %d int\'s\n",
  1296.                     YYdots );
  1297.         buzzer(2);
  1298.     }
  1299.  
  1300.     for (row = 0; row <= YYdots; row++)
  1301.         verhulst_array[row] = 0;
  1302.  
  1303.     mono = 0;
  1304.     if ( colors == 2 ) mono = 1;
  1305.     if ( mono )    {
  1306.         if ( inside ) {
  1307.             outside = 0;
  1308.             inside = 1;
  1309.         }
  1310.         else outside = 1;
  1311.     }
  1312.  
  1313.     if ((filter_cycles = parm.x) == 0) filter_cycles = DEFAULTFILTER;
  1314.    
  1315.     init.y = dy0[YYdots - 1]; /* Y value of bottom visible pixel */
  1316.     
  1317.     for (column = 0; column < XXdots; column++)
  1318.    {
  1319.         verhulst( dx0[column] );     /* calculate array once per column */
  1320.  
  1321.         for (row = YYdots; row > 0; row--)
  1322.       {
  1323.          int color;
  1324.          color = verhulst_array[ row ];
  1325.             if ( color && mono ) color = inside;
  1326.             else if ( (!color) && mono ) color = outside;
  1327.             verhulst_array[ row ] = 0;
  1328.             (*plot)(column,row,color);
  1329.       }
  1330.       if(check_key()) 
  1331.       {
  1332.             farmemfree( (char far *)verhulst_array );
  1333.             return(-1);
  1334.       }
  1335.    }
  1336.     farmemfree( (char far *)verhulst_array );
  1337.    return(0);
  1338. }
  1339.  
  1340. static void verhulst( double rate )  /* P. F. Verhulst (1845) */
  1341. {
  1342.     double    OldPopulation, NewPopulation;
  1343.     unsigned    int pixel_row, counter;
  1344.  
  1345.     OldPopulation = SEED;
  1346.  
  1347.     for ( counter = 0; counter < filter_cycles; counter++ )
  1348.     {
  1349.         NewPopulation = ((1 + rate) * OldPopulation ) -
  1350.             ( rate * OldPopulation * OldPopulation );
  1351.         OldPopulation = NewPopulation;
  1352.     }
  1353.     for ( counter = 0; counter < maxit; counter++ )
  1354.     {
  1355.         NewPopulation = ((1 + rate) * OldPopulation ) -
  1356.             ( rate * OldPopulation * OldPopulation );
  1357.  
  1358.         OldPopulation = NewPopulation;
  1359.  
  1360.         /* assign population value to Y coordinate in pixels */
  1361.         pixel_row = YYdots - (int)((NewPopulation - init.y) / deltaY);
  1362.  
  1363.         /* if it's visible on the screen, save it in the column array */
  1364.         if ((pixel_row >= 0) && (pixel_row <= YYdots))    
  1365.                             verhulst_array[ pixel_row ] ++;
  1366.     }
  1367. }
  1368.  
  1369. /* END Phil Wilson's Code */
  1370.  
  1371.  
  1372. /* -------------------------------------------------------------------- */
  1373. /*        Bailout Routines Macros                                                  */
  1374. /* -------------------------------------------------------------------- */
  1375.  
  1376. #define FLOATBAILOUT()   \
  1377.    tempsqrx = sqr(new.x);\
  1378.    tempsqry = sqr(new.y);\
  1379.    if((tempsqrx + tempsqry) >= rqlim) return(1);\
  1380.    old = new;
  1381.  
  1382. #define LONGBAILOUT()   \
  1383.    ltempsqrx = lsqr(lnew.x);\
  1384.    ltempsqry = lsqr(lnew.y);\
  1385.    lmagnitud = ltempsqrx + ltempsqry;\
  1386.    if (lmagnitud >= llimit || lmagnitud < 0 || labs(lnew.x) > llimit2\
  1387.          || labs(lnew.y) > llimit2) \
  1388.       return(1);\
  1389.    lold = lnew;
  1390.  
  1391. #define FLOATTRIGBAILOUT()  \
  1392.    if (fabs(old.y) >= rqlim2) return(1);
  1393.  
  1394. #define LONGTRIGBAILOUT()  \
  1395.    if(labs(lold.y) >= llimit2) return(1);
  1396.  
  1397. #define FLOATHTRIGBAILOUT()  \
  1398.    if (fabs(old.x) >= rqlim2) return(1);
  1399.  
  1400. #define LONGHTRIGBAILOUT()  \
  1401.    if(labs(lold.x) >= llimit2) return(1);
  1402.  
  1403. #define TRIG16CHECK(X)  \
  1404.       if(labs((X)) > l16triglim) return(1);
  1405.  
  1406. #if 0
  1407. /* this define uses usual trig instead of fast trig */
  1408. #define FPUsincos(px,psinx,pcosx) \
  1409.    *(psinx) = sin(*(px));\
  1410.    *(pcosx) = cos(*(px));
  1411.    
  1412. #define FPUsinhcosh(px,psinhx,pcoshx) \
  1413.    *(psinhx) = sinh(*(px));\
  1414.    *(pcoshx) = cosh(*(px));
  1415. #endif
  1416.  
  1417. #define LTRIGARG(X)    \
  1418.    if(labs((X)) > l16triglim)\
  1419.    {\
  1420.       double tmp;\
  1421.       tmp = (X);\
  1422.       tmp /= fudge;\
  1423.       tmp = fmod(tmp,twopi);\
  1424.       tmp *= fudge;\
  1425.       (X) = tmp;\
  1426.    }\
  1427.     
  1428. /* -------------------------------------------------------------------- */
  1429. /*        Fractal (once per iteration) routines            */
  1430. /* -------------------------------------------------------------------- */
  1431. static double xt, yt, t2;
  1432.  
  1433. /* Raise complex number (base) to the (exp) power, storing the result
  1434. ** in complex (result).
  1435. */
  1436. cpower(struct complex *base, int exp, struct complex *result)
  1437. {
  1438.     xt = base->x;   yt = base->y;
  1439.  
  1440.     if (exp & 1) 
  1441.     {  
  1442.        result->x = xt;     
  1443.        result->y = yt;     
  1444.     }
  1445.     else         
  1446.     {  
  1447.        result->x = 1.0;    
  1448.        result->y = 0.0;    
  1449.     }
  1450.  
  1451.     exp >>= 1;
  1452.     while (exp) 
  1453.     {
  1454.         t2 = xt * xt - yt * yt;
  1455.         yt = 2 * xt * yt;
  1456.         xt = t2;
  1457.  
  1458.         if (exp & 1) 
  1459.         {
  1460.             t2 = xt * result->x - yt * result->y;
  1461.             result->y = result->y * xt + yt * result->x;
  1462.             result->x = t2;
  1463.         }
  1464.         exp >>= 1;
  1465.     }
  1466. }
  1467. /* long version */
  1468. static long lxt, lyt, lt1, lt2;
  1469. lcpower(struct lcomplex *base, int exp, struct lcomplex *result, int bitshift)
  1470. {
  1471.     static long maxarg;
  1472.     maxarg = 64L<<bitshift;
  1473.  
  1474.     overflow = 0;    
  1475.     lxt = base->x;   lyt = base->y;
  1476.  
  1477.     if (exp & 1) 
  1478.     {  
  1479.        result->x = lxt;     
  1480.        result->y = lyt;     
  1481.     }
  1482.     else         
  1483.     {  
  1484.        result->x = 1L<<bitshift;    
  1485.        result->y = 0L;    
  1486.     }
  1487.  
  1488.     exp >>= 1;
  1489.     while (exp) 
  1490.     {
  1491.         /*
  1492.         if(labs(lxt) >= maxarg || labs(lyt) >= maxarg)
  1493.            return(-1);
  1494.         */    
  1495.         lt2 = multiply(lxt, lxt, bitshift) - multiply(lyt,lyt,bitshift);
  1496.         lyt = multiply(lxt,lyt,bitshift)<<1;
  1497.         if(overflow)
  1498.            return(overflow);
  1499.         lxt = lt2;
  1500.  
  1501.         if (exp & 1) 
  1502.         {
  1503.             lt2 = multiply(lxt,result->x, bitshift) - multiply(lyt,result->y,bitshift);
  1504.             result->y = multiply(result->y,lxt,bitshift) + multiply(lyt,result->x,bitshift);
  1505.             result->x = lt2;
  1506.         }
  1507.         exp >>= 1;
  1508.     }
  1509.     if(result->x == 0 && result->y == 0) 
  1510.        overflow = 1;
  1511.     return(overflow);
  1512. }
  1513.  
  1514. z_to_the_z(struct complex *z, struct complex *out)
  1515. {
  1516.     static struct complex tmp1,tmp2;
  1517.     /* raises complex z to the z power */
  1518.     extern int errno;
  1519.     errno = 0;
  1520.  
  1521.     if(fabs(z->x) < DBL_EPSILON) return(-1);
  1522.  
  1523.     /* log(x + iy) = 1/2(log(x*x + y*y) + i(arc_tan(y/x)) */
  1524.     tmp1.x = .5*log(sqr(z->x)+sqr(z->y));
  1525.     
  1526.     /* the labs in next line added to prevent discontinuity in image */
  1527.     tmp1.y = atan(labs(z->y/z->x));
  1528.  
  1529.     /* log(z)*z */
  1530.     tmp2.x = tmp1.x * z->x - tmp1.y * z->y;
  1531.     tmp2.y = tmp1.x * z->y + tmp1.y * z->x;
  1532.  
  1533.     /* z*z = e**(log(z)*z) */
  1534.     /* e**(x + iy) =  e**x * (cos(y) + isin(y)) */
  1535.  
  1536.     tmpexp = exp(tmp2.x);
  1537.  
  1538.     FPUsincos(&tmp2.y,&siny,&cosy);
  1539.     out->x = tmpexp*cosy;
  1540.     out->y = tmpexp*siny;
  1541.     return(errno);
  1542. }
  1543.  
  1544. /* Distance of complex z from unit circle */
  1545. #define DIST1(z) (((z).x-1.0)*((z).x-1.0)+((z).y)*((z).y))
  1546. #define LDIST1(z) (lsqr((((z).x)-fudge)) + lsqr(((z).y)))
  1547.  
  1548. #ifdef NEWTON
  1549.  
  1550. int NewtonFractal()
  1551. {
  1552.     static char start=1;
  1553.     if(start)
  1554.     {
  1555.        printf("c version");
  1556.        start = 0;
  1557.     }
  1558.     cpower(&old, degree-1, &tmp);
  1559.     complex_mult(tmp, old, &new);
  1560.  
  1561.     if (DIST1(new) < threshold) 
  1562.     {
  1563.        if(fractype==NEWTBASIN)
  1564.        {
  1565.           int tmpcolor;
  1566.           int i;
  1567.           tmpcolor = -1;
  1568.           /* this code determines which degree-th root of root the        
  1569.              Newton formula converges to. The roots of a 1 are 
  1570.              distributed on a circle of radius 1 about the origin. */
  1571.           for(i=0;i<degree;i++)
  1572.              /* color in alternating shades with iteration according to
  1573.                 which root of 1 it converged to */
  1574.               if( distance(roots[i],old) < threshold)
  1575.               { 
  1576.                /*    tmpcolor = 1+(i&7)+((color&1)<<3); */
  1577.                    tmpcolor = 1+i;
  1578.                    break;
  1579.               }     
  1580.            if(tmpcolor == -1)
  1581.               color = maxcolor;
  1582.            else
  1583.               color = tmpcolor;
  1584.        }
  1585.        return(1);
  1586.     }   
  1587.     new.x = d1overd * new.x + roverd;
  1588.     new.y *= d1overd;
  1589.  
  1590.     /* Watch for divide underflow */
  1591.     if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN) 
  1592.       return(1);
  1593.     else 
  1594.     {
  1595.         t2 = 1.0 / t2;
  1596.         old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
  1597.         old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
  1598.     }
  1599.     return(0);
  1600. }
  1601.  
  1602. #endif
  1603.  
  1604. complex_mult(arg1,arg2,pz)
  1605. struct complex arg1,arg2,*pz;
  1606. {
  1607.    pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
  1608.    pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
  1609.    return(0);
  1610. }
  1611.  
  1612. complex_div(numerator,denominator,pout)
  1613. struct complex numerator,denominator,*pout;
  1614. {  
  1615.    double mod;
  1616.    if((mod = modulus(denominator)) < FLT_MIN)
  1617.       return(1);
  1618.    conjugate(&denominator);
  1619.    complex_mult(numerator,denominator,pout);
  1620.    pout->x = pout->x/mod;
  1621.    pout->y = pout->y/mod;
  1622.    return(0);
  1623. }
  1624.  
  1625.  
  1626. lcomplex_mult(arg1,arg2,pz,bitshift)
  1627. struct lcomplex arg1,arg2,*pz;
  1628. int bitshift;
  1629. {
  1630.    overflow = 0;
  1631.    pz->x = multiply(arg1.x,arg2.x,bitshift) - multiply(arg1.y,arg2.y,bitshift);
  1632.    pz->y = multiply(arg1.x,arg2.y,bitshift) + multiply(arg1.y,arg2.x,bitshift);
  1633.    return(overflow);
  1634. }
  1635.  
  1636. #define MPCmod(x) (pMPadd(pMPmul((x).r, (x).r), pMPmul((x).i, (x).i)))
  1637. struct MPC mpcold, mpcnew, mpctmp, mpctmp1;
  1638. struct MP mproverd, mpd1overd, mpthreshold,sqrmpthreshold;
  1639. struct MP mpt2;
  1640. struct MP mpone;
  1641. extern struct MPC MPCone;
  1642. extern int MPOverflow;
  1643.  
  1644. int MPCNewtonFractal()
  1645. {
  1646.     MPOverflow = 0;
  1647.     mpctmp   = MPCpow(mpcold,degree-1); 
  1648.  
  1649.     mpcnew.r = pMPsub(pMPmul(mpctmp.r,mpcold.r),pMPmul(mpctmp.i,mpcold.i));
  1650.     mpcnew.i = pMPadd(pMPmul(mpctmp.r,mpcold.i),pMPmul(mpctmp.i,mpcold.r));
  1651.  
  1652.     mpctmp1.r = pMPsub(mpcnew.r, MPCone.r);
  1653.     mpctmp1.i = pMPsub(mpcnew.i, MPCone.i);
  1654.  
  1655.     if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
  1656.     {
  1657.       if(fractype==MPNEWTBASIN)
  1658.       {
  1659.          int tmpcolor;
  1660.          int i;
  1661.          tmpcolor = -1;
  1662.          for(i=0;i<degree;i++)
  1663.              if( pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
  1664.              { 
  1665.                   tmpcolor = 1+i;
  1666.                   break;
  1667.              }     
  1668.           if(tmpcolor == -1)
  1669.              color = maxcolor;
  1670.           else
  1671.              color = tmpcolor;
  1672.       }
  1673.        return(1);
  1674.     }   
  1675.  
  1676.     mpcnew.r = pMPadd(pMPmul(mpd1overd,mpcnew.r),mproverd);
  1677.     mpcnew.i = pMPmul(mpcnew.i,mpd1overd);
  1678.  
  1679.     mpt2 = MPCmod(mpctmp);
  1680.     mpt2 = pMPdiv(mpone,mpt2);
  1681.     
  1682.     mpcold.r = pMPmul(mpt2,(pMPadd(pMPmul(mpcnew.r,mpctmp.r),pMPmul(mpcnew.i,mpctmp.i))));
  1683.     mpcold.i = pMPmul(mpt2,(pMPsub(pMPmul(mpcnew.i,mpctmp.r),pMPmul(mpcnew.r,mpctmp.i))));
  1684.  
  1685.     new.x = *pMP2d(mpcold.r);
  1686.     new.y = *pMP2d(mpcold.i);
  1687.     return(MPOverflow);
  1688. }
  1689.  
  1690. JuliaFractal() 
  1691. {
  1692.    /* used for C prototype of fast integer math routines for classic
  1693.       Mandelbrot and Julia */ 
  1694.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  1695.    lnew.y = (multiply(lold.x, lold.y, bitshift) << 1) + longparm->y;
  1696.    LONGBAILOUT();
  1697.    return(0);
  1698. }
  1699.  
  1700. Barnsley1Fractal() 
  1701. {
  1702.    /* Barnsley's Mandelbrot type M1 from "Fractals
  1703.    Everywhere" by Michael Barnsley, p. 322 */ 
  1704.  
  1705.    /* calculate intermediate products */
  1706.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  1707.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  1708.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  1709.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  1710.    /* orbit calculation */                  
  1711.    if(lold.x >= 0) 
  1712.    {
  1713.       lnew.x = (oldxinitx - longparm->x - oldyinity);               
  1714.       lnew.y = (oldyinitx - longparm->y + oldxinity);
  1715.    }
  1716.    else 
  1717.    {
  1718.       lnew.x = (oldxinitx + longparm->x - oldyinity);               
  1719.       lnew.y = (oldyinitx + longparm->y + oldxinity);               
  1720.    }
  1721.    LONGBAILOUT();
  1722.    return(0);
  1723. }
  1724.  
  1725. Barnsley2Fractal() 
  1726. {
  1727.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1728.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  1729.     
  1730.    /* calculate intermediate products */
  1731.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  1732.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  1733.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  1734.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  1735.  
  1736.    /* orbit calculation */                  
  1737.    if(oldxinity + oldyinitx >= 0) 
  1738.    {
  1739.       lnew.x = oldxinitx - longparm->x - oldyinity;               
  1740.       lnew.y = oldyinitx - longparm->y + oldxinity;               
  1741.    }
  1742.    else 
  1743.    {
  1744.       lnew.x = oldxinitx + longparm->x - oldyinity;               
  1745.       lnew.y = oldyinitx + longparm->y + oldxinity;               
  1746.    }
  1747.    LONGBAILOUT();
  1748.    return(0);
  1749. }
  1750.  
  1751. JuliafpFractal() 
  1752. {
  1753.    /* floating point version of classical Mandelbrot/Julia */
  1754.    new.x = tempsqrx - tempsqry + floatparm->x;
  1755.    new.y = 2.0 * old.x * old.y + floatparm->y;
  1756.    FLOATBAILOUT();
  1757.    return(0);
  1758. }
  1759.  
  1760. LambdaFractal() 
  1761. {
  1762.    /* variation of classical Mandelbrot/Julia */
  1763.    
  1764.    /* in complex math) temp = Z * (1-Z) */
  1765.    ltempsqrx = lold.x  - multiply(lold.x, lold.x, bitshift)
  1766.                  + multiply(lold.y, lold.y, bitshift);
  1767.    ltempsqry = lold.y
  1768.                  - (multiply(lold.y, lold.x, bitshift)  << 1);
  1769.    /* (in complex math) Z = Lambda * Z */
  1770.    lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
  1771.         - multiply(longparm->y, ltempsqry, bitshift);
  1772.    lnew.y = multiply(longparm->x, ltempsqry, bitshift)
  1773.         + multiply(longparm->y, ltempsqrx, bitshift);
  1774.    LONGBAILOUT();
  1775.    return(0);
  1776. }
  1777.  
  1778. SierpinskiFractal() 
  1779. {
  1780.    /* following code translated from basic - see "Fractals
  1781.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */ 
  1782.    lnew.x = (lold.x << 1);        /* new.x = 2 * old.x  */
  1783.    lnew.y = (lold.y << 1);        /* new.y = 2 * old.y  */
  1784.    if( lold.y > ltemp.y)        /* if old.y > .5 */
  1785.       lnew.y = lnew.y - ltemp.x;    /* new.y = 2 * old.y - 1 */
  1786.    else if(lold.x > ltemp.y)    /* if old.x > .5 */
  1787.       lnew.x = lnew.x - ltemp.x;    /* new.x = 2 * old.x - 1 */
  1788.    /* end barnsley code */
  1789.    LONGBAILOUT();
  1790.    return(0);
  1791. }
  1792.  
  1793. static long lcosx, lcoshx, lsinx, lsinhx;
  1794. static long lcosy, lcoshy, lsiny, lsinhy;
  1795.  
  1796. LongLambdasineFractal() 
  1797. {
  1798.    LONGTRIGBAILOUT();
  1799.  
  1800.    SinCos086(lold.x, &lsinx, &lcosx);
  1801.    SinhCosh086(lold.y, &lsinhy, &lcoshy);
  1802.  
  1803.    ltemp.x = multiply(lsinx,        lcoshy,   bitshift);
  1804.    ltemp.y = multiply(lcosx,        lsinhy,   bitshift);
  1805.  
  1806.    lnew.x  = multiply(longparm->x, ltemp.x, bitshift)
  1807.            - multiply(longparm->y, ltemp.y, bitshift);
  1808.    lnew.y  = multiply(longparm->x, ltemp.y, bitshift)
  1809.            + multiply(longparm->y, ltemp.x, bitshift);
  1810.    lold = lnew;
  1811.    return(0);
  1812. }
  1813.  
  1814. LongLambdacosFractal() 
  1815. {
  1816.    LONGTRIGBAILOUT();
  1817.  
  1818.    SinCos086(lold.x, &lsinx, &lcosx);
  1819.    SinhCosh086(lold.y, &lsinhy, &lcoshy);
  1820.  
  1821.    ltemp.x = multiply(lcosx,        lcoshy,   bitshift);
  1822.    ltemp.y = multiply(lsinx,        lsinhy,   bitshift);
  1823.  
  1824.    lnew.x  = multiply(longparm->x, ltemp.x, bitshift)
  1825.            - multiply(longparm->y, ltemp.y, bitshift);
  1826.    lnew.y  = multiply(longparm->x, ltemp.y, bitshift)
  1827.            + multiply(longparm->y, ltemp.x, bitshift);
  1828.    lold = lnew;
  1829.    return(0);
  1830. }
  1831.  
  1832. LongLambdasinhFractal() 
  1833. {
  1834.    LONGHTRIGBAILOUT();
  1835.  
  1836.    SinCos086  (lold.y, &lsiny,  &lcosy);
  1837.    SinhCosh086(lold.x, &lsinhx, &lcoshx);
  1838.  
  1839.    ltemp.x = multiply(lcosy,        lsinhx,   bitshift);
  1840.    ltemp.y = multiply(lsiny,        lcoshx,   bitshift);
  1841.  
  1842.    lnew.x  = multiply(longparm->x, ltemp.x, bitshift)
  1843.            - multiply(longparm->y, ltemp.y, bitshift);
  1844.    lnew.y  = multiply(longparm->x, ltemp.y, bitshift)
  1845.            + multiply(longparm->y, ltemp.x, bitshift);
  1846.    lold = lnew;
  1847.    return(0);
  1848. }
  1849.  
  1850. LongLambdacoshFractal() 
  1851. {
  1852.    LONGHTRIGBAILOUT();
  1853.  
  1854.    SinCos086(lold.y, &lsiny, &lcosy);
  1855.    SinhCosh086(lold.x, &lsinhx, &lcoshx);
  1856.  
  1857.    ltemp.x = multiply(lcosy,        lcoshx,   bitshift);
  1858.    ltemp.y = multiply(lsiny,        lsinhx,   bitshift);
  1859.  
  1860.    lnew.x  = multiply(longparm->x, ltemp.x, bitshift)
  1861.            - multiply(longparm->y, ltemp.y, bitshift);
  1862.    lnew.y  = multiply(longparm->x, ltemp.y, bitshift)
  1863.            + multiply(longparm->y, ltemp.x, bitshift);
  1864.    lold = lnew;
  1865.    return(0);
  1866. }
  1867.  
  1868. LambdasineFractal() 
  1869. {
  1870.    /* found this in  "Science of Fractal Images" */ 
  1871.  
  1872.    FLOATTRIGBAILOUT();
  1873.  
  1874.    /* calculate sin(z) */
  1875.    FPUsincos  (&old.x,&sinx,&cosx);
  1876.    FPUsinhcosh(&old.y,&sinhy,&coshy);
  1877.    tmp.x = sinx*coshy;
  1878.    tmp.y = cosx*sinhy;
  1879.  
  1880.    /*multiply by lamda */
  1881.    new.x = floatparm->x * tmp.x - floatparm->y * tmp.y;
  1882.    new.y = floatparm->y * tmp.x + floatparm->x * tmp.y;
  1883.    old = new;
  1884.    return(0);
  1885. }
  1886.  
  1887. LambdacosineFractal() 
  1888. {
  1889.    /* found this in  "Science of Fractal Images" */ 
  1890.  
  1891.    FLOATTRIGBAILOUT();
  1892.  
  1893.    /* calculate cos(z) */
  1894.    FPUsincos  (&old.x,&sinx,&cosx);
  1895.    FPUsinhcosh(&old.y,&sinhy,&coshy);
  1896.    tmp.x = cosx*coshy;
  1897.    tmp.y = sinx*sinhy;
  1898.  
  1899.    /*multiply by lamda */
  1900.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  1901.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  1902.    old = new;
  1903.    return(0);
  1904. }
  1905.  
  1906. LambdasinhFractal() 
  1907. {
  1908.    FLOATHTRIGBAILOUT();
  1909.  
  1910.    /* calculate sinh(z) */
  1911.  
  1912.    FPUsincos  (&old.y,&siny,&cosy);
  1913.    FPUsinhcosh(&old.x,&sinhx,&coshx);
  1914.    tmp.x = sinhx*cosy;
  1915.    tmp.y = coshx*siny;
  1916.    
  1917.    /*multiply by lamda */
  1918.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  1919.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  1920.    old = new;
  1921.    return(0);
  1922. }
  1923.  
  1924. LambdacoshFractal() 
  1925. {
  1926.    FLOATHTRIGBAILOUT();
  1927.  
  1928.    /* calculate cosh(z) */
  1929.    FPUsincos  (&old.y,&siny,&cosy);
  1930.    FPUsinhcosh(&old.x,&sinhx,&coshx);
  1931.    tmp.x = coshx*cosy;
  1932.    tmp.y = sinhx*siny;
  1933.  
  1934.    /*multiply by lamda */
  1935.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  1936.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  1937.    old = new;
  1938.    return(0);
  1939. }
  1940.  
  1941. LambdaexponentFractal() 
  1942. {
  1943.    /* found this in  "Science of Fractal Images" */ 
  1944.  
  1945.    /* calculate exp(z) */
  1946.    if (fabs(old.y) >= 1.0e8) return(1); /* TLOSS  errors */
  1947.    if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
  1948.  
  1949.    FPUsincos  (&old.y,&siny,&cosy);
  1950.    
  1951.    if (old.x >= rqlim && cosy >= 0.0) return(1);
  1952.    tmpexp = exp(old.x);
  1953.    tmp.x = tmpexp*cosy;
  1954.    tmp.y = tmpexp*siny;
  1955.  
  1956.    /*multiply by lamda */
  1957.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  1958.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  1959.    old = new;
  1960.    return(0);
  1961. }
  1962.  
  1963. LongLambdaexponentFractal() 
  1964. {
  1965.    /* found this in  "Science of Fractal Images" */ 
  1966.  
  1967.    /* calculate exp(z) */
  1968.    if (labs(lold.y) >= 1000L<<bitshift) return(1); /* TLOSS  errors */
  1969.    if (labs(lold.x) >=  8L<<bitshift) return(1); /* DOMAIN errors */
  1970.  
  1971.    SinCos086  (lold.y, &lsiny,  &lcosy);
  1972.  
  1973.    if (lold.x >= llimit && lcosy >= 0L) return(1);
  1974.    ltmp1 = Exp086(lold.x);
  1975.  
  1976.    ltemp.x = multiply(ltmp1,       lcosy,   bitshift);
  1977.    ltemp.y = multiply(ltmp1,       lsiny,   bitshift);
  1978.  
  1979.    lnew.x  = multiply(longparm->x, ltemp.x, bitshift)
  1980.            - multiply(longparm->y, ltemp.y, bitshift);
  1981.    lnew.y  = multiply(longparm->x, ltemp.y, bitshift)
  1982.            + multiply(longparm->y, ltemp.x, bitshift);
  1983.    lold = lnew;
  1984.    return(0);
  1985. }
  1986.  
  1987. FloatSinexponentFractal() 
  1988. {
  1989.    /* another Scientific American biomorph type */
  1990.    /* z(n+1) = e**z(n) + sin(z(n)) + C */
  1991.  
  1992.    if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
  1993.    tmpexp = exp(old.x);
  1994.    
  1995.    FPUsincos  (&old.y,&siny,&cosy);
  1996.    FPUsincos  (&old.x,&sinx,&cosx);
  1997.    FPUsinhcosh(&old.y,&sinhy,&coshy);
  1998.  
  1999.    /*new = e**old      + sin(old)   + C  */     
  2000.    new.x = tmpexp*cosy + sinx*coshy + floatparm->x;
  2001.    new.y = tmpexp*siny + cosx*sinhy + floatparm->y;
  2002.  
  2003.    FLOATBAILOUT();
  2004.    return(0);
  2005. }
  2006.  
  2007.  
  2008. LongSinexponentFractal() 
  2009. {
  2010.    /* calculate exp(z) */
  2011.    
  2012.    /* domain check for fast transcendental functions */
  2013.    TRIG16CHECK(lold.x);
  2014.    TRIG16CHECK(lold.y);
  2015.    
  2016.    ltmp1 = Exp086(lold.x);
  2017.    SinCos086  (lold.x, &lsinx,  &lcosx);
  2018.    SinCos086  (lold.y, &lsiny,  &lcosy);
  2019.    SinhCosh086(lold.y, &lsinhy, &lcoshy);
  2020.  
  2021.    lnew.x = multiply(ltmp1,       lcosy,   bitshift) +
  2022.             multiply(lsinx,       lcoshy,  bitshift) + longparm->x;
  2023.     
  2024.    lnew.y = multiply(ltmp1,       lsiny,   bitshift);
  2025.             multiply(lcosx,       lsinhy,  bitshift) + longparm->y;
  2026.    lold = lnew;
  2027.    LONGBAILOUT();
  2028.    return(0);
  2029. }
  2030.  
  2031. MarksLambdaFractal() 
  2032. {
  2033.    /* Mark Peterson's variation of "lambda" function */
  2034.  
  2035.    /* Z1 = (C^(exp-1) * Z**2) + C */
  2036.    ltemp.x = ltempsqrx - ltempsqry;
  2037.    ltemp.y = multiply(lold.x ,lold.y ,bitshift)<<1;
  2038.  
  2039.    lnew.x = multiply(lcoefficient.x, ltemp.x, bitshift)
  2040.         - multiply(lcoefficient.y, ltemp.y, bitshift) + longparm->x;
  2041.    lnew.y = multiply(lcoefficient.x, ltemp.y, bitshift)
  2042.         + multiply(lcoefficient.y, ltemp.x, bitshift) + longparm->y;
  2043.  
  2044.    LONGBAILOUT();
  2045.    return(0);
  2046. }
  2047.  
  2048. UnityFractal() 
  2049. {
  2050.    /* brought to you by Mark Peterson - you won't find this in any fractal
  2051.       books unless they saw it here first - Mark invented it! */
  2052.  
  2053.    XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
  2054.    if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delx))
  2055.       return(1);
  2056.    lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
  2057.    lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
  2058.    return(0);
  2059. }
  2060.  
  2061. Mandel4Fractal() 
  2062. {
  2063.    /* By writing this code, Bert has left behind the excuse "don't
  2064.       know what a fractal is, just know how to make'em go fast".
  2065.       Bert is hereby declared a bonafide fractal expert! Supposedly
  2066.       this routine calculates the Mandelbrot/Julia set based on the
  2067.       polynomial z**4 + lambda, but I wouldn't know -- can't follow
  2068.       all that integer math speedup stuff - Tim */
  2069.  
  2070.    /* first, compute (x + iy)**2 */
  2071.    lnew.x  = ltempsqrx - ltempsqry;
  2072.    lnew.y = (multiply(lold.x, lold.y, bitshift) << 1);
  2073.    LONGBAILOUT();
  2074.  
  2075.    /* then, compute ((x + iy)**2)**2 + lambda */
  2076.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  2077.    lnew.y = (multiply(lold.x, lold.y, bitshift) << 1) + longparm->y;
  2078.    LONGBAILOUT();
  2079.    return(0);
  2080. }
  2081.  
  2082. floatZtozpluszpwrFractal() 
  2083. {
  2084.    cpower(&old,5,&new);
  2085.    z_to_the_z(&old,&old);
  2086.    new.x = new.x + old.x +floatparm->x;
  2087.    new.y = new.y + old.y +floatparm->y;
  2088.    FLOATBAILOUT();
  2089.    return(0);
  2090. }
  2091.  
  2092. longZpowerFractal()
  2093. {
  2094.    if(lcpower(&lold,c_exp,&lnew,bitshift))
  2095.       lnew.x = lnew.y = 8L<<bitshift;
  2096.    lnew.x += longparm->x;   
  2097.    lnew.y += longparm->y;   
  2098.    LONGBAILOUT();
  2099.    return(0);
  2100. }
  2101. floatZpowerFractal()
  2102. {
  2103.    cpower(&old,c_exp,&new);
  2104.    new.x += floatparm->x;   
  2105.    new.y += floatparm->y;   
  2106.    FLOATBAILOUT();
  2107.    return(0);
  2108. }
  2109.  
  2110. Barnsley3Fractal() 
  2111. {
  2112.    /* An unnamed Mandelbrot/Julia function from "Fractals
  2113.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  2114.  
  2115.    /* calculate intermediate products */
  2116.    oldxinitx   = multiply(lold.x, lold.x, bitshift);
  2117.    oldyinity   = multiply(lold.y, lold.y, bitshift);
  2118.    oldxinity   = multiply(lold.x, lold.y, bitshift);
  2119.    
  2120.    /* orbit calculation */                  
  2121.    if(lold.x > 0) 
  2122.    {
  2123.       lnew.x = oldxinitx   - oldyinity - fudge;
  2124.       lnew.y = oldxinity << 1;               
  2125.    }
  2126.    else 
  2127.    {
  2128.       lnew.x = oldxinitx - oldyinity - fudge
  2129.            + multiply(longparm->x,lold.x,bitshift);
  2130.       lnew.y = oldxinity <<1;
  2131.  
  2132.       /* This term added by Tim Wegner to make dependent on the 
  2133.          imaginary part of the parameter. (Otherwise Mandelbrot 
  2134.          is uninteresting. */
  2135.       lnew.y += multiply(longparm->y,lold.x,bitshift); 
  2136.    }
  2137.    LONGBAILOUT();
  2138.    return(0);
  2139. }
  2140.  
  2141. SinZsquaredFractal() 
  2142. {
  2143.    /* From Scientific American, July 1989 */
  2144.    /* A Biomorph                          */
  2145.    /* z(n+1) = sin(z(n))+z(n)**2+C        */
  2146.  
  2147.    SinCos086  (lold.x, &lsinx,  &lcosx);
  2148.    SinhCosh086(lold.y, &lsinhy, &lcoshy);
  2149.  
  2150.    lnew.x = multiply(lsinx,        lcoshy,   bitshift);
  2151.    lnew.y = multiply(lcosx,        lsinhy,   bitshift);
  2152.  
  2153.    lnew.x += ltempsqrx - ltempsqry + longparm->x;
  2154.    lnew.y += (multiply(lold.x, lold.y, bitshift) << 1) + longparm->y;
  2155.    LONGBAILOUT();
  2156.    return(0);
  2157. }
  2158.  
  2159. SinZsquaredfpFractal() 
  2160. {
  2161.    /* From Scientific American, July 1989 */
  2162.    /* A Biomorph                          */
  2163.    /* z(n+1) = sin(z(n))+z(n)**2+C        */
  2164.  
  2165.    FPUsincos  (&old.x,&sinx,&cosx);
  2166.    FPUsinhcosh(&old.y,&sinhy,&coshy);
  2167.    new.x = tempsqrx - tempsqry + floatparm->x + sinx*coshy;
  2168.    new.y = 2.0 * old.x * old.y + floatparm->y + cosx*sinhy;
  2169.    FLOATBAILOUT();
  2170.    return(0);
  2171. }
  2172.  
  2173. PopcornFractal()
  2174. {
  2175.    extern int row;
  2176.    tmp = old;
  2177.    tmp.x *= 3.0;
  2178.    tmp.y *= 3.0;
  2179.    FPUsincos(&tmp.x,&sinx,&cosx);
  2180.    FPUsincos(&tmp.y,&siny,&cosy);
  2181.    tmp.x = sinx/cosx + old.x;
  2182.    tmp.y = siny/cosy + old.y;
  2183.    FPUsincos(&tmp.x,&sinx,&cosx);
  2184.    FPUsincos(&tmp.y,&siny,&cosy);
  2185.    new.x = old.x - .05*siny;
  2186.    new.y = old.y - .05*sinx;
  2187.    /*
  2188.    new.x = old.x - .05*sin(old.y+tan(3*old.y));
  2189.    new.y = old.y - .05*sin(old.x+tan(3*old.x));
  2190.    */
  2191.    if(plot == noplot)
  2192.    {
  2193.       plot_orbit(new.x,new.y,1+row%colors);
  2194.       old = new;
  2195.    }
  2196.    else
  2197.       FLOATBAILOUT();
  2198.    return(0);
  2199. }
  2200.  
  2201. LPopcornFractal()
  2202. {
  2203.    static long O5 = (long)(.05*(1L<<16));
  2204.    extern int row;
  2205.    ltemp = lold;
  2206.    ltemp.x *= 3L;
  2207.    ltemp.y *= 3L;
  2208.    LTRIGARG(ltemp.x);
  2209.    LTRIGARG(ltemp.y);
  2210.    SinCos086(ltemp.x,&lsinx,&lcosx);
  2211.    SinCos086(ltemp.y,&lsiny,&lcosy);
  2212.    ltemp.x = divide(lsinx,lcosx,bitshift) + lold.x;
  2213.    ltemp.y = divide(lsiny,lcosy,bitshift) + lold.y;
  2214.    LTRIGARG(ltemp.x);
  2215.    LTRIGARG(ltemp.y);
  2216.    SinCos086(ltemp.x,&lsinx,&lcosx);
  2217.    SinCos086(ltemp.y,&lsiny,&lcosy);
  2218.    lnew.x = lold.x - multiply(O5,lsiny,bitshift);
  2219.    lnew.y = lold.y - multiply(O5,lsinx,bitshift);
  2220.    if(plot == noplot)
  2221.    {
  2222.       iplot_orbit(lnew.x,lnew.y,1+row%colors);
  2223.       lold = lnew;
  2224.    }
  2225.    else
  2226.       LONGBAILOUT();
  2227.    return(0);
  2228. }
  2229.  
  2230. NotImplementedFractal() /* should never get HERE!! */
  2231. {        
  2232.    /*
  2233.    This fractal found in the classic "Table of Even Prime Numbers", the
  2234.    abridged edition, by Vacuous E. Void. */
  2235.    return(1);
  2236. }
  2237.  
  2238. /* -------------------------------------------------------------------- */
  2239. /*        Initialization (once per pixel) routines                        */
  2240. /* -------------------------------------------------------------------- */
  2241.  
  2242. #if 0
  2243. /* this code translated to asm - lives in newton.asm */
  2244. /* transform points with reciprocal function */
  2245. void invertz1(struct complex *z)
  2246. {
  2247.    z->x = dx0[col];
  2248.    z->y = dy0[row];
  2249.    z->x -= f_xcenter; z->y -= f_ycenter;  /* Normalize values to center of circle */
  2250.  
  2251.    tempsqrx = sqr(z->x) + sqr(z->y);  /* Get old radius */
  2252.    if(fabs(tempsqrx) > FLT_MIN)
  2253.       tempsqrx = f_radius / tempsqrx;
  2254.    else
  2255.       tempsqrx = FLT_MAX;   /* a big number, but not TOO big */
  2256.    z->x *= tempsqrx;      z->y *= tempsqrx;      /* Perform inversion */
  2257.    z->x += f_xcenter; z->y += f_ycenter; /* Renormalize */
  2258. }
  2259. #endif
  2260.  
  2261. void long_julia_per_pixel()
  2262. {
  2263.    /* integer julia types */
  2264.    /* lambda */
  2265.    /* barnsleyj1 */  
  2266.    /* barnsleyj2 */  
  2267.    /* sierpinski */
  2268.    if(invert)
  2269.    {
  2270.       /* invert */
  2271.       invertz2(&old);
  2272.       
  2273.       /* watch out for overflow */
  2274.       if(sqr(old.x)+sqr(old.y) >= 127)
  2275.       {
  2276.          old.x = 8;  /* value to bail out in one iteration */
  2277.          old.y = 8;
  2278.       }
  2279.  
  2280.       /* convert to fudged longs */
  2281.       lold.x = old.x*fudge;
  2282.       lold.y = old.y*fudge;
  2283.    }
  2284.    else
  2285.    {
  2286.       lold.x = lx0[col];
  2287.       lold.y = ly0[row];
  2288.    }
  2289. }
  2290.  
  2291. void long_mandel_per_pixel()
  2292. {
  2293.    /* integer mandel types */
  2294.    /* barnsleym1 */
  2295.    /* barnsleym2 */
  2296.    linit.x = lx0[col];
  2297.  
  2298.    if(invert)
  2299.    {
  2300.       /* invert */
  2301.       invertz2(&init);
  2302.  
  2303.       /* watch out for overflow */
  2304.       if(sqr(init.x)+sqr(init.y) >= 127)
  2305.       {
  2306.          init.x = 8;  /* value to bail out in one iteration */
  2307.          init.y = 8;
  2308.       }
  2309.        
  2310.       /* convert to fudged longs */
  2311.       linit.x = init.x*fudge;
  2312.       linit.y = init.y*fudge;
  2313.    }
  2314.  
  2315.    lold.x = linit.x + llambda.x;   /* initial pertubation of parameters set */
  2316.    lold.y = linit.y + llambda.y;
  2317.  
  2318. }
  2319. void julia_per_pixel()
  2320. {
  2321.    /* julia */
  2322.  
  2323.    if(invert)
  2324.    {
  2325.       /* invert */
  2326.       invertz2(&old);
  2327.        
  2328.       /* watch out for overflow */
  2329.       if(bitshift <= 24)
  2330.          if (sqr(old.x)+sqr(old.y) >= 127)
  2331.          {
  2332.             old.x = 8;  /* value to bail out in one iteration */
  2333.             old.y = 8;
  2334.          }
  2335.       if(bitshift >  24)
  2336.          if (sqr(old.x)+sqr(old.y) >= 4.0)
  2337.          {
  2338.             old.x = 2;  /* value to bail out in one iteration */
  2339.             old.y = 2;
  2340.          }
  2341.        
  2342.       /* convert to fudged longs */
  2343.       lold.x = old.x*fudge;
  2344.       lold.y = old.y*fudge;
  2345.    }
  2346.    else
  2347.    {
  2348.       lold.x = lx0[col];
  2349.       lold.y = ly0[row];
  2350.    }
  2351.  
  2352.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2353.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2354. }
  2355.  
  2356. void mandel_per_pixel()
  2357. {
  2358.    /* mandel */
  2359.  
  2360.    if(invert)
  2361.    {
  2362.       invertz2(&init);
  2363.        
  2364.       /* watch out for overflow */
  2365.       if(bitshift <= 24)
  2366.          if (sqr(init.x)+sqr(init.y) >= 127)
  2367.          {
  2368.             init.x = 8;  /* value to bail out in one iteration */
  2369.             init.y = 8;
  2370.          }
  2371.       if(bitshift >  24)
  2372.          if (sqr(init.x)+sqr(init.y) >= 4)
  2373.          {
  2374.             init.x = 2;  /* value to bail out in one iteration */
  2375.             init.y = 2;
  2376.          }
  2377.        
  2378.       /* convert to fudged longs */
  2379.       linit.x = init.x*fudge;
  2380.       linit.y = init.y*fudge;
  2381.    }
  2382.    else
  2383.       linit.x = lx0[col];
  2384.    
  2385.    lold.x = linit.x + llambda.x; /* initial pertubation of parameters set */
  2386.    lold.y = linit.y + llambda.y;
  2387.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2388.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2389. }
  2390.  
  2391.  
  2392. void marksmandel_per_pixel()
  2393. {
  2394.    /* marksmandel */
  2395.  
  2396.    if(invert)
  2397.    {
  2398.       invertz2(&init);
  2399.        
  2400.       /* watch out for overflow */
  2401.       if(sqr(init.x)+sqr(init.y) >= 127)
  2402.       {
  2403.          init.x = 8;  /* value to bail out in one iteration */
  2404.          init.y = 8;
  2405.       }
  2406.        
  2407.       /* convert to fudged longs */
  2408.       linit.x = init.x*fudge;
  2409.       linit.y = init.y*fudge;
  2410.    }
  2411.    else
  2412.       linit.x = lx0[col];
  2413.    
  2414.    lold.x = linit.x + llambda.x; /* initial pertubation of parameters set */
  2415.    lold.y = linit.y + llambda.y;
  2416.  
  2417.    if(c_exp > 3)
  2418.       lcpower(&lold,c_exp-1,&lcoefficient,bitshift);
  2419.    else if(c_exp == 3) {
  2420.       lcoefficient.x = multiply(lold.x, lold.x, bitshift)
  2421.          - multiply(lold.y, lold.y, bitshift);
  2422.       lcoefficient.y = multiply(lold.x, lold.y, bitshift) << 1;
  2423.    }
  2424.    else if(c_exp == 2)
  2425.       lcoefficient = lold;
  2426.    else if(c_exp < 2) {
  2427.       lcoefficient.x = 1L << bitshift;
  2428.       lcoefficient.y = 0L;
  2429.    }
  2430.  
  2431.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2432.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2433. }
  2434.  
  2435.  
  2436. void mandelfp_per_pixel()
  2437. {
  2438.    /* floating point mandelbrot */
  2439.    /* mandelfp */
  2440.  
  2441.    if(invert)
  2442.       invertz2(&init);
  2443.    else
  2444.       init.x = dx0[col];
  2445.    old.x = init.x + parm.x; /* initial pertubation of parameters set */
  2446.    old.y = init.y + parm.y;
  2447.  
  2448.    tempsqrx = sqr(old.x);  /* precalculated value for regular Mandelbrot */
  2449.    tempsqry = sqr(old.y);
  2450. }
  2451.  
  2452. void juliafp_per_pixel()
  2453. {
  2454.    /* floating point julia */
  2455.    /* juliafp */
  2456.    if(invert)
  2457.       invertz2(&old);
  2458.    else
  2459.    {
  2460.      old.x = dx0[col];
  2461.      old.y = dy0[row];
  2462.    }   
  2463.    tempsqrx = sqr(old.x);  /* precalculated value for regular Julia */
  2464.    tempsqry = sqr(old.y);
  2465. }
  2466. void MPCjulia_per_pixel()
  2467. {
  2468.    /* floating point julia */
  2469.    /* juliafp */
  2470.    if(invert)
  2471.       invertz2(&old);
  2472.    else
  2473.    {
  2474.      old.x = dx0[col];
  2475.      old.y = dy0[row];
  2476.    }   
  2477.    mpcold.r = pd2MP(old.x);
  2478.    mpcold.i = pd2MP(old.y);
  2479. }
  2480.  
  2481. void othermandelfp_per_pixel()
  2482. {
  2483.    /* mandelsine */
  2484.    /* mandelcos  */
  2485.    /* mandelexp  */
  2486.    if(invert)
  2487.       invertz2(&init);
  2488.    else
  2489.       init.x = dx0[col];
  2490.    old.x = init.x + parm.x; /* initial pertubation of parameters set */
  2491.    old.y = init.y + parm.y;
  2492. }
  2493.  
  2494. void otherjuliafp_per_pixel()
  2495. {
  2496.    /* lambdasine */
  2497.    /* lambdacos  */
  2498.    /* lambdaexp  */
  2499.    /* mandelsine */
  2500.    /* mandelcos  */
  2501.    /* mandelexp  */
  2502.    if(invert)
  2503.       invertz2(&old);
  2504.    else
  2505.    {
  2506.       old.x = dx0[col];
  2507.       old.y = dy0[row];
  2508.    }   
  2509. }
  2510.  
  2511. /* -------------------------------------------------------------------- */
  2512. /*        Setup (once per fractal image) routines            */
  2513. /* -------------------------------------------------------------------- */
  2514.  
  2515. MandelSetup()        /* Mandelbrot Routine */
  2516. {
  2517.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0 
  2518.           && forcesymmetry == 999 && biomorph == -1)
  2519.       calctype = calcmand; /* the normal case - use CALCMAND */
  2520.    else
  2521.    {
  2522.       /* special case: use the main processing loop */
  2523.       calctype = StandardFractal;
  2524.       longparm = &linit;
  2525.    }
  2526.    return(1);
  2527. }
  2528. JuliaSetup()        /* Julia Routine */
  2529. {
  2530.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0 
  2531.           && forcesymmetry == 999 && biomorph == -1)
  2532.       calctype = calcmand; /* the normal case - use CALCMAND */
  2533.    else
  2534.    {
  2535.       /* special case: use the main processing loop */
  2536.       calctype = StandardFractal;
  2537.       longparm = &llambda;
  2538.    }
  2539.    return(1);
  2540. }
  2541.  
  2542. NewtonSetup()        /* Newton/NewtBasin Routines */
  2543. {
  2544.    int i;
  2545.    extern int basin; 
  2546.    extern int fpu;
  2547.  
  2548.    if(fpu != 0 && debugflag != 1010)
  2549.    {
  2550.       if(fractype == MPNEWTON)
  2551.          fractype = NEWTON;
  2552.       else if(fractype == MPNEWTBASIN)
  2553.          fractype = NEWTBASIN;   
  2554.    }
  2555.  
  2556.    if(fpu == 0 && debugflag != 1010)
  2557.    {
  2558.       if(fractype == NEWTON)
  2559.          fractype = MPNEWTON;
  2560.       else if(fractype == NEWTBASIN)
  2561.          fractype = MPNEWTBASIN;   
  2562.    }
  2563.  
  2564.    /* set up table of roots of 1 along unit circle */
  2565.    degree = (int)parm.x;
  2566.    if(degree < 2)
  2567.       degree = 3;   /* defaults to 3, but 2 is possible */
  2568.    root = 1;
  2569.    
  2570.    /* precalculated values */
  2571.    roverd       = (double)root / (double)degree;
  2572.    d1overd      = (double)(degree - 1) / (double)degree;
  2573.    maxcolor     = 0;
  2574.    threshold    = .3*PI/degree; /* less than half distance between roots */
  2575.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN) {
  2576.       mproverd     = pd2MP(roverd);
  2577.       mpd1overd    = pd2MP(d1overd);
  2578.       mpthreshold  = pd2MP(threshold);
  2579.       mpone        = pd2MP(1.0);
  2580.       }
  2581.  
  2582.    floatmin = FLT_MIN;
  2583.    floatmax = FLT_MAX;
  2584.  
  2585.    basin = 0;
  2586.    if(roots != staticroots) {
  2587.      free(roots);
  2588.      roots = staticroots;
  2589.      }
  2590.          
  2591.    if (fractype==NEWTBASIN)
  2592.    {
  2593.       basin = 1;
  2594.       if(degree > 16)
  2595.       {
  2596.          if((roots=(struct complex *)malloc(degree*sizeof(struct complex)))==NULL)
  2597.          {
  2598.             roots = staticroots;
  2599.             degree = 16;
  2600.          }
  2601.       }
  2602.       else
  2603.          roots = staticroots;
  2604.             
  2605.       /* list of roots to discover where we converged for newtbasin */
  2606.       for(i=0;i<degree;i++)
  2607.       {
  2608.          roots[i].x = cos(i*PI*2.0/(double)degree);
  2609.          roots[i].y = sin(i*PI*2.0/(double)degree);
  2610.       }
  2611.    }
  2612.    else if (fractype==MPNEWTBASIN)
  2613.    {
  2614.       basin = 1;
  2615.       if(degree > 16)
  2616.       {
  2617.         if((MPCroots=(struct MPC *)malloc(degree*sizeof(struct MPC)))==NULL)
  2618.          {
  2619.             MPCroots = (struct MPC *)staticroots;
  2620.             degree = 16;
  2621.          }
  2622.       }
  2623.       else
  2624.          MPCroots = (struct MPC *)staticroots;
  2625.    
  2626.       /* list of roots to discover where we converged for newtbasin */
  2627.       for(i=0;i<degree;i++)
  2628.       {
  2629.          MPCroots[i].r = pd2MP(cos(i*PI*2.0/(double)degree));
  2630.          MPCroots[i].i = pd2MP(sin(i*PI*2.0/(double)degree));
  2631.       }
  2632.    }
  2633.  
  2634.    if (degree%4 == 0)
  2635.       setsymmetry(XYAXIS);
  2636.    else
  2637.       setsymmetry(XAXIS);
  2638.  
  2639.    calctype=StandardFractal;
  2640.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN)
  2641.       setMPfunctions();
  2642.    return(1);
  2643. }
  2644.  
  2645.  
  2646. StandaloneSetup()
  2647. {
  2648.       timer(fractalspecific[fractype].calctype,0);
  2649.       return(0);        /* effectively disable solid-guessing */
  2650. }
  2651.  
  2652. UnitySetup()
  2653. {
  2654.    periodicitycheck = 0;        /* disable periodicity checks */
  2655.    FgOne = (1L << bitshift);
  2656.    FgTwo = FgOne + FgOne;
  2657.    return(1);
  2658. }
  2659. MandelfpSetup()
  2660. {
  2661.    c_exp = param[2];
  2662.    if(fractype==FPMANDELZPOWER && c_exp < 1)
  2663.       c_exp = 1;
  2664.    if(fractype==FPMANDELZPOWER && c_exp & 1) /* odd exponents */
  2665.       setsymmetry(XYAXIS_NOPARM);         
  2666.    
  2667.    floatparm = &init;
  2668.    return(1);
  2669. }
  2670. JuliafpSetup()
  2671. {
  2672.    c_exp = param[2];
  2673.    if(fractype==FPJULIAZPOWER && c_exp < 1)
  2674.       c_exp = 1;
  2675.    if(fractype==FPJULIAZPOWER && c_exp & 1) /* odd exponents */
  2676.       setsymmetry(NOSYM);         
  2677.    floatparm = &parm;
  2678.    return(1);
  2679. }
  2680. MandellongSetup()
  2681. {
  2682.    c_exp = param[2];
  2683.    
  2684.    if(fractype==LMANDELZPOWER && c_exp < 1) 
  2685.       c_exp = 1;
  2686.    if((fractype==MARKSMANDEL   && !(c_exp & 1)) ||
  2687.       (fractype==LMANDELZPOWER && c_exp & 1))
  2688.       setsymmetry(XYAXIS_NOPARM);    /* odd exponents */      
  2689.   longparm = &linit;
  2690.    return(1);
  2691. }
  2692. JulialongSetup()
  2693. {
  2694.    c_exp = param[2];
  2695.    if(fractype==LJULIAZPOWER && c_exp < 1)
  2696.       c_exp = 1;
  2697.    if(fractype==LJULIAZPOWER && c_exp & 1) /* odd exponents */
  2698.       setsymmetry(NOSYM);         
  2699.    longparm = &llambda;
  2700.    return(1);
  2701. }
  2702.  
  2703. MarksJuliaSetup()
  2704. {
  2705.    c_exp = param[2];
  2706.    longparm = &llambda;
  2707.    lold = *longparm;
  2708.    if(c_exp > 2)
  2709.       lcpower(&lold,c_exp,&lcoefficient,bitshift);
  2710.    else if(c_exp == 2)
  2711.    {
  2712.       lcoefficient.x = multiply(lold.x,lold.x,bitshift) - multiply(lold.y,lold.y,bitshift);
  2713.       lcoefficient.y = multiply(lold.x,lold.y,bitshift) << 1;
  2714.    }
  2715.    else if(c_exp < 2)
  2716.       lcoefficient = lold;
  2717.    return(1);
  2718. }
  2719.  
  2720. SierpinskiSetup()
  2721. {
  2722.    /* sierpinski */
  2723.    periodicitycheck = 0;        /* disable periodicity checks */
  2724.    ltemp.x = 1; ltemp.x = ltemp.x << bitshift;    /* ltemp.x = 1 */
  2725.    ltemp.y = ltemp.x >> 1;            /* ltemp.y = .5 */
  2726.    return(1);
  2727. }
  2728.  
  2729. StandardSetup()
  2730. {
  2731.    return(1);
  2732. }
  2733. /* parameter descriptions */
  2734. /* for Mandelbrots */
  2735. static char realz0[] = "Real Portion of Z(0)";
  2736. static char imagz0[] = "Imaginary Portion of Z(0)";
  2737.  
  2738. /* for Julias */
  2739. static char realparm[] = "Real Part of Parameter";
  2740. static char imagparm[] = "Imaginary Part of Parameter";
  2741.  
  2742. /* for Newtons */
  2743. static char newtdegree[] = "Polynomial Degree ( > 2 )";
  2744.  
  2745. /* for MarksMandel/Julia */
  2746. static char exponent[] = "Parameter Exponent (  > 0 )";
  2747.  
  2748. /* for Complex Newton */
  2749. static char realroot[]   = "Real part of Root";
  2750. static char imagroot[]   = "Imag part of Root";
  2751. static char realdegree[] = "Real part of Degree";
  2752. static char imagdegree[] = "Imag part of Degree";
  2753.  
  2754. /* bailout defines */
  2755. #define FTRIGBAILOUT 2500.0
  2756. #define LTRIGBAILOUT   64.0
  2757. #define STDBAILOUT         4.0
  2758.  
  2759. struct fractalspecificstuff fractalspecific[] = 
  2760. {
  2761.    /*
  2762.      fractal name and parameters (text strings)
  2763.      xmin  xmax  ymin  ymax int tojulia   tomandel symmetry 
  2764.    |------|-----|-----|-----|--|--------|---------|--------|
  2765.      orbit fnct     per_pixel fnct  per_image fnct  calctype fcnt    bailout
  2766.    |---------------|---------------|---------------|----------------|-------|
  2767.    */
  2768.  
  2769.    "mandel",      realz0, imagz0,"","",
  2770.    -2.5,  1.5, -1.5,  1.5, 1, JULIA,     NOFRACTAL, MANDELFP, XAXIS_NOPARM, 
  2771.    JuliaFractal,  mandel_per_pixel,MandelSetup,    calcmand,        STDBAILOUT,
  2772.  
  2773.    "julia",       realparm, imagparm,"","",
  2774.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, MANDEL, JULIAFP,  ORIGIN, 
  2775.    JuliaFractal,   julia_per_pixel, JuliaSetup,    calcmand,        STDBAILOUT,
  2776.  
  2777.    "*newtbasin",   newtdegree,"", "","",
  2778.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, MPNEWTBASIN,   NOSYM, 
  2779.    NewtonFractal2, otherjuliafp_per_pixel,  NewtonSetup, StandardFractal,0.0,
  2780.  
  2781.    "lambda",      realparm, imagparm,"","",
  2782.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, MANDELLAMBDA, NOFRACTAL,  NOSYM, 
  2783.    LambdaFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,STDBAILOUT,
  2784.    
  2785.    "*mandel",    realz0, imagz0,"","",
  2786.    -2.5,  1.5, -1.5,  1.5, 0, JULIAFP,   NOFRACTAL, MANDEL,  XAXIS_NOPARM, 
  2787.    JuliafpFractal,mandelfp_per_pixel, MandelfpSetup,StandardFractal,STDBAILOUT,
  2788.    
  2789.    "*newton",      newtdegree,"", "","",
  2790.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, MPNEWTON,   XAXIS, 
  2791.    NewtonFractal2, otherjuliafp_per_pixel,  NewtonSetup, StandardFractal,0.0,
  2792.  
  2793.    "*julia",     realparm, imagparm,"","",
  2794.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, MANDELFP, JULIA,  ORIGIN, 
  2795.    JuliafpFractal, juliafp_per_pixel,  JuliafpSetup,StandardFractal,STDBAILOUT,
  2796.  
  2797.    "plasma",      "Graininess Factor (.1 to 50, default is 2)","","","",
  2798.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM,
  2799.    NULL,           NULL,   StandaloneSetup,      plasma,          0.0,
  2800.    
  2801.    "*lambdasine",  realparm, imagparm,"","",
  2802.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, MANDELSINE, LLAMBDASINE, PI_SYM, 
  2803.    LambdasineFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,FTRIGBAILOUT,
  2804.    
  2805.    "*lambdacos",   realparm, imagparm,"","",
  2806.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, MANDELCOS, LLAMBDACOS,  PI_SYM, 
  2807.    LambdacosineFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,FTRIGBAILOUT,
  2808.  
  2809.    "*lambdaexp",   realparm, imagparm,"","",
  2810.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, MANDELEXP, LLAMBDAEXP,  XAXIS, 
  2811.    LambdaexponentFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,FTRIGBAILOUT,
  2812.  
  2813.    "test",        "(testpt Param #1)","(testpt param #2)","(testpt param #3)", "(testpt param #4)",
  2814.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  2815.    NULL,          NULL,             StandaloneSetup, test,    STDBAILOUT,
  2816.  
  2817.   "sierpinski",  "","","","",
  2818.    -0.9,  1.7, -0.9,  1.7, 1, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  2819.    SierpinskiFractal,long_julia_per_pixel, SierpinskiSetup,StandardFractal,127.0,
  2820.  
  2821.   "barnsleym1",  realz0, imagz0,"","",
  2822.    -2.0,  2.0, -1.5,  1.5, 1, BARNSLEYJ1,NOFRACTAL, NOFRACTAL,  XYAXIS_NOPARM, 
  2823.    Barnsley1Fractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2824.  
  2825.   "barnsleyj1",  realparm, imagparm,"","",
  2826.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, BARNSLEYM1, NOFRACTAL,  ORIGIN, 
  2827.    Barnsley1Fractal,long_julia_per_pixel,JulialongSetup,StandardFractal,STDBAILOUT,
  2828.  
  2829.    "barnsleym2",  realz0, imagz0,"","",
  2830.    -2.0,  2.0, -1.5,  1.5, 1, BARNSLEYJ2,NOFRACTAL, NOFRACTAL,  YAXIS_NOPARM, 
  2831.    Barnsley2Fractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2832.    
  2833.    "barnsleyj2",  realparm, imagparm,"","",
  2834.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, BARNSLEYM2, NOFRACTAL,  ORIGIN, 
  2835.    Barnsley2Fractal,long_julia_per_pixel,JulialongSetup,StandardFractal,STDBAILOUT,
  2836.    
  2837.    "*mandelsine",  realz0, imagz0,"","",
  2838.    -8.0,  8.0, -6.0,  6.0, 0, LAMBDASINE,NOFRACTAL, LMANDELSINE, XYAXIS_NOPARM, 
  2839.    LambdasineFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,FTRIGBAILOUT,
  2840.  
  2841.    "*mandelcos",   realz0, imagz0,"","",
  2842.    -8.0,  8.0, -6.0,  6.0, 0, LAMBDACOS, NOFRACTAL, LMANDELCOS, XYAXIS_NOPARM, 
  2843.    LambdacosineFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,FTRIGBAILOUT,
  2844.  
  2845.    "*mandelexp",   realz0, imagz0,"","",
  2846.    -4.0,  4.0, -3.0,  3.0, 0, LAMBDAEXP, NOFRACTAL, LMANDELEXP,   0, 
  2847.    LambdaexponentFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,FTRIGBAILOUT,
  2848.  
  2849.    "mandellambda",realz0, imagz0,"","",
  2850.    -2.0,  2.0, -1.5,  1.5, 1, LAMBDA,    NOFRACTAL, NOFRACTAL,  XAXIS_NOPARM, 
  2851.    LambdaFractal,mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2852.  
  2853.    "marksmandel", realz0, imagz0, exponent,"",
  2854.    -2.0,  2.0, -1.5,  1.5, 1, MARKSJULIA, NOFRACTAL, NOFRACTAL,  NOSYM,
  2855.    MarksLambdaFractal,marksmandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2856.  
  2857.    "marksjulia", realparm, imagparm, exponent,"",
  2858.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, MARKSMANDEL, NOFRACTAL,   ORIGIN,
  2859.    MarksLambdaFractal,julia_per_pixel,MarksJuliaSetup,StandardFractal,STDBAILOUT,
  2860.  
  2861.    "unity",       "","","","",
  2862.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, NOFRACTAL, NOFRACTAL,   XYAXIS, 
  2863.    UnityFractal, long_julia_per_pixel,UnitySetup,StandardFractal,0.0,
  2864.    
  2865.    "mandel4",      realz0, imagz0,"","",
  2866.    -2.0,  2.0, -1.5,  1.5, 1, JULIA4,     NOFRACTAL, NOFRACTAL,  XAXIS_NOPARM, 
  2867.    Mandel4Fractal,  mandel_per_pixel, MandellongSetup, StandardFractal,  STDBAILOUT,
  2868.  
  2869.    "julia4",       realparm, imagparm,"","",
  2870.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, MANDEL4, NOFRACTAL, ORIGIN, 
  2871.    Mandel4Fractal,   julia_per_pixel, JulialongSetup,StandardFractal,    STDBAILOUT,
  2872.  
  2873.    "ifs",        "","","", "",
  2874.    -8.0,  8.0, -1.0, 11.0, 1, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  2875.    NULL,          NULL,             StandaloneSetup, ifs,    0.0,
  2876.  
  2877.    "ifs3d",        "","","", "",
  2878.    -11.0,  11.0, -11.0, 11.0, 1, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  2879.    NULL,          NULL,      StandaloneSetup, ifs3d,    0.0,
  2880.  
  2881.    "barnsleym3",  realz0, imagz0,"","",
  2882.    -2.0,  2.0, -1.5,  1.5, 1, BARNSLEYJ3,NOFRACTAL, NOFRACTAL,  XAXIS_NOPARM, 
  2883.    Barnsley3Fractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2884.    
  2885.    "barnsleyj3",  realparm, imagparm,"","",
  2886.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, BARNSLEYM3, NOFRACTAL,  XAXIS, 
  2887.    Barnsley3Fractal,long_julia_per_pixel,JulialongSetup,StandardFractal,STDBAILOUT,
  2888.  
  2889.    "demm",    realz0, imagz0,"","",
  2890.    -2.5,  1.5, -1.5,  1.5, 0, DEMJ,   NOFRACTAL, NOFRACTAL,   XAXIS_NOPARM, 
  2891.    NULL,                NULL,         StandaloneSetup, DistanceEstimatorMethod,STDBAILOUT,
  2892.  
  2893.    "demj",     realparm, imagparm,"","",
  2894.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, DEMM, NOFRACTAL,    ORIGIN, 
  2895.    NULL,                NULL,         StandaloneSetup, DistanceEstimatorMethod,STDBAILOUT,
  2896.  
  2897.    "bifurcation",     "", "","","",
  2898.    1.9,  3.0, 0,  1.34, 0, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  2899.    NULL,                NULL,         StandaloneSetup, Bifurcation,STDBAILOUT,
  2900.  
  2901.    "*mandelsinh",  realz0, imagz0,"","",
  2902.    -8.0,  8.0, -6.0,  6.0, 0, LAMBDASINH,NOFRACTAL, LMANDELSINH, XYAXIS_NOPARM, 
  2903.    LambdasinhFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,FTRIGBAILOUT,
  2904.  
  2905.    "*lambdasinh",  realparm, imagparm,"","",
  2906.    -4.0,  4.0, -3.0,  3.0, 0, NOFRACTAL, MANDELSINH, LLAMBDASINH, PI_SYM, 
  2907.    LambdasinhFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,FTRIGBAILOUT,
  2908.  
  2909.    "*mandelcosh",   realz0, imagz0,"","",
  2910.    -8.0,  8.0, -6.0,  6.0, 0, LAMBDACOSH, NOFRACTAL, LMANDELCOSH, XYAXIS_NOPARM, 
  2911.    LambdacoshFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,FTRIGBAILOUT,
  2912.    
  2913.    "*lambdacosh",   realparm, imagparm,"","",
  2914.    -4.0,  4.0, -3.0,  3.0, 0, NOFRACTAL, MANDELCOSH, LLAMBDACOSH, PI_SYM, 
  2915.    LambdacoshFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,FTRIGBAILOUT,
  2916.  
  2917.    "mandelsine",realz0, imagz0,"","",
  2918.    -8.0,  8.0, -6.0,  6.0, 16, LLAMBDASINE, NOFRACTAL, MANDELSINE, XYAXIS_NOPARM, 
  2919.    LongLambdasineFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,LTRIGBAILOUT,
  2920.  
  2921.    "lambdasine",      realparm, imagparm,"","",
  2922.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANDELSINE, LAMBDASINE,PI_SYM, 
  2923.    LongLambdasineFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,LTRIGBAILOUT,
  2924.  
  2925.    "mandelcos",realz0, imagz0,"","",
  2926.    -8.0,  8.0, -6.0,  6.0, 16, LLAMBDACOS, NOFRACTAL, MANDELCOS, XYAXIS_NOPARM, 
  2927.    LongLambdacosFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,LTRIGBAILOUT,
  2928.  
  2929.    "lambdacos",      realparm, imagparm,"","",
  2930.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANDELCOS, LAMBDACOS, PI_SYM, 
  2931.    LongLambdacosFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,LTRIGBAILOUT,
  2932.  
  2933.    "mandelsinh",realz0, imagz0,"","",
  2934.    -8.0,  8.0, -6.0,  6.0, 16, LLAMBDASINH, NOFRACTAL, MANDELSINH, XYAXIS_NOPARM, 
  2935.    LongLambdasinhFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,LTRIGBAILOUT,
  2936.  
  2937.    "lambdasinh",      realparm, imagparm,"","",
  2938.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANDELSINH, LAMBDASINH, ORIGIN, 
  2939.    LongLambdasinhFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,LTRIGBAILOUT,
  2940.  
  2941.    "mandelcosh",realz0, imagz0,"","",
  2942.    -8.0,  8.0, -6.0,  6.0, 16, LLAMBDACOSH,    NOFRACTAL, MANDELCOSH,  XYAXIS_NOPARM, 
  2943.    LongLambdacoshFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,LTRIGBAILOUT,
  2944.  
  2945.    "lambdacosh",      realparm, imagparm,"","",
  2946.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANDELCOSH, LAMBDACOSH, ORIGIN, 
  2947.    LongLambdacoshFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,LTRIGBAILOUT,
  2948.  
  2949.    "mansinzsqrd",      realz0, imagz0,"","",
  2950.    -2.5,  1.5, -1.5,  1.5, 16, LJULSINZSQRD,  NOFRACTAL, FPMANSINZSQRD, XAXIS_NOPARM, 
  2951.    SinZsquaredFractal,mandel_per_pixel,MandellongSetup,StandardFractal, STDBAILOUT,
  2952.  
  2953.    "julsinzsqrd",       realparm, imagparm,"","",
  2954.    -2.0,  2.0, -1.5,  1.5, 16, NOFRACTAL, LMANSINZSQRD, FPJULSINZSQRD,  NOSYM, 
  2955.    SinZsquaredFractal,julia_per_pixel, JulialongSetup,StandardFractal, STDBAILOUT,
  2956.  
  2957.    "*mansinzsqrd",    realz0, imagz0,"","",
  2958.    -2.5,  1.5, -1.5,  1.5, 0, FPJULSINZSQRD,   NOFRACTAL, LMANSINZSQRD, XAXIS_NOPARM, 
  2959.    SinZsquaredfpFractal,mandelfp_per_pixel, MandelfpSetup,StandardFractal, STDBAILOUT,
  2960.  
  2961.    "*julsinzsqrd",     realparm, imagparm,"","",
  2962.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, FPMANSINZSQRD, LJULSINZSQRD,   NOSYM, 
  2963.    SinZsquaredfpFractal, juliafp_per_pixel,  JuliafpSetup,StandardFractal, STDBAILOUT,
  2964.  
  2965.    "mandelexp",realz0, imagz0,"","",
  2966.    -4.0,  4.0, -3.0,  3.0, 16, LLAMBDAEXP,    NOFRACTAL,  MANDELEXP, XAXIS, 
  2967.    LongLambdaexponentFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,LTRIGBAILOUT,
  2968.  
  2969.    "lambdaexp",      realparm, imagparm,"","",
  2970.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANDELEXP, LAMBDAEXP, XAXIS, 
  2971.    LongLambdaexponentFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,LTRIGBAILOUT,
  2972.  
  2973.    "manzpower", realz0, imagz0, exponent,"",
  2974.    -2.0,  2.0, -1.5,  1.5, 1, LJULIAZPOWER, NOFRACTAL, FPMANDELZPOWER,  XAXIS,
  2975.    longZpowerFractal,mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  2976.  
  2977.    "julzpower", realparm, imagparm, exponent,"",
  2978.    -2.0,  2.0, -1.5,  1.5, 1, NOFRACTAL, LMANDELZPOWER, FPJULIAZPOWER,   ORIGIN,
  2979.    longZpowerFractal,julia_per_pixel,JulialongSetup,StandardFractal,STDBAILOUT,
  2980.  
  2981.    "*manzpower",    realz0, imagz0, exponent,"",
  2982.    -2.5,  1.5, -1.5,  1.5, 0, FPJULIAZPOWER,   NOFRACTAL, LMANDELZPOWER,  XAXIS_NOPARM, 
  2983.    floatZpowerFractal,othermandelfp_per_pixel, MandelfpSetup,StandardFractal,STDBAILOUT,
  2984.  
  2985.    "*julzpower",     realparm, imagparm, exponent,"",
  2986.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, FPMANDELZPOWER, LJULIAZPOWER,  ORIGIN, 
  2987.    floatZpowerFractal, otherjuliafp_per_pixel,  JuliafpSetup,StandardFractal,STDBAILOUT,
  2988.  
  2989.    "manzzpwr",    realz0, imagz0, exponent,"",
  2990.    -2.5,  1.5, -1.5,  1.5, 0, FPJULZTOZPLUSZPWR,   NOFRACTAL, NOFRACTAL,  XAXIS_NOPARM, 
  2991.    floatZtozpluszpwrFractal,othermandelfp_per_pixel, MandelfpSetup,StandardFractal,STDBAILOUT,
  2992.  
  2993.    "julzzpwr",     realparm, imagparm, exponent,"",
  2994.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, FPMANZTOZPLUSZPWR, NOFRACTAL,  NOSYM, 
  2995.    floatZtozpluszpwrFractal, otherjuliafp_per_pixel,  JuliafpSetup,StandardFractal,STDBAILOUT,
  2996.  
  2997.    "mansinexp",realz0, imagz0,"","",
  2998.    -8.0,  8.0, -6.0,  6.0, 16, LJULSINEXP,    NOFRACTAL,  FPMANSINEXP, XAXIS_NOPARM, 
  2999.    LongSinexponentFractal,long_mandel_per_pixel,MandellongSetup,StandardFractal,STDBAILOUT,
  3000.  
  3001.    "julsinexp",      realparm, imagparm,"","",
  3002.    -4.0,  4.0, -3.0,  3.0, 16, NOFRACTAL, LMANSINEXP,FPJULSINEXP, NOSYM,
  3003.    LongSinexponentFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,STDBAILOUT,
  3004.  
  3005.    "*mansinexp",   realz0, imagz0,"","",
  3006.    -8.0,  8.0, -6.0,  6.0, 0, FPJULSINEXP, NOFRACTAL, LMANSINEXP,   XAXIS_NOPARM, 
  3007.    FloatSinexponentFractal,othermandelfp_per_pixel,MandelfpSetup,StandardFractal,STDBAILOUT,
  3008.  
  3009.    "*julsinexp",   realparm, imagparm,"","",
  3010.    -4.0,  4.0, -3.0,  3.0, 0, NOFRACTAL, FPMANSINEXP, LJULSINEXP,   NOSYM, 
  3011.    FloatSinexponentFractal,otherjuliafp_per_pixel,JuliafpSetup,StandardFractal,STDBAILOUT,
  3012.  
  3013.    "*popcorn", "", "", "","",
  3014.    -3.0,  3.0, -2.2,  2.2, 0, NOFRACTAL, NOFRACTAL, LPOPCORN,  NOPLOT, 
  3015.    PopcornFractal, otherjuliafp_per_pixel,  JuliafpSetup,StandardFractal,STDBAILOUT,
  3016.  
  3017.    "popcorn", "", "", "","",
  3018.    -3.0,  3.0, -2.2,  2.2, 16, NOFRACTAL, NOFRACTAL, FPPOPCORN,  NOPLOT, 
  3019.    LPopcornFractal,   long_julia_per_pixel, JulialongSetup,  StandardFractal,STDBAILOUT,
  3020.  
  3021.    "*lorenz","Time Step","a","b", "c",
  3022.    -15,  15, 0, 30, 0, NOFRACTAL, NOFRACTAL, LLORENZ,   NOSYM, 
  3023.    NULL,          NULL,             StandaloneSetup, floatlorenz,    0.0,
  3024.  
  3025.    "lorenz","Time Step","a","b", "c",
  3026.    -15,  15, 0, 30, 16, NOFRACTAL, NOFRACTAL, FPLORENZ,   NOSYM, 
  3027.    NULL,          NULL,             StandaloneSetup, longlorenz,    0.0,
  3028.  
  3029.    "lorenz3d","Time Step","a","b", "c",
  3030.    -30.0,  30.0,  -30.0,   30.0, 16, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  3031.    NULL,          NULL,      StandaloneSetup, lorenz3dlong,    0.0,
  3032.  
  3033.    "newton",      newtdegree,"", "","",
  3034.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, NEWTON,   XAXIS, 
  3035.    MPCNewtonFractal, MPCjulia_per_pixel,  NewtonSetup, StandardFractal,0.0,
  3036.  
  3037.    "newtbasin",      newtdegree,"", "","",
  3038.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, NEWTBASIN,   NOSYM, 
  3039.    MPCNewtonFractal, MPCjulia_per_pixel,  NewtonSetup, StandardFractal,0.0,
  3040.  
  3041.    "complexnewton", realdegree, imagdegree, realroot, imagroot,
  3042.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM,
  3043.    ComplexNewton, otherjuliafp_per_pixel,  ComplexNewtonSetup, StandardFractal,0.0,
  3044.  
  3045.    "complexbasin", realdegree, imagdegree, realroot, imagroot,
  3046.    -2.0,  2.0, -1.5,  1.5, 0, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM,
  3047.    ComplexBasin, otherjuliafp_per_pixel,  ComplexNewtonSetup,  StandardFractal, 0.0,
  3048.  
  3049.    NULL, NULL, NULL, NULL, NULL,    /* marks the END of the list */
  3050.    0,  0, 0,  0, 0, NOFRACTAL, NOFRACTAL, NOFRACTAL,   NOSYM, 
  3051.    NULL, NULL, NULL, NULL,0
  3052. };
  3053.