fracsubr.c

frasr200的win 版本源码(18.21),使用make文件,使用的vc版本较低,在我的环境下编译有问题! 很不错的分形程序代码!

/*
FRACSUBR.C contains subroutines which belong primarily to CALCFRAC.C and
FRACTALS.C, i.e. which are non-fractal-specific fractal engine subroutines.
*/

#include <stdio.h>
#ifndef XFRACT
#include <stdarg.h>
#else
#include <varargs.h>
#endif
#include <float.h>
#include <sys/types.h>
#include <sys/timeb.h>
#include <stdlib.h>
#include "fractint.h"
#include "fractype.h"
#include "mpmath.h"
#include "prototyp.h"

/* routines in this module	*/

static long   _fastcall fudgetolong(double d);
static double _fastcall fudgetodouble(long l);
static void   _fastcall adjust_to_limits(double);
static void   _fastcall smallest_add(double *);
static int    _fastcall ratio_bad(double,double);
static void   _fastcall plotdorbit(double,double,int);
static int    _fastcall combine_worklist(void);


extern int    calc_status;		/* status of calculations */
extern char far *resume_info;		/* pointer to resume info if allocated */
       int    resume_len;		/* length of resume info */
static int    resume_offset;		/* offset in resume info gets */
extern double plotmx1,plotmx2,plotmy1,plotmy2; /* real->screen conversion */
extern int    orbit_ptr;		/* pointer into save_orbit array */
extern int orbit_delay; 		/* orbit delay value */
extern int far *save_orbit;		/* array to save orbit values */
extern int    orbit_color;		/* XOR color */
extern int    num_worklist;		/* resume worklist for standard engine */
extern int    fractype; 		/* fractal type */
extern char   stdcalcmode;		/* '1', '2', 'g', 'b'       */
extern char   floatflag;		/* floating-point fractals? */
extern int    integerfractal;		/* TRUE if fractal uses integer math */
extern struct workliststuff worklist[MAXCALCWORK];
extern int    sxdots,sydots;		/* # of dots on the physical screen    */
extern int    sxoffs,syoffs;		/* physical top left of logical screen */
extern int    xdots, ydots;		/* # of dots on the logical screen     */
extern int    colors;			/* maximum colors available */
extern long   fudge;			/* fudge factor (2**n) */
extern int    bitshift; 		/* bit shift for fudge */
extern int    inside;			/* inside color: 1=blue     */
extern int    outside;			/* outside color    */
extern double xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* corners */
extern long   xmin, xmax, ymin, ymax, x3rd, y3rd;  /* integer equivs */
extern int    maxit;			/* try this many iterations */
extern int    attractors;		/* number of finite attractors	*/
extern _CMPLX  attr[];		/* finite attractor vals (f.p)	*/
extern _LCMPLX lattr[]; 	/* finite attractor vals (int)	*/
extern int    attrperiod[]; 		/* finite attractor period	*/
extern _CMPLX  old,new;
extern _LCMPLX lold,lnew;
extern double tempsqrx,tempsqry;
extern long   ltempsqrx,ltempsqry;
extern int    xxstart,xxstop;		/* these are same as worklist, */
extern int    yystart,yystop,yybegin;	/* declared as separate items  */
extern int    periodicitycheck;
extern int    basin;
extern int    finattract;
extern int    pixelpi;			/* value of pi in pixels */
extern double closenuff;
extern long   lclosenuff;
extern int    ixstart, ixstop, iystart, iystop;
extern int    color;
extern int    decomp[];
extern double potparam[3];		/* three potential parameters*/
extern int    distest;			/* non-zero if distance estimator */
extern double param[];  		/* parameters */
extern int    invert;			/* non-zero if inversion active */
extern int    biomorph,usr_biomorph;
extern int    debugflag; /* internal use only - you didn't see this */
extern long   creal, cimag;		/* for calcmand */
extern long   delx, dely;		/* screen pixel increments  */
extern long   delx2, dely2;		/* screen pixel increments  */
extern double delxx, delyy;		/* screen pixel increments  */
extern double delxx2, delyy2;		/* screen pixel increments  */
extern long   delmin;			/* for calcfrac/calcmand    */
extern double ddelmin;			/* same as a double	    */
extern int    potflag;			/* continuous potential flag */
extern int    bailout;
extern double rqlim;
extern double  dxsize, dysize;		/* xdots-1, ydots-1	    */
extern int soundflag;
extern int basehertz;

extern long   far *lx0, far *ly0;	/* x, y grid		    */
extern long   far *lx1, far *ly1;	/* adjustment for rotate    */
/* note that lx1 & ly1 values can overflow into sign bit; since     */
/* they're used only to add to lx0/ly0, 2s comp straightens it out  */
extern double far *dx0, far *dy0;	/* floating pt equivs */
extern double far *dx1, far *dy1;

extern char   usr_floatflag;
extern char   usr_stdcalcmode;
extern int    usr_periodicitycheck;
extern int    usr_distest;
extern double zwidth;

extern int neworbittype;

#define FUDGEFACTOR	29	/* fudge all values up by 2**this */
#define FUDGEFACTOR2	24	/* (or maybe this)		  */


void calcfracinit() /* initialize a *pile* of stuff for fractal calculation */
{
   int i;
   double ftemp;

   floatflag = usr_floatflag;
   if (calc_status == 2) /* on resume, ensure floatflag correct */
      if (curfractalspecific->isinteger)
	 floatflag = 0;
      else
	 floatflag = 1;
   /* if floating pt only, set floatflag for TAB screen */
   if (!curfractalspecific->isinteger && curfractalspecific->tofloat == NOFRACTAL)
      floatflag = 1;

init_restart:

   /* the following variables may be forced to a different setting due to
      calc routine constraints;  usr_xxx is what the user last said is wanted,
      xxx is what we actually do in the current situation */
   stdcalcmode	    = usr_stdcalcmode;
   periodicitycheck = usr_periodicitycheck;
   distest	    = usr_distest;
   biomorph	    = usr_biomorph;

   potflag = 0;
   if (potparam[0] != 0.0
     && colors >= 64
     && (curfractalspecific->calctype == StandardFractal
	 || curfractalspecific->calctype == calcmand
	 || curfractalspecific->calctype == calcmandfp)) {
      potflag = 1;
      distest = 0;    /* can't do distest too */
      }

   if (distest)
      floatflag = 1;  /* force floating point for dist est */

   if (floatflag) { /* ensure type matches floatflag */
      if (curfractalspecific->isinteger != 0
	&& curfractalspecific->tofloat != NOFRACTAL)
	 fractype = curfractalspecific->tofloat;
      }
   else {
      if (curfractalspecific->isinteger == 0
	&& curfractalspecific->tofloat != NOFRACTAL)
	 fractype = curfractalspecific->tofloat;
      }
   /* match Julibrot with integer mode of orbit */   
   if(fractype == JULIBROTFP && fractalspecific[neworbittype].isinteger)
   {
      int i;
      if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
         neworbittype = i;
      else
         fractype = JULIBROT;   
   }
   else if(fractype == JULIBROT && fractalspecific[neworbittype].isinteger==0)
   {
      int i;
      if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
         neworbittype = i;
      else
         fractype = JULIBROTFP;   
   }
      
   curfractalspecific = &fractalspecific[fractype];

   integerfractal = curfractalspecific->isinteger;

/*   if (fractype == JULIBROT)
      rqlim = 4;
   else */ if (potflag && potparam[2] != 0.0)
      rqlim = potparam[2];
/* else if (decomp[0] > 0 && decomp[1] > 0)
      rqlim = (double)decomp[1]; */
   else if (bailout) /* user input bailout */
      rqlim = bailout;
   else if (biomorph != -1) /* biomorph benefits from larger bailout */
      rqlim = 100;
   else
      rqlim = curfractalspecific->orbit_bailout;
   if (integerfractal) /* the bailout limit mustn't be too high here */
      if (rqlim > 127.0) rqlim = 127.0;

   if ((curfractalspecific->flags&NOROTATE) != 0) {
      /* ensure min<max and unrotated rectangle */
      if (xxmin > xxmax) { ftemp = xxmax; xxmax = xxmin; xxmin = ftemp; }
      if (yymin > yymax) { ftemp = yymax; yymax = yymin; yymin = ftemp; }
      xx3rd = xxmin; yy3rd = yymin;
      }

   /* set up bitshift for integer math */
   bitshift = FUDGEFACTOR2; /* by default, the smaller shift */
   if (integerfractal > 1)  /* use specific override from table */
      bitshift = integerfractal;
   if (integerfractal == 0) /* float? */
      if ((i = curfractalspecific->tofloat) != NOFRACTAL) /* -> int? */
	 if (fractalspecific[i].isinteger > 1) /* specific shift? */
	    bitshift = fractalspecific[i].isinteger;

/* We want this code if we're using the assembler calcmand */
   if (fractype == MANDEL || fractype == JULIA) { /* adust shift bits if.. */
      if (potflag == 0				  /* not using potential */
	&& (param[0] > -2.0 && param[0] < 2.0)  /* parameters not too large */
	&& (param[1] > -2.0 && param[1] < 2.0)
	&& !invert				  /* and not inverting */
	&& biomorph == -1			  /* and not biomorphing */
	&& rqlim <= 4.0 			  /* and bailout not too high */
	&& (outside > -2 || outside < -5)	  /* and no funny outside stuff */
	&& debugflag != 1234)			  /* and not debugging */
	 bitshift = FUDGEFACTOR;		  /* use the larger bitshift */
      }

   fudge = 1L << bitshift;

   l_at_rad = fudge/32768L;
   f_at_rad = 1.0/32768L;

   /* now setup arrays of real coordinates corresponding to each pixel */

   adjust_to_limits(1.0); /* make sure all corners in valid range */

   delxx  = (xxmax - xx3rd) / dxsize; /* calculate stepsizes */
   delyy  = (yymax - yy3rd) / dysize;
   delxx2 = (xx3rd - xxmin) / dysize;
   delyy2 = (yy3rd - yymin) / dxsize;

  if(fractype != CELLULAR) { /* fudgetolong fails w >10 digits in double */
   creal = fudgetolong(param[0]); /* integer equivs for it all */
   cimag = fudgetolong(param[1]);
   xmin  = fudgetolong(xxmin);
   xmax  = fudgetolong(xxmax);
   x3rd  = fudgetolong(xx3rd);
   ymin  = fudgetolong(yymin);
   ymax  = fudgetolong(yymax);
   y3rd  = fudgetolong(yy3rd);
   delx  = fudgetolong(delxx);
   dely  = fudgetolong(delyy);
   delx2 = fudgetolong(delxx2);
   dely2 = fudgetolong(delyy2);
  }

   if (fractype != PLASMA) { /* skip this if plasma to avoid 3d problems */
      if (integerfractal && !invert) {
	 if (	(delx  == 0 && delxx  != 0.0)
	     || (delx2 == 0 && delxx2 != 0.0)
	     || (dely  == 0 && delyy  != 0.0)
	     || (dely2 == 0 && delyy2 != 0.0) )
	    goto expand_retry;
	 lx0[0] = xmin; 		/* fill up the x, y grids */
	 ly0[0] = ymax;
	 lx1[0] = ly1[0] = 0;
	 for (i = 1; i < xdots; i++ ) {
	    lx0[i] = lx0[i-1] + delx;
	    ly1[i] = ly1[i-1] - dely2;
	    }
	 for (i = 1; i < ydots; i++ ) {
	    ly0[i] = ly0[i-1] - dely;
	    lx1[i] = lx1[i-1] + delx2;
	    }
	 /* past max res?  check corners within 10% of expected */
	 if (	ratio_bad((double)lx0[xdots-1]-xmin,(double)xmax-x3rd)
	     || ratio_bad((double)ly0[ydots-1]-ymax,(double)y3rd-ymax)
	     || ratio_bad((double)lx1[(ydots>>1)-1],((double)x3rd-xmin)/2)
	     || ratio_bad((double)ly1[(xdots>>1)-1],((double)ymin-y3rd)/2) ) {
expand_retry:
	    if (integerfractal		/* integer fractal type? */
	      && curfractalspecific->tofloat != NOFRACTAL)
	       floatflag = 1;		/* switch to floating pt */
	    else
	       adjust_to_limits(2.0);	/* double the size */
	    if (calc_status == 2)	/* due to restore of an old file? */
	       calc_status = 0; 	/*   whatever, it isn't resumable */
	    goto init_restart;
	    }
	 /* re-set corners to match reality */
	 xmax = lx0[xdots-1] + lx1[ydots-1];
	 ymin = ly0[ydots-1] + ly1[xdots-1];
	 x3rd = xmin + lx1[ydots-1];
	 y3rd = ly0[ydots-1];
	 xxmin = fudgetodouble(xmin);
	 xxmax = fudgetodouble(xmax);
	 xx3rd = fudgetodouble(x3rd);
	 yymin = fudgetodouble(ymin);
	 yymax = fudgetodouble(ymax);
	 yy3rd = fudgetodouble(y3rd);
	 }
      else {
	 /* set up dx0 and dy0 analogs of lx0 and ly0 */
	 /* put fractal parameters in doubles */
	 dx0[0] = xxmin;		/* fill up the x, y grids */
	 dy0[0] = yymax;
	 dx1[0] = dy1[0] = 0;
	 for (i = 1; i < xdots; i++ ) {
	    dx0[i] = dx0[i-1] + delxx;
	    dy1[i] = dy1[i-1] - delyy2;
	    }
	 for (i = 1; i < ydots; i++ ) {
	    dy0[i] = dy0[i-1] - delyy;
	    dx1[i] = dx1[i-1] + delxx2;
	    }
	 if (	ratio_bad(dx0[xdots-1]-xxmin,xxmax-xx3rd)
	     || ratio_bad(dy0[ydots-1]-yymax,yy3rd-yymax)
	     || ratio_bad(dx1[ydots-1],xx3rd-xxmin)
	     || ratio_bad(dy1[xdots-1],yymin-yy3rd))
	    goto expand_retry;
	 /* re-set corners to match reality */
	 xxmax = dx0[xdots-1] + dx1[ydots-1];
	 yymin = dy0[ydots-1] + dy1[xdots-1];
	 xx3rd = xxmin + dx1[ydots-1];
	 yy3rd = dy0[ydots-1];
	 }
      }

   /* for periodicity close-enough, and for unity: */
   /*	  min(max(delx,delx2),max(dely,dely2)	   */