calcfrac.c

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

	 if (pot_startdisk() < 0)
	 {
	    pot16bit = 0;	/* startdisk failed or cancelled */
	    stdcalcmode = tmpcalcmode;	/* maybe we can carry on??? */
	 }
   }
   if (stdcalcmode == 'b' && (curfractalspecific->flags & NOTRACE))
      stdcalcmode = '1';
   if (stdcalcmode == 'g' && (curfractalspecific->flags & NOGUESS))
      stdcalcmode = '1';

   /* default setup a new worklist */
   num_worklist = 1;
   worklist[0].xxstart = 0;
   worklist[0].yystart = worklist[0].yybegin = 0;
   worklist[0].xxstop = xdots - 1;
   worklist[0].yystop = ydots - 1;
   worklist[0].pass = worklist[0].sym = 0;
   if (resuming) /* restore worklist, if we can't the above will stay in place */
   {
      start_resume();
      get_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
      end_resume();
   }

   if (distest) /* setup stuff for distance estimator */
   {
      double ftemp,ftemp2;
      dem_delta = sqr(delxx) + sqr(delyy2);
      if ((ftemp = sqr(delyy) + sqr(delxx2)) > dem_delta)
	 dem_delta = ftemp;
      if (distestwidth == 0)
	 distestwidth = 71;
      ftemp = distestwidth;
      dem_delta *= sqr(ftemp)/10000; /* multiply by thickness desired */
      dem_width = ( sqrt( sqr(xxmax-xxmin) + sqr(xx3rd-xxmin) ) * ydots/xdots
	  + sqrt( sqr(yymax-yymin) + sqr(yy3rd-yymin) ) ) / distest;
      ftemp = (rqlim < DEM_BAILOUT) ? DEM_BAILOUT : rqlim;
      ftemp += 3; /* bailout plus just a bit */
      ftemp2 = log(ftemp);
      dem_toobig = sqr(ftemp) * sqr(ftemp2) * 4 / dem_delta;
   }

   while (num_worklist > 0)
   {
      calctype = curfractalspecific->calctype; /* per_image can override */
      symmetry = curfractalspecific->symmetry; /*   calctype & symmetry  */
      plot = putcolor; /* defaults when setsymmetry not called or does nothing */

      /* pull top entry off worklist */
      ixstart = xxstart = worklist[0].xxstart;
      ixstop  = xxstop	= worklist[0].xxstop;
      iystart = yystart = worklist[0].yystart;
      iystop  = yystop	= worklist[0].yystop;
      yybegin  = worklist[0].yybegin;
      workpass = worklist[0].pass;
      worksym  = worklist[0].sym;
      --num_worklist;
      for (i=0; i<num_worklist; ++i)
	 worklist[i] = worklist[i+1];

      calc_status = 1; /* mark as in-progress */

      curfractalspecific->per_image();

      /* some common initialization for escape-time pixel level routines */
      closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
      closenuff /= fudge;
      lclosenuff = closenuff * fudge;	/* "close enough" value */
      kbdcount=max_kbdcount;
      /* savedmask is for calcmand's periodicity checking */
      savedmask = 0xC0000000; /* top 2 bits on */
      tmplong = (delmin >> abs(periodicitycheck)) | 1;
      while (tmplong > 0) /* while top bit not on */
      {
	 tmplong <<= 1;
	 savedmask = (savedmask >> 1) | 0x80000000;
      }

      setsymmetry(symmetry,1);

/* added for testing autologmap() */
      if (!(resuming)&&(abs(LogFlag) ==2))
       {  /* calculate round screen edges to work out best start for logmap */
         LogFlag = ( autologmap() * (LogFlag / abs(LogFlag)));
         SetupLogTable();
       }

      /* call the appropriate escape-time engine */
      switch (stdcalcmode)
      {
      case 't':
	 tesseral();
	 break;
      case 'b':
	 bound_trace_main();
	 break;
      case 'g':
	 solidguess();
	 break;
      default:
	 OneOrTwoPass();
      }

      if (check_key()) /* interrupted? */
	 break;
   }

   if (num_worklist > 0)
   {  /* interrupted, resumable */
      alloc_resume(sizeof(worklist)+10,1);
      put_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
   }
   else
      calc_status = 4; /* completed */
}

static int OneOrTwoPass()
{
   int i;
   totpasses = 1;
   if (stdcalcmode == '2') totpasses = 2;
   if (stdcalcmode == '2' && workpass == 0) /* do 1st pass of two */
   {
      if (StandardCalc(1) == -1)
      {
	 add_worklist(xxstart,xxstop,yystart,yystop,row,0,worksym);
	 return(-1);
      }
      if (num_worklist > 0) /* worklist not empty, defer 2nd pass */
      {
	 add_worklist(xxstart,xxstop,yystart,yystop,yystart,1,worksym);
	 return(0);
      }
      workpass = 1;
      yybegin = yystart;
   }
   /* second or only pass */
   if (StandardCalc(2) == -1)
   {
      i = yystop;
      if (iystop != yystop) /* must be due to symmetry */
	 i -= row - iystart;
      add_worklist(xxstart,xxstop,row,i,row,workpass,worksym);
      return(-1);
   }
   return(0);
}

static int _fastcall StandardCalc(int passnum)
{
   got_status = 0;
   curpass = passnum;
   row = yybegin;
   while (row <= iystop)
   {
      currow = row;
      reset_periodicity = 1;
      col = ixstart;
      while (col <= ixstop)
      {
         if(showdot>0)
            (*plot) (col, row, showdot&(colors-1));

	 /* on 2nd pass of two, skip even pts */
	 if (passnum == 1 || stdcalcmode == '1' || (row&1) != 0 || (col&1) != 0)
	 {
	    if ((*calctype)() == -1) /* StandardFractal(), calcmand() or calcmandfp() */
	       return(-1); /* interrupted */
	    reset_periodicity = 0;
	    if (passnum == 1) /* first pass, copy pixel and bump col */
	    {
	       if ((row&1) == 0 && row < iystop)
	       {
		  (*plot)(col,row+1,color);
		  if ((col&1) == 0 && col < ixstop)
		     (*plot)(col+1,row+1,color);
	       }
	       if ((col&1) == 0 && col < ixstop)
		  (*plot)(++col,row,color);
	    }
	 }
	 ++col;
      }
      if (passnum == 1 && (row&1) == 0)
	 ++row;
      ++row;
   }
   return(0);
}

int calcmand()		/* fast per pixel 1/2/b/g, called with row & col set */
{
   /* setup values from far array to avoid using es reg in calcmand.asm */
   linitx = lx0[col] + lx1[row];
   linity = ly0[row] + ly1[col];
   if (calcmandasm() >= 0)
   {
      if (LogTable /* map color, but not if maxit & adjusted for inside,etc */
      && (realcolor < maxit || (inside < 0 && color == maxit)))
	 color = LogTable[min(color, maxit)];
      if (color >= colors) /* don't use color 0 unless from inside/outside */
	 if (colors < 16)
	    color &= andcolor;
	 else
	    color = ((color - 1) % andcolor) + 1;  /* skip color zero */
      if(debugflag != 470)
         if(color <= 0 && stdcalcmode == 'b' )   /* fix BTM bug */
            color = 1;
      (*plot) (col, row, color);
   }
   return (color);
}

/************************************************************************/
/* added by Wes Loewer - sort of a floating point version of calcmand() */
/* can also handle invert, any rqlim, potflag, zmag, epsilon cross,     */
/* and all the current outside options    -Wes Loewer 11/03/91          */
/************************************************************************/
int calcmandfp()
{
   if(invert)
      invertz2(&init);
   else
   {
      init.y = dy0[row]+dy1[col];
      init.x = dx0[col]+dx1[row];
   }
   if (calcmandfpasm() >= 0)
   {
      if (potflag)
	 color = potential(magnitude, realcolor);
      if (LogTable /* map color, but not if maxit & adjusted for inside,etc */
		  && (realcolor < maxit || (inside < 0 && color == maxit)))
	    color = LogTable[min(color, maxit)];
	 if (color >= colors) /* don't use color 0 unless from inside/outside */
	    if (colors < 16)
	       color &= andcolor;
	    else
	       color = ((color - 1) % andcolor) + 1;  /* skip color zero */
         if(debugflag != 470)
            if(color == 0 && stdcalcmode == 'b' )   /* fix BTM bug */
               color = 1;
      (*plot) (col, row, color);
   }
   return (color);
}

#define green 2
#define yellow 6
int StandardFractal()	/* per pixel 1/2/b/g, called with row & col set */
{
   double tantable[16];
   int hooper;
   double close;
   long lclose;
   int cyclelen = -1;
   int savedcolor;
   int caught_a_cycle;
   int savedand, savedincr;	/* for periodicity checking */
   _LCMPLX lsaved;
   int i, attracted;
   _LCMPLX lat;
   _CMPLX  at;
   _CMPLX deriv;
   int dem_color;
   _CMPLX dem_new;
   close = .01;
   lclose = close*fudge;

   if(inside == STARTRAIL)
   {
      int i;
      for(i=0;i<16;i++)
	 tantable[i] = 0.0;
   }
   else if(inside == EPSCROSS)
   {
      close = .01;
      lclose = close*fudge;
   }
   if (periodicitycheck == 0 || inside == ZMAG || inside == STARTRAIL)
      oldcolor = 32767; 	/* don't check periodicity at all */
   else if (inside == PERIOD)	/* for display-periodicity */
      oldcolor = maxit*4/5;	/* don't check until nearly done */
   else if (reset_periodicity)
      oldcolor = 250;		/* don't check periodicity 1st 250 iterations */

   /* really fractal specific, but we'll leave it here */
   if (!integerfractal)
   {
      if (useinitorbit == 1)
	 saved = initorbit;
      else {
	 saved.x = 0;
	 saved.y = 0;
      }
      init.y = dy0[row] + dy1[col];
      if (distest)
      {
	 rqlim = rqlim_save;		  /* start with regular bailout */
	 if (distest != 1 || colors == 2) /* not doing regular outside colors */
	    if (rqlim < DEM_BAILOUT)	  /* so go straight for dem bailout */
	       rqlim = DEM_BAILOUT;
	 deriv.x = 1;
	 deriv.y = 0;
	 magnitude = 0;
	 dem_color = -1;
      }
   }
   else
   {
      if (useinitorbit == 1)
	 lsaved = linitorbit;
      else {
	 lsaved.x = 0;
	 lsaved.y = 0;
      }
      linit.y = ly0[row] + ly1[col];
   }
   orbit_ptr = 0;
   color = 0;
   if(fractype==JULIAFP || fractype==JULIA)
      color = -1;
   caught_a_cycle = 0;
   if (inside == PERIOD) {
       savedand = 16;		/* begin checking every 16th cycle */
   } else {
       savedand = 1;		/* begin checking every other cycle */
   }
   savedincr = 1;		/* start checking the very first time */

   if (inside <= BOF60 && inside >= BOF61)
   {
      magnitude = lmagnitud = 0;
      min_orbit = 100000.0;
   }
   overflow = 0;		/* reset integer math overflow flag */

   curfractalspecific->per_pixel(); /* initialize the calculations */

   attracted = FALSE;
   while (++color < maxit)
   {
      if(showdot>0)
         (*plot) (col, row, showdot&(colors-1));

      /* calculation of one orbit goes here */
      /* input in "old" -- output in "new" */

      if (distest)
      {
	 double ftemp;
	 /* Distance estimator for points near Mandelbrot set */
	 /* Original code by Phil Wilson, hacked around by PB */
	 /* Algorithms from Peitgen & Saupe, Science of Fractal Images, p.198 */
	 ftemp	 = 2 * (old.x * deriv.x - old.y * deriv.y) + 1;
	 deriv.y = 2 * (old.y * deriv.x + old.x * deriv.y);
	 deriv.x = ftemp;
	 if (sqr(deriv.x)+sqr(deriv.y) > dem_toobig)
	    break;
	 /* if above exit taken, the later test vs dem_delta will place this
		    point on the boundary, because mag(old)<bailout just now */
	 if (curfractalspecific->orbitcalc())
	 {
	    if (dem_color < 0)	      /* note "regular" color for later */
	    {
	       dem_color = color;
	       dem_new = new;
	    }
	    if (rqlim >= DEM_BAILOUT  /* exit if past real bailout */
	    || magnitude >= (rqlim = DEM_BAILOUT) /* reset to real bailout */
		|| magnitude == 0)	 /* exit if type doesn't "floatbailout" */
	       break;
	    old = new;		      /* carry on till past real bailout */
	 }
      }
      else /* the usual case */
	 if (curfractalspecific->orbitcalc() && inside != STARTRAIL)
	    break;
      if (show_orbit)
	 if (!integerfractal)
	    plot_orbit(new.x, new.y, -1);
	 else
	    iplot_orbit(lnew.x, lnew.y, -1);
      if(inside == STARTRAIL)
      {
	 if(0 < color && color < 16)
	 {
	 if (integerfractal)
	 {
	       new.x = lnew.x;
	       new.x /= fudge;
	       new.y = lnew.y;
	       new.y /= fudge;
	    }
	    tantable[color-1] = new.y/(new.x+.000001);
	 }
      }
      else if(inside == EPSCROSS)
      {
	 hooper = 0;
	 if(integerfractal)
	 {
	    if(labs(lnew.x) < lclose)
	    {
	       hooper = 1; /* close to y axis */
	       goto plot_inside;
	    }
	    else if(labs(lnew.y) < lclose)
	    {
	       hooper = 2; /* close to x axis */
	       goto plot_inside;
	    }
	 }
	 else
	 {
	    if(fabs(new.x) < close)
	    {
	       hooper = 1; /* close to y axis */
	       goto plot_inside;
	    }
	    else if(fabs(new.y) < close)
	    {
	       hooper = 2; /* close to x axis */
	       goto plot_inside;
	    }
	 }
      }
      else if (inside <= BOF60 && inside >= BOF61)
      {
	 if (integerfractal)
	 {
	    if (lmagnitud == 0)
	       lmagnitud = lsqr(lnew.x) + lsqr(lnew.y);
	    magnitude = lmagnitud;
	    magnitude = magnitude / fudge;
	 }
	 else
	    if (magnitude == 0.0)