miscfrac.c

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

   if (border <= 0)
      border = 10;

   srand(rseed);
   if (!rflag) ++rseed;

   if (mode == 0) {
      xmax = xdots / 2 + border;  /* Initial box */
      xmin = xdots / 2 - border;
      ymax = ydots / 2 + border;
      ymin = ydots / 2 - border;
   }
   if (mode == 1) {
      xmax = xdots / 2 + border;  /* Initial box */
      xmin = xdots / 2 - border;
      ymin = ydots - 20;
   }
   if (mode == 2) {
      if (xdots>ydots)
         radius = ydots - 5;
      else
         radius = xdots - 5;
   }
   if (resuming) /* restore worklist, if we can't the above will stay in place */
   {
      start_resume();
      if (mode != 2)
         get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
             sizeof(int),&ymin,0);
      else
         get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
             sizeof(int),&radius,0);
      end_resume();
   }

   if (mode==0)
      putcolor(xdots / 2, ydots / 2,RANDOM(colors-1)+1);  /* Seed point */

   if (mode==1)
      for (i=0;i<=xdots;i++)
         putcolor(i,ydots-1,colors-1);

   if (mode==2){
      if (xdots>ydots){
         for (i=0;i<ydots;i++){
            putcolor(xdots/2-ydots/2 , i , colors-1);
            putcolor(xdots/2+ydots/2 , i , colors-1);
            putcolor(xdots/2-ydots/2+i , 0 , colors-1);
            putcolor(xdots/2-ydots/2+i , ydots-1 , colors-1);
         }
      }
      else {
         for (i=0;i<xdots;i++){
            putcolor(0 , ydots/2-xdots/2+i , colors-1);
            putcolor(xdots-1 , ydots/2-xdots/2+i , colors-1);
            putcolor(i , ydots/2-xdots/2 , colors-1);
            putcolor(i , ydots/2+xdots/2 , colors-1);
         }
      }
   }

   while (1)
   {
      /* Release new point on circle just inside the box */

      if (mode==0)
      {
         if (floatflag)
         {
            angle=2*(double)rand()/(RAND_MAX/PI);
            FPUsincos(&angle,&sine,&cosine);
            x = cosine*(xmax-xmin) + xdots;
            y = sine  *(ymax-ymin) + ydots;
         }
         else
         {
            SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
            x = (lcosine*(long)(xmax-xmin) >> 16) + xdots;
            y = (lsine  *(long)(ymax-ymin) >> 16) + ydots;
         }

         x = x >> 1; /* divide by 2 */
         y = y >> 1;
      }
      if (mode==1)
      {
         y=ymin;
         x=RANDOM(xmax-xmin) + (xdots-xmax+xmin)/2;
      }
      if (mode==2)
      {
         if (floatflag)
         {
            angle=2*(double)rand()/(RAND_MAX/PI);
            FPUsincos(&angle,&sine,&cosine);
            x = cosine*radius + xdots;
            y = sine  *radius + ydots;
         }
         else
         {
            SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
            x = (lcosine*(long)(radius) >> 16) + xdots;
            y = (lsine  *(long)(radius) >> 16) + ydots;
         }
         x = x >> 1;
         y = y >> 1;
      }
      /* Loop as long as the point (x,y) is surrounded by color 0 */
      /* on all eight sides                                       */
      while((getcolor(x+1,y+1) == 0) && (getcolor(x+1,y) == 0) &&
          (getcolor(x+1,y-1) == 0) && (getcolor(x  ,y+1) == 0) &&
          (getcolor(x  ,y-1) == 0) && (getcolor(x-1,y+1) == 0) &&
          (getcolor(x-1,y) == 0) && (getcolor(x-1,y-1) == 0))
      {
         /* Erase moving point */
         if (show_orbit)
            putcolor(x,y,0);

         /* Make sure point is inside the box (if mode==0)*/
         if (mode==0){
            if (x==xmax)
               x--;
            else if (x==xmin)
               x++;
            if (y==ymax)
               y--;
            else if (y==ymin)
               y++;
         }

         if (mode==1)
         {
            if (x>xdots-2)
               x--;
            else if (x<1)
               x++;
            if (y<ymin)
               y++;
         }

         /* Take one random step */
         x += RANDOM(3) - 1;
         y += RANDOM(3) - 1;

         /* Check keyboard */
         if ((++plasma_check & 0x7f) == 1)
            if(check_key())
            {
               alloc_resume(20,1);
               if (mode!=2)
                  put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
                      sizeof(int),&ymin,0);
               else
                  put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
                      sizeof(int),&radius,0);

               plasma_check--;
               return 1;
            }

         /* Show the moving point */
         if (show_orbit)
            putcolor(x,y,RANDOM(colors-1)+1);

      }
      putcolor(x,y,RANDOM(colors-1)+1);

      /* Is point too close to the edge? */

      if (mode==0)
      {
         if (((x+border)>xmax) || ((x-border)<xmin)
             || ((y-border)<ymin) || ((y+border)>ymax))
         {
            /* Increase box size, but not past the edge of the screen */
            if (ymin != 1)
            {
               ymin--;
               ymax++;
            }
            if (xmin != 1)
            {
               xmin--;
               xmax++;
            }
            if ((ymin==1) || (xmin==1))
               return 0;
         }
      }
      if (mode==1)
      {
         if (y < ymin+5)
            ymin = y - 5;
         if (ymin<2)
            return 0;
      }
      if (mode==2)
      {
         if (abs(x-xdots/2)<5 && abs(y-ydots/2)<5)
            return 0;

         r = (x-xdots/2)*(x-xdots/2)+(y-ydots/2)*(y-ydots/2);
         fSqrt14(r,f_tmp);
         r = 2 * f_tmp;
         if (r < radius)
            radius = r;
      }
   }
}

/************ standalone engine for "bifurcation" types ***************/

/***************************************************************/
/* The following code now forms a generalised Fractal Engine   */
/* for Bifurcation fractal typeS.  By rights it now belongs in */
/* CALCFRACT.C, but it's easier for me to leave it here !      */

/* Original code by Phil Wilson, hacked around by Kev Allen.   */

/* Besides generalisation, enhancements include Periodicity    */
/* Checking during the plotting phase (AND halfway through the */
/* filter cycle, if possible, to halve calc times), quicker    */
/* floating-point calculations for the standard Verhulst type, */
/* and new bifurcation types (integer bifurcation, f.p & int   */
/* biflambda - the real equivalent of complex Lambda sets -    */
/* and f.p renditions of bifurcations of r*sin(Pi*p), which    */
/* spurred Mitchel Feigenbaum on to discover his Number).      */

/* To add further types, extend the fractalspecific[] array in */
/* usual way, with Bifurcation as the engine, and the name of  */
/* the routine that calculates the next bifurcation generation */
/* as the "orbitcalc" routine in the fractalspecific[] entry.  */

/* Bifurcation "orbitcalc" routines get called once per screen */
/* pixel column.  They should calculate the next generation    */
/* from the doubles Rate & Population (or the longs lRate &    */
/* lPopulation if they use integer math), placing the result   */
/* back in Population (or lPopulation).  They should return 0  */
/* if all is ok, or any non-zero value if calculation bailout  */
/* is desirable (eg in case of errors, or the series tending   */
/* to infinity).		Have fun !		       */
/***************************************************************/

#define DEFAULTFILTER 1000     /* "Beauty of Fractals" recommends using 5000
			       (p.25), but that seems unnecessary. Can
			       override this value with a nonzero param1 */

#define SEED 0.66		/* starting value for population */

static int far *verhulst_array;
unsigned int filter_cycles;
static unsigned int half_time_check;
static long   lPopulation, lRate;
double Population,  Rate;
static int    mono, outside_x;
static long   LPI;

int Bifurcation(void)
{
   unsigned long array_size;
   int row, column;
   column = 0;
   if (resuming)
   {
      start_resume();
      get_resume(sizeof(int),&column,0);
      end_resume();
   }
   array_size =  (iystop + 1) * sizeof(int); /* should be iystop + 1 */
   if ((verhulst_array = (int far *) farmemalloc(array_size)) == NULL)
   {
      static char far msg[]={"Insufficient free memory for calculation."};
      stopmsg(0,msg);
      return(-1);
   }

   LPI = PI * fudge;

   for (row = 0; row <= iystop; row++) /* should be iystop */
      verhulst_array[row] = 0;

   mono = 0;
   if (colors == 2)
      mono = 1;
   if (mono)
   {
      if (inside)
      {
	 outside_x = 0;
	 inside = 1;
      }
      else
	 outside_x = 1;
   }

   filter_cycles = (parm.x <= 0) ? DEFAULTFILTER : parm.x;
   half_time_check = FALSE;
   if (periodicitycheck && maxit < filter_cycles)
   {
      filter_cycles = (filter_cycles - maxit + 1) / 2;
      half_time_check = TRUE;
   }

   if (integerfractal)
      linit.y = ly0[iystop];   /* Y-value of	*/
   else
      init.y = dy0[iystop];   /* bottom pixels */

   while (column <= ixstop)
   {
      if(check_key())
      {
	 farmemfree((char far *)verhulst_array);
	 alloc_resume(10,1);
	 put_resume(sizeof(int),&column,0);
	 return(-1);
      }
      if (integerfractal)
	 lRate = lx0[column];
      else
	 Rate = dx0[column];
      verhulst();	 /* calculate array once per column */

      for (row = iystop; row >= 0; row--) /* should be iystop & >=0 */
      {
	 int color;
	 color = verhulst_array[row];
	 if(color && mono)
	    color = inside;
	 else if((!color) && mono)
	    color = outside_x;
	 else if (color>=colors)
	    color = colors-1;
	 verhulst_array[row] = 0;
	 (*plot)(column,row,color); /* was row-1, but that's not right? */
      }
      column++;
   }
   farmemfree((char far *)verhulst_array);
   return(0);
}

static void verhulst()		/* P. F. Verhulst (1845) */
{
   unsigned int pixel_row, counter, errors;

    if (integerfractal)
       lPopulation = (parm.y == 0) ? SEED * fudge : parm.y * fudge;
    else
       Population = (parm.y == 0 ) ? SEED : parm.y;

   errors = overflow = FALSE;

   for (counter=0 ; counter < filter_cycles ; counter++)
   {
      errors = (*(curfractalspecific->orbitcalc))();
      if (errors)
	 return;
   }
   if (half_time_check) /* check for periodicity at half-time */
   {
      Bif_Period_Init();
      for (counter=0 ; counter < maxit ; counter++)
      {
	 errors = (*(curfractalspecific->orbitcalc))();
	 if (errors) return;
	 if (periodicitycheck && Bif_Periodic(counter)) break;
      }
      if (counter >= maxit)   /* if not periodic, go the distance */
      {
	 for (counter=0 ; counter < filter_cycles ; counter++)
	 {
	    errors = (*(curfractalspecific->orbitcalc))();
	    if (errors) return;
	 }
      }
   }

   if (periodicitycheck) Bif_Period_Init();
   for (counter=0 ; counter < maxit ; counter++)
   {
      errors = (*(curfractalspecific->orbitcalc))();
      if (errors) return;

      /* assign population value to Y coordinate in pixels */
      if (integerfractal)
	 pixel_row = iystop - (lPopulation - linit.y) / dely; /* iystop */
      else
	 pixel_row = iystop - (int)((Population - init.y) / deltaY);

      /* if it's visible on the screen, save it in the column array */
      if (pixel_row <= iystop) /* JCO 6/6/92 */
	 verhulst_array[ pixel_row ] ++;
      if (periodicitycheck && Bif_Periodic(counter))
      {
	 if (pixel_row <= iystop) /* JCO 6/6/92 */
	    verhulst_array[ pixel_row ] --;
	 break;
      }
   }
}
static	long	lBif_closenuf, lBif_savedpop;	/* poss future use  */
static	double	 Bif_closenuf,	Bif_savedpop;
static	int	 Bif_savedinc,	Bif_savedand;

static void Bif_Period_Init()
{
   Bif_savedinc = 1;
   Bif_savedand = 1;
   if (integerfractal)
   {
      lBif_savedpop = -1;
      lBif_closenuf = dely / 8;
   }
   else
   {
      Bif_savedpop = -1.0;
      Bif_closenuf = deltaY / 8.0;
   }
}

static int _fastcall Bif_Periodic (time)  /* Bifurcation Population Periodicity Check */
int time;		/* Returns : 1 if periodicity found, else 0 */
{
   if ((time & Bif_savedand) == 0)	/* time to save a new value */
   {
      if (integerfractal) lBif_savedpop = lPopulation;
      else		     Bif_savedpop =  Population;
      if (--Bif_savedinc == 0)
      {
	 Bif_savedand = (Bif_savedand << 1) + 1;
	 Bif_savedinc = 4;
      }
   }
   else 			/* check against an old save */
   {
      if (integerfractal)
      {
	 if (labs(lBif_savedpop-lPopulation) <= lBif_closenuf)
	    return(1);
      }
      else
      {
	 if (fabs(Bif_savedpop-Population) <= Bif_closenuf)
	    return(1);
      }
   }
   return(0);
}

/**********************************************************************/
/*												      */
/* The following are Bifurcation "orbitcalc" routines...              */
/*												      */
/**********************************************************************/
#ifdef XFRACT
int BifurcLambda() /* Used by lyanupov */
  {
    Population = Rate * Population * (1 - Population);
    return (fabs(Population) > BIG);
  }
#endif

/* Modified formulas below to generalize bifurcations. JCO 7/3/92 */

#define LCMPLXtrig0(arg,out) Arg1->l = (arg); ltrig0(); (out)=Arg1->l
#define  CMPLXtrig0(arg,out) Arg1->d = (arg); dtrig0(); (out)=Arg1->d

int BifurcVerhulstTrig()
  {
/*  Population = Pop + Rate * fn(Pop) * (1 - fn(Pop)) */
    tmp.x = Population;
    tmp.y = 0;
    CMPLXtrig0(tmp, tmp);
    Population += Rate * tmp.x * (1 - tmp.x);
    return (fabs(Population) > BIG);
  }

int LongBifurcVerhulstTrig()
  {
#ifndef XFRACT
    ltmp.x = lPopulation;
    ltmp.y = 0;
    LCMPLXtrig0(ltmp, ltmp);
    ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
    lPopulation += multiply(lRate,ltmp.y,bitshift);
    return (overflow);
#endif
  }

int BifurcStewartTrig()
  {
/*  Population = (Rate * fn(Population) * fn(Population)) - 1.0 */
    tmp.x = Population;
    tmp.y = 0;
    CMPLXtrig0(tmp, tmp);
    Population = (Rate * tmp.x * tmp.x) - 1.0;
    return (fabs(Population) > BIG);
  }

int LongBifurcStewartTrig()
  {
#ifndef XFRACT
    ltmp.x = lPopulation;
    ltmp.y = 0;
    LCMPLXtrig0(ltmp, ltmp);
    lPopulation = multiply(ltmp.x,ltmp.x,bitshift);
    lPopulation = multiply(lPopulation,lRate,	   bitshift);
    lPopulation -= fudge;
    return (overflow);
#endif
  }

int BifurcSetTrigPi()
  {
    tmp.x = Population * PI;
    tmp.y = 0;
    CMPLXtrig0(tmp, tmp);
    Population = Rate * tmp.x;
    return (fabs(Population) > BIG);
  }

int LongBifurcSetTrigPi()
  {
#ifndef XFRACT
    ltmp.x = multiply(lPopulation,LPI,bitshift);
    ltmp.y = 0;
    LCMPLXtrig0(ltmp, ltmp);
    lPopulation = multiply(lRate,ltmp.x,bitshift);
    return (overflow);
#endif
  }

int BifurcAddTrigPi()
  {
    tmp.x = Population * PI;
    tmp.y = 0;
    CMPLXtrig0(tmp, tmp);
    Population += Rate * tmp.x;
    return (fabs(Population) > BIG);
  }