calcfrac.c

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

/*
CALCFRAC.C contains the high level ("engine") code for calculating the
fractal images (well, SOMEBODY had to do it!).
Original author Tim Wegner, but just about ALL the authors have contributed
SOME code to this routine at one time or another, or contributed to one of
the many massive restructurings.
This module is linked as an overlay, use ENTER_OVLY and EXIT_OVLY.
The following modules work very closely with CALCFRAC.C:
  FRACTALS.C	the fractal-specific code for escape-time fractals.
  FRACSUBR.C	assorted subroutines belonging mainly to calcfrac.
  CALCMAND.ASM	fast Mandelbrot/Julia integer implementation
Additional fractal-specific modules are also invoked from CALCFRAC:
  LORENZ.C	engine level and fractal specific code for attractors.
  JB.C		julibrot logic
  PARSER.C	formula fractals
  and more
 -------------------------------------------------------------------- */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <float.h>
#include <limits.h>
#ifndef XFRACT
#include <dos.h>
#endif
#include <limits.h>
#include "fractint.h"
#include "fractype.h"
#include "mpmath.h"
#include "targa_lc.h"
#include "prototyp.h"

/* routines in this module	*/

static void perform_worklist(void);
static int  OneOrTwoPass(void);
static int  _fastcall StandardCalc(int);
static int  _fastcall potential(double,int);
static void decomposition(void);
static int  bound_trace_main(void);
static void step_col_row(void);
static int  _fastcall boundary_trace(int,int);
static int  _fastcall calc_xy(int,int);
static int  _fastcall fillseg1(int,int,int,int);
static int  _fastcall fillseg(int,int,int,int);
static void _fastcall reverse_string(BYTE *,BYTE *,int);
static int  solidguess(void);
static int  _fastcall guessrow(int,int,int);
static void _fastcall plotblock(int,int,int,int);
static void _fastcall setsymmetry(int,int);
static int  _fastcall xsym_split(int,int);
static int  _fastcall ysym_split(int,int);
static void set_Plasma_palette(void);
static U16 _fastcall adjust(int xa,int ya,int x,int y,int xb,int yb);
static void _fastcall subDivide(int x1,int y1,int x2,int y2);
static int _fastcall new_subD (int x1,int y1,int x2,int y2, int recur);
static void verhulst(void);
static void Bif_Period_Init(void);
static int  _fastcall Bif_Periodic(int);
static void set_Cellular_palette(void);

/**CJLT new function prototypes: */
static int tesseral(void);
static int _fastcall tesschkcol(int,int,int);
static int _fastcall tesschkrow(int,int,int);
static int _fastcall tesscol(int,int,int);
static int _fastcall tessrow(int,int,int);
/* added for testing autologmap() */
static int autologmap(void);

extern _CMPLX initorbit;
extern char useinitorbit;
_LCMPLX linitorbit;
extern int showdot;
extern unsigned int decoderline[];
extern int overflow;
long lmagnitud, llimit, llimit2, lclosenuff, l16triglim;
_CMPLX init,tmp,old,new,saved;
extern int biomorph,usr_biomorph;
extern _LCMPLX linit;
extern int basin;
extern int cpu;
extern char savename[80];   /* save files using this name */
extern int resave_flag;
extern int started_resaves;
extern int dotmode;
extern int save_release;    /* use this to maintain compatibility */

int color, oldcolor, row, col, passes;
int realcolor;
int iterations, invert;
double f_radius,f_xcenter, f_ycenter; /* for inversion */
extern double far *dx0, far *dy0;
extern double far *dx1, far *dy1;

extern int LogFlag;
extern BYTE far *LogTable;
extern int rangeslen;
extern int far *ranges;

void (_fastcall *plot)(int,int,int) = putcolor;
typedef void (_fastcall *PLOT)(int,int,int);

extern double inversion[];	    /* inversion radius, f_xcenter, f_ycenter */
extern int	xdots, ydots;	    /* coordinates of dots on the screen  */
extern int	sxdots,sydots;
extern int	sxoffs,syoffs;
extern int	colors; 	    /* maximum colors available */
extern int	andcolor;	    /* colors-1 		*/
extern int	inside; 	    /* "inside" color to use    */
extern int	fillcolor; 	    /* "fillcolor" color to use    */
extern int	outside;	    /* "outside" color to use   */
extern int	finattract;
double		min_orbit;	    /* orbit value closest to origin */
int		min_index;	    /* iteration of min_orbit */
extern int	maxit;		    /* try this many iterations */
extern int	fractype;	    /* fractal type */
extern char	stdcalcmode;	    /* '1', '2', 'g', 'b' */
extern int	debugflag;	    /* for debugging purposes */
extern	int	diskvideo;	    /* for disk-video klooges */
extern int	calc_status;	    /* status of calculations */
extern long	calctime;	    /* total calc time for image */

extern int rflag, rseed;
extern int decomp[];
extern int distest,distestwidth;

extern double	param[];	    /* parameters */
extern int	potflag;	    /* potential enabled? */
extern double	potparam[];	    /* potential parameters */
extern int	pot16bit;	    /* store 16bit continuous potential */
extern long	far *lx0, far *ly0; /* X, Y points */
extern long	far *lx1, far *ly1; /* X, Y points */
extern long	fudge;		    /* fudge factor (2**n) */
extern int	bitshift;	    /* bit shift for fudge */
extern long	delmin; 	    /* for periodicity checking */

extern double xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* corners */
extern long   xmin, xmax, ymin, ymax;		   /* integer equivs */
extern long   delx,dely;			   /* X, Y increments */
extern double delxx,delxx2,delyy,delyy2;
double deltaX, deltaY;
double magnitude, rqlim, rqlim2, rqlim_save;
extern _CMPLX parm,parm2;
int (*calctype)();
double closenuff;
int pixelpi; /* value of pi in pixels */
unsigned long lm;		/* magnitude limit (CALCMAND) */
extern long linitx,linity;	/* in calcmand */
extern unsigned long savedmask; /* in calcmand */
extern int TranspSymmetry;

/* ORBIT variables */
int	show_orbit;			/* flag to turn on and off */
int	orbit_ptr;			/* pointer into save_orbit array */
int far *save_orbit;			/* array to save orbit values */
int	orbit_color=15; 		/* XOR color */

int	ixstart, ixstop, iystart, iystop;	/* start, stop here */
int	symmetry;	   /* symmetry flag */
int	reset_periodicity; /* nonzero if escape time pixel rtn to reset */
int	kbdcount, max_kbdcount;    /* avoids checking keyboard too often */

extern	int	integerfractal; 	/* TRUE if fractal uses integer math */

char far *resume_info = NULL;		/* pointer to resume info if allocated */
int resuming;				/* nonzero if resuming after interrupt */
int num_worklist;			/* resume worklist for standard engine */
struct workliststuff worklist[MAXCALCWORK];
int xxstart,xxstop;			/* these are same as worklist, */
int yystart,yystop,yybegin;		/* declared as separate items  */
int workpass,worksym;			/* for the sake of calcmand    */

extern long timer_interval;		/* timer(...) total */

VOIDFARPTR typespecific_workarea = NULL;

static double dem_delta, dem_width;	/* distance estimator variables */
static double dem_toobig;
#define DEM_BAILOUT 535.5  /* (pb: not sure if this is special or arbitrary) */

/* variables which must be visible for tab_display */
int got_status; /* -1 if not, 0 for 1or2pass, 1 for ssg, 2 for btm, 3 for 3d */
int curpass,totpasses;
int currow,curcol;

/* static vars for solidguess & its subroutines */
static int maxblock,halfblock;
static int guessplot;			/* paint 1st pass row at a time?   */
static int right_guess,bottom_guess;
#define maxyblk 7    /* maxxblk*maxyblk*2 <= 4096, the size of "prefix" */
#define maxxblk 202  /* each maxnblk is oversize by 2 for a "border" */
		     /* maxxblk defn must match fracsubr.c */
/* next has a skip bit for each maxblock unit;
   1st pass sets bit  [1]... off only if block's contents guessed;
   at end of 1st pass [0]... bits are set if any surrounding block not guessed;
   bits are numbered [..][y/16+1][x+1]&(1<<(y&15)) */
extern unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
/* size of next puts a limit of MAXPIXELS pixels across on solid guessing logic */
#ifndef XFRACT
extern char dstack[4096];               /* common temp, two put_line calls */
extern char suffix[4096];             /* cellular, put_line/get_line calls */
extern unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
#else
BYTE dstack[4096];              /* common temp, two put_line calls */
unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
#endif

int nxtscreenflag; /* for cellular next screen generation */
int	attractors;		    /* number of finite attractors  */
_CMPLX	attr[N_ATTR];	    /* finite attractor vals (f.p)  */
_LCMPLX lattr[N_ATTR];	    /* finite attractor vals (int)  */
int    attrperiod[N_ATTR];	    /* period of the finite attractor */

/***** vars for new btm *****/
enum direction {North,East,South,West};
enum direction going_to;
int trail_row, trail_col;

#ifndef sqr
#define sqr(x) ((x)*(x))
#endif

#ifndef lsqr
#define lsqr(x) (multiply((x),(x),bitshift))
#endif

/* -------------------------------------------------------------------- */
/*		These variables are external for speed's sake only      */
/* -------------------------------------------------------------------- */
extern _LCMPLX lold,lnew,lparm,lparm2;	 /* added "lold" */

int periodicitycheck;
extern long ltempsqrx,ltempsqry;
extern double tempsqrx,tempsqry;
extern _LCMPLX ltmp;
extern int display3d;


void calcfrac_overlay() { }	/* for restore_active_ovly */
extern int first_err;   /* flag for math errors */

/******* calcfract - the top level routine for generating an image *******/

int calcfract()
{
   ENTER_OVLY(OVLY_CALCFRAC);

   first_err = 1;
   attractors = 0;	    /* default to no known finite attractors  */
   display3d = 0;
   basin = 0;

   init_misc();  /* set up some variables in parser.c */

   /* following delta values useful only for types with rotation disabled */
   /* currently used only by bifurcation */
   if (integerfractal)
   {
      distest = 0;
      deltaX = (double)lx0[	 1] / fudge - xxmin;
      deltaY = yymax - (double)ly0[	 1] / fudge;
   }
   else
   {
      deltaX = dx0[	 1] - xxmin;
      deltaY = yymax - dy0[	 1];
   }

   parm.x   = param[0];
   parm.y   = param[1];
   parm2.x  = param[2];
   parm2.y  = param[3];

   if (LogFlag && colors < 16) {
      static char far msg[]={"Need at least 16 colors to use logmap"};
      stopmsg(0,msg);
      LogFlag = 0;
      }
   if (LogFlag || rangeslen)
      if (!(LogTable = farmemalloc((long)maxit + 1))) {
	  static char far msg[]={"Insufficient memory for logmap/ranges with this maxiter"};
      stopmsg(0,msg);
      }
      else if (rangeslen) {
	 int i,k,l,m,numval,flip,altern;
	 i = k = l = 0;
	 while (i < rangeslen) {
	    m = flip = 0;
	    altern = 32767;
	    if ((numval = ranges[i++]) < 0) {
	       altern = ranges[i++];	/* sub-range iterations */
	       numval = ranges[i++];
	       }
	    if (numval > maxit || i >= rangeslen)
	       numval = maxit;
	    while (l <= numval)  {
	       LogTable[l++] = k + flip;
	       if (++m >= altern) {
		  flip ^= 1;		/* Alternate colors */
		  m = 0;
		  }
	       }
	    ++k;
	    if (altern != 32767) ++k;
	    }
	 }
      else
	 SetupLogTable();

   lm = 4L << bitshift; 		/* CALCMAND magnitude limit */

   /* ORBIT stuff */
   save_orbit = (int far *)((double huge *)dx0 + 4*MAXPIXELS);
   show_orbit = 0;
   orbit_ptr = 0;
   orbit_color = 15;
   if(colors < 16)
      orbit_color = 1;

   if(inversion[0] != 0.0)
   {
      f_radius	  = inversion[0];
      f_xcenter   = inversion[1];
      f_ycenter   = inversion[2];

      if (inversion[0] == AUTOINVERT)  /*  auto calc radius 1/6 screen */
	 inversion[0] = f_radius = min(fabs(xxmax - xxmin),
	     fabs(yymax - yymin)) / 6.0;

      if (invert < 2 || inversion[1] == AUTOINVERT)  /* xcenter not already set */
      {
	 inversion[1] = f_xcenter = (xxmin + xxmax) / 2.0;
	 if (fabs(f_xcenter) < fabs(xxmax-xxmin) / 100)
	    inversion[1] = f_xcenter = 0.0;
      }

      if (invert < 3 || inversion[2] == AUTOINVERT)  /* ycenter not already set */
      {
	 inversion[2] = f_ycenter = (yymin + yymax) / 2.0;
	 if (fabs(f_ycenter) < fabs(yymax-yymin) / 100)
	    inversion[2] = f_ycenter = 0.0;
      }

      invert = 3; /* so values will not be changed if we come back */
   }

   closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
   closenuff /= fudge;
   rqlim_save = rqlim;
   rqlim2 = sqrt(rqlim);
   if (integerfractal)		/* for integer routines (lambda) */
   {
      lparm.x = parm.x * fudge;    /* real portion of Lambda */
      lparm.y = parm.y * fudge;    /* imaginary portion of Lambda */
      lparm2.x = parm2.x * fudge;  /* real portion of Lambda2 */
      lparm2.y = parm2.y * fudge;  /* imaginary portion of Lambda2 */
      llimit = rqlim * fudge;	   /* stop if magnitude exceeds this */
      if (llimit <= 0) llimit = 0x7fffffff; /* klooge for integer math */
      llimit2 = rqlim2 * fudge;    /* stop if magnitude exceeds this */
      lclosenuff = closenuff * fudge;	/* "close enough" value */
      l16triglim = 8L<<16;	   /* domain limit of fast trig functions */
      linitorbit.x = initorbit.x * fudge;
      linitorbit.y = initorbit.y * fudge;
   }
   resuming = (calc_status == 2);
   if (!resuming) /* free resume_info memory if any is hanging around */
   {
      end_resume();
      if (resave_flag) {
	 updatesavename(savename); /* do the pending increment */
	 resave_flag = started_resaves = 0;
	 }
      calctime = 0;
   }

   if (curfractalspecific->calctype != StandardFractal
       && curfractalspecific->calctype != calcmand
       && curfractalspecific->calctype != calcmandfp
       && curfractalspecific->calctype != lyapunov
       && curfractalspecific->calctype != calcfroth)
   {
      calctype = curfractalspecific->calctype; /* per_image can override */
      symmetry = curfractalspecific->symmetry; /*   calctype & symmetry  */
      plot = putcolor; /* defaults when setsymmetry not called or does nothing */
      iystart = ixstart = yystart = xxstart = yybegin = 0;
      iystop = yystop = ydots -1;
      ixstop = xxstop = xdots -1;
      calc_status = 1; /* mark as in-progress */
      distest = 0; /* only standard escape time engine supports distest */
      /* per_image routine is run here */
      if (curfractalspecific->per_image())
      { /* not a stand-alone */
	 /* next two lines in case periodicity changed */
	 closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
	 closenuff /= fudge;
	 lclosenuff = closenuff * fudge;	/* "close enough" value */
	 setsymmetry(symmetry,0);
	 timer(0,calctype); /* non-standard fractal engine */
      }
      if (check_key())
      {
	 if (calc_status == 1) /* calctype didn't set this itself, */
	    calc_status = 3;   /* so mark it interrupted, non-resumable */
      }
      else
	 calc_status = 4; /* no key, so assume it completed */
   }
   else /* standard escape-time engine */
      timer(0,(int (*)())perform_worklist);
   calctime += timer_interval;

   if(LogTable)
   {
      farmemfree(LogTable);
      LogTable = NULL;
   }
   if(typespecific_workarea)
   {
      free_workarea();
   }

   EXIT_OVLY;
   return((calc_status == 4) ? 0 : -1);
}

/**************** general escape-time engine routines *********************/

static void perform_worklist()
{
   int i;
   long tmplong; /* this temp must be signed */

   if (potflag && pot16bit)
   {
      int tmpcalcmode = stdcalcmode;

      stdcalcmode = '1'; /* force 1 pass */
      if (resuming == 0)