miscfrac.c

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

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

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

int LongBifurcLambdaTrig()
  {
#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
  }

#define LCMPLXpwr(arg1,arg2,out)    Arg2->l = (arg1); Arg1->l = (arg2);\
	 lStkPwr(); Arg1++; Arg2++; (out) = Arg2->l

long beta;

int BifurcMay()
  { /* X = (lambda * X) / (1 + X)^beta, from R.May as described in Pickover,
            Computers, Pattern, Chaos, and Beauty, page 153 */
    tmp.x = 1.0 + Population;
    tmp.x = pow(tmp.x, -beta); /* pow in math.h included with mpmath.h */
    Population = (Rate * Population) * tmp.x;
    return (fabs(Population) > BIG);
  }

int LongBifurcMay()
  {
#ifndef XFRACT
    ltmp.x = lPopulation + fudge;
    ltmp.y = 0;
    lparm2.x = beta * fudge;
    LCMPLXpwr(ltmp, lparm2, ltmp);
    lPopulation = multiply(lRate,lPopulation,bitshift);
    lPopulation = divide(lPopulation,ltmp.x,bitshift);
    return (overflow);
#endif
  }

int BifurcMaySetup()
  {

   beta = (long)param[2];
   if(beta < 2)
      beta = 2;
   param[2] = (double)beta;

   timer(0,curfractalspecific->calctype);
   return(0);
  }

/* Here Endeth the Generalised Bifurcation Fractal Engine   */

/* END Phil Wilson's Code (modified slightly by Kev Allen et. al. !) */


/******************* standalone engine for "popcorn" ********************/

int popcorn()	/* subset of std engine */
{
   int start_row;
   start_row = 0;
   if (resuming)
   {
      start_resume();
      get_resume(sizeof(int),&start_row,0);
      end_resume();
   }
   kbdcount=max_kbdcount;
   plot = noplot;
   tempsqrx = ltempsqrx = 0; /* PB added this to cover weird BAILOUTs */
   for (row = start_row; row <= iystop; row++)
   {
      reset_periodicity = 1;
      for (col = 0; col <= ixstop; col++)
      {
	 if (StandardFractal() == -1) /* interrupted */
	 {
	    alloc_resume(10,1);
	    put_resume(sizeof(int),&row,0);
	    return(-1);
	 }
	 reset_periodicity = 0;
      }
   }
   calc_status = 4;
   return(0);
}

/******************* standalone engine for "lyapunov" *********************/
/*** Roy Murphy [76376,721]						***/
/*** revision history:							***/
/*** initial version: Winter '91					***/
/***    Fall '92 integration of Nicholas Wilt's ASM speedups		***/
/***	Jan 93' integration with calcfrac() yielding boundary tracing,	***/
/***	tesseral, and solid guessing, and inversion, inside=nnn		***/
/*** save_release behavior:						***/
/***    1730 & prior: ignores inside=, calcmode='1', (a,b)->(x,y)	***/
/***	1731: other calcmodes and inside=nnn				***/
/***	1732: the infamous axis swap: (b,a)->(x,y),			***/
/***		the order parameter becomes a long int			***/
/**************************************************************************/
int lyaLength, lyaSeedOK;
int lyaRxy[34];

#define WES 1   /* define WES to be 0 to use Nick's lyapunov.obj */
#if WES
extern int lyapunov_cycles(double, double);
#else
extern int lyapunov_cycles(int, double, double, double);
#endif

int lyapunov_cycles_in_c(int, double, double);

lyapunov () {
    double a, b;

    if (check_key()) {
	return -1;
    	}
    overflow=FALSE;
    if (param[1]==1) Population = (1.0+rand())/(2.0+RAND_MAX);
    else if (param[1]==0) {
	if (fabs(Population)>BIG || Population==0 || Population==1)
	    Population = (1.0+rand())/(2.0+RAND_MAX);
	}
    else Population = param[1];
    (*plot)(col, row, 1);
    if (invert) {
	invertz2(&init);
	a = init.y;
	b = init.x;
	}
    else {
	a = dy0[row]+dy1[col];
	b = dx0[col]+dx1[row];
	}
#ifndef XFRACT
    /*  the assembler routines don't work for a & b outside the
	ranges 0 < a < 4 and 0 < b < 4. So, fall back on the C
	routines if part of the image sticks out.
	*/
#if WES
        color=lyapunov_cycles(a, b);
#else
    if (lyaSeedOK && a>0 && b>0 && a<=4 && b<=4)
	color=lyapunov_cycles(filter_cycles, Population, a, b);
    else
	color=lyapunov_cycles_in_c(filter_cycles, a, b);
#endif
#else
    color=lyapunov_cycles_in_c(filter_cycles, a, b);
#endif
    if (inside>0 && color==0)
	color = inside;
    else if (color>=colors)
	color = colors-1;
    (*plot)(col, row, color);
    return color;
}


lya_setup () {
    /* This routine sets up the sequence for forcing the Rate parameter
	to vary between the two values.  It fills the array lyaRxy[] and
	sets lyaLength to the length of the sequence.

	The sequence is coded in the bit pattern in an integer.
	Briefly, the sequence starts with an A the leading zero bits
	are ignored and the remaining bit sequence is decoded.  The
	sequence ends with a B.  Not all possible sequences can be
	represented in this manner, but every possible sequence is
	either represented as itself, as a rotation of one of the
	representable sequences, or as the inverse of a representable
	sequence (swapping 0s and 1s in the array.)  Sequences that
	are the rotation and/or inverses of another sequence will generate
	the same lyapunov exponents.

	A few examples follow:
	    number    sequence
  		0	ab
  		1	aab
  		2	aabb
  		3	aaab
  		4	aabbb
  		5	aabab
  		6	aaabb (this is a duplicate of 4, a rotated inverse)
  		7	aaaab
  		8	aabbbb	etc.
	 */

    long i;
    int t;

    if ((filter_cycles=param[2])==0)
	filter_cycles=maxit/2;
    lyaSeedOK = param[1]>0 && param[1]<=1 && debugflag!=90;
    lyaLength = 1;

    i = param[0];
#ifndef XFRACT
    if (save_release<1732) i &= 0x0FFFFL; /* make it a short to reporduce prior stuff*/
#endif
    lyaRxy[0] = 1;
    for (t=31; t>=0; t--)
	if (i & (1<<t)) break;
    for (; t>=0; t--)
	lyaRxy[lyaLength++] = (i & (1<<t)) != 0;
    lyaRxy[lyaLength++] = 0;
    if (save_release<1732)		/* swap axes prior to 1732 */
	for (t=lyaLength; t>=0; t--)
	    lyaRxy[t] = !lyaRxy[t];
    if (save_release<1731) {		/* ignore inside=, stdcalcmode */
    	stdcalcmode='1';
	if (inside==1) inside = 0;
	}
    if (inside<0) {
        static char far msg[]=
	    {"Sorry, inside options other than inside=nnn are not supported by the lyapunov"};
        stopmsg(0,(char far *)msg);
        inside=1;
	}
    return 1;
}

int lyapunov_cycles_in_c(int filter_cycles, double a, double b) {
    int color, count, i, lnadjust;
    double lyap, total, temp;
    /* e10=22026.4657948  e-10=0.0000453999297625 */

    for (i=0; i<filter_cycles; i++) {
	for (count=0; count<lyaLength; count++) {
	    Rate = lyaRxy[count] ? a : b;
	    if (curfractalspecific->orbitcalc()) {
		overflow = TRUE;
		goto jumpout;
		}
	    }
	}
    total = 1.0;
    lnadjust = 0;
    for (i=0; i < maxit/2; i++) {
	for (count = 0; count < lyaLength; count++) {
	    Rate = lyaRxy[count] ? a : b;
	    if (curfractalspecific->orbitcalc()) {
		overflow = TRUE;
		goto jumpout;
		}
	    temp = fabs(Rate-2.0*Rate*Population);
		if ((total *= (temp))==0) {
		overflow = TRUE;
		goto jumpout;
		}
	    }
	while (total > 22026.4657948) {
	    total *= 0.0000453999297625;
	    lnadjust += 10;
	    }
	while (total < 0.0000453999297625) {
	    total *= 22026.4657948;
	    lnadjust -= 10;
	    }
	}

jumpout:
    if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0)
	color = 0;
    else {
	if (LogFlag)
	lyap = -temp/((double) lyaLength*i);
    else
	lyap = 1 - exp(temp/((double) lyaLength*i));
	color = 1 + (int)(lyap * (colors-1));
	}
    return color;
}


/******************* standalone engine for "cellular" ********************/
/* Originally coded by Ken Shirriff.
   Modified beyond recognition by Jonathan Osuch.
     Original or'd the neighborhood, changed to sum the neighborhood
     Changed prompts and error messages
     Added CA types
     Set the palette to some standard? CA colors
     Changed *cell_array to near and used dstack so put_line and get_line
       could be used all the time
     Made space bar generate next screen
     Increased string/rule size to 16 digits and added CA types 9/20/92
*/

#define BAD_T         1
#define BAD_MEM       2
#define STRING1       3
#define STRING2       4
#define TABLEK        5
#define TYPEKR        6
#define RULELENGTH    7

#define CELLULAR_DONE 10

#ifndef XFRACT
static BYTE *cell_array[2];
#else
static BYTE far *cell_array[2];
#endif

S16 r, k_1, rule_digits;
int lstscreenflag;

void abort_cellular(int err, int t)
{
   int i;
   switch (err)
   {
      case BAD_T:
         {
         char msg[30];
         sprintf(msg,"Bad t=%d, aborting\n", t);
         stopmsg(0,(char far *)msg);
         }
         break;
      case BAD_MEM:
         {
         static char far msg[]={"Insufficient free memory for calculation" };
         stopmsg(0,msg);
         }
         break;
      case STRING1:
         {
         static char far msg[]={"String can be a maximum of 16 digits" };
         stopmsg(0,msg);
         }
         break;
      case STRING2:
         {
         static char far msg[]={"Make string of 0's through  's" };
         msg[27] = k_1 + 48; /* turn into a character value */
         stopmsg(0,msg);
         }
         break;
      case TABLEK:
         {
         static char far msg[]={"Make Rule with 0's through  's" };
         msg[27] = k_1 + 48; /* turn into a character value */
         stopmsg(0,msg);
         }
         break;
      case TYPEKR:
         {
         static char far msg[]={"Type must be 21, 31, 41, 51, 61, 22, 32, \
42, 23, 33, 24, 25, 26, 27" };
         stopmsg(0,msg);
         }
         break;
      case RULELENGTH:
         {
         static char far msg[]={"Rule must be    digits long" };
         i = rule_digits / 10;
         if(i==0)
            msg[14] = rule_digits + 48;
         else {
            msg[13] = i;
            msg[14] = rule_digits % 10 + 48;
         }
         stopmsg(0,msg);
         }
         break;
      case CELLULAR_DONE:
         break;
   }
   if(cell_array[0] != NULL)
#ifndef XFRACT
      cell_array[0] = NULL;
#else
      farmemfree((char far *)cell_array[0]);
#endif
   if(cell_array[1] != NULL)
#ifndef XFRACT
      cell_array[1] = NULL;
#else
      farmemfree((char far *)cell_array[1]);
#endif
}

cellular () {
   S16 start_row;
   S16 filled, notfilled;
   U16 cell_table[32];
   U16 init_string[16];
   U16 kr, k;
   U32 lnnmbr;
   U16 i,j,l;
   S16 t, t2;
   S32 randparam;
   double n;
   char buf[30];

   set_Cellular_palette();

   randparam = param[0];
   lnnmbr = param[3];
   kr = param[2];
   switch (kr) {
     case 21:
     case 31:
     case 41:
     case 51:
     case 61:
     case 22:
     case 32:
     case 42:
     case 23:
     case 33:
     case 24:
     case 25:
     case 26:
     case 27:
        break;
     default:
        abort_cellular(TYPEKR, 0);
        return -1;
        break;
   }

   r = kr % 10; /* Number of nearest neighbors to sum */
   k = kr / 10; /* Number of different states, k=3 has states 0,1,2 */
   k_1 = k - 1; /* Highest state value, k=3 has highest state value of 2 */
   rule_digits = (r * 2 + 1) * k_1 + 1; /* Number of digits in the rule */

   if ((!rflag) && randparam == -1)
       --rseed;
   if (randparam != 0 && randparam != -1) {
      n = param[0];
      sprintf(buf,"%.16g",n); /* # of digits in initial string */
      t = strlen(buf);
      if (t>16 || t <= 0) {
         abort_cellular(STRING1, 0);
         return -1;
      }
      for (i=0;i<16;i++)
         init_string[i] = 0; /* zero the array */
      t2 = (S16) ((16 - t)/2);
      for (i=0;i<t;i++) { /* center initial string in array */
         init_string[i+t2] = buf[i] - 48; /* change character to number */
         if (init_string[i+t2]>k_1) {
            abort_cellular(STRING2, 0);
            return -1;
         }
      }
   }

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

/* generate rule table from parameter 1 */
#ifndef XFRACT
   n = param[1];
#else
   /* gcc can't manage to convert a big double to an unsigned long properly. */
   if (param[1]>0x7fffffff) {
       n = (param[1]-0x7fffffff);
       n += 0x7fffffff;
   } else {
       n = param[1];
   }
#endif
   if (n == 0) { /* calculate a random rule */
      n = rand15()%k;
      for (i=1;i<rule_digits;i++) {
         n *= 10;
         n += rand15()%k;
      }
      param[1] = n;
   }
   sprintf(buf,"%.*g",rule_digits ,n);
   t = strlen(buf);
   if (rule_digits < t || t < 0) { /* leading 0s could make t smaller */
      abort_cellular(RULELENGTH, 0);
      return -1;
   }
   for (i=0;i<rule_digits;i++) /* zero the table */
      cell_table[i] = 0;
   for (i=0;i<t;i++) { /* reverse order */
      cell_table[i] = buf[t-i-1] - 48; /* change character to number */
      if (cell_table[i]>k_1) {
         abort_cellular(TABLEK, 0);
         return -1;
      }
   }


   start_row = 0;
#ifndef XFRACT
  /* two 4096 byte arrays, at present at most 2024 + 1 bytes should be */
  /* needed in each array (max screen width + 1) */
   cell_array[0] = (BYTE *)&dstack[0]; /* dstack is in general.asm */
   cell_array[1] = (BYTE *)&boxy[0]; /* boxy is in general.asm */
#else
   cell_array[0] = (BYTE far *)farmemalloc(ixstop+1);
   cell_array[1] = (BYTE far *)farmemalloc(ixstop+1);
#endif
   if (cell_array[0]==NULL || cell_array[1]==NULL) {
      abort_cellular(BAD_MEM, 0);
      return(-1);
   }

/* nxtscreenflag toggled by space bar in fractint.c, 1 for continuous */
/* 0 to stop on next screen */

   filled = 0;
   notfilled = 1-filled;
   if (resuming && !nxtscreenflag && !lstscreenflag) {
      start_resume();
      get_resume(sizeof(int),&start_row,0);
      end_resume();
      get_line(start_row,0,ixstop,cell_array[filled]);
   }
   else if (nxtscreenflag && !lstscreenflag) {
      start_resume();
      end_resume();
      get_line(iystop,0,ixstop,cell_array[filled]);