fractals.c

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

   /* then, compute ((x + iy)**2)**2 + lambda */
   new.x  = tempsqrx - tempsqry + floatparm->x;
   new.y =  old.x*old.y*2 + floatparm->y;
   return(floatbailout());
}

floatZtozPluszpwrFractal()
{
   cpower(&old,(int)param[2],&new);
   old = ComplexPower(old,old);
   new.x = new.x + old.x +floatparm->x;
   new.y = new.y + old.y +floatparm->y;
   return(floatbailout());
}

longZpowerFractal()
{
#ifndef XFRACT
   if(lcpower(&lold,c_exp,&lnew,bitshift))
      lnew.x = lnew.y = 8L<<bitshift;
   lnew.x += longparm->x;
   lnew.y += longparm->y;
   return(longbailout());
#endif
}

longCmplxZpowerFractal()
{
#ifndef XFRACT
   _CMPLX x, y;

   x.x = (double)lold.x / fudge;
   x.y = (double)lold.y / fudge;
   y.x = (double)lparm2.x / fudge;
   y.y = (double)lparm2.y / fudge;
   x = ComplexPower(x, y);
   if(fabs(x.x) < fgLimit && fabs(x.y) < fgLimit) {
      lnew.x = (long)(x.x * fudge);
      lnew.y = (long)(x.y * fudge);
   }
   else
      overflow = 1;
   lnew.x += longparm->x;
   lnew.y += longparm->y;
   return(longbailout());
#endif
}

floatZpowerFractal()
{
   cpower(&old,c_exp,&new);
   new.x += floatparm->x;
   new.y += floatparm->y;
   return(floatbailout());
}

floatCmplxZpowerFractal()
{
   new = ComplexPower(old, parm2);
   new.x += floatparm->x;
   new.y += floatparm->y;
   return(floatbailout());
}

Barnsley3Fractal()
{
   /* An unnamed Mandelbrot/Julia function from "Fractals
   Everywhere" by Michael Barnsley, p. 292, example 4.1 */

   /* calculate intermediate products */
   oldxinitx   = multiply(lold.x, lold.x, bitshift);
   oldyinity   = multiply(lold.y, lold.y, bitshift);
   oldxinity   = multiply(lold.x, lold.y, bitshift);

   /* orbit calculation */
   if(lold.x > 0)
   {
      lnew.x = oldxinitx   - oldyinity - fudge;
      lnew.y = oldxinity << 1;
   }
   else
   {
      lnew.x = oldxinitx - oldyinity - fudge
	   + multiply(longparm->x,lold.x,bitshift);
      lnew.y = oldxinity <<1;

      /* This term added by Tim Wegner to make dependent on the
	 imaginary part of the parameter. (Otherwise Mandelbrot
	 is uninteresting. */
      lnew.y += multiply(longparm->y,lold.x,bitshift);
   }
   return(longbailout());
}

Barnsley3FPFractal()
{
   /* An unnamed Mandelbrot/Julia function from "Fractals
   Everywhere" by Michael Barnsley, p. 292, example 4.1 */


   /* calculate intermediate products */
   foldxinitx  = old.x * old.x;
   foldyinity  = old.y * old.y;
   foldxinity  = old.x * old.y;

   /* orbit calculation */
   if(old.x > 0)
   {
      new.x = foldxinitx - foldyinity - 1.0;
      new.y = foldxinity * 2;
   }
   else
   {
      new.x = foldxinitx - foldyinity -1.0 + floatparm->x * old.x;
      new.y = foldxinity * 2;

      /* This term added by Tim Wegner to make dependent on the
	 imaginary part of the parameter. (Otherwise Mandelbrot
	 is uninteresting. */
      new.y += floatparm->y * old.x;
   }
   return(floatbailout());
}

TrigPlusZsquaredFractal()
{
#ifndef XFRACT
   /* From Scientific American, July 1989 */
   /* A Biomorph			  */
   /* z(n+1) = trig(z(n))+z(n)**2+C	  */
   LCMPLXtrig0(lold,lnew);
   lnew.x += ltempsqrx - ltempsqry + longparm->x;
   lnew.y += multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
   return(longbailout());
#endif
}

TrigPlusZsquaredfpFractal()
{
   /* From Scientific American, July 1989 */
   /* A Biomorph			  */
   /* z(n+1) = trig(z(n))+z(n)**2+C	  */

   CMPLXtrig0(old,new);
   new.x += tempsqrx - tempsqry + floatparm->x;
   new.y += 2.0 * old.x * old.y + floatparm->y;
   return(floatbailout());
}

Richard8fpFractal()
{
   /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
   CMPLXtrig0(old,new);
/*   CMPLXtrig1(*floatparm,tmp); */
   new.x += tmp.x;
   new.y += tmp.y;
   return(floatbailout());
}

Richard8Fractal()
{
#ifndef XFRACT
   /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
   LCMPLXtrig0(lold,lnew);
/*   LCMPLXtrig1(*longparm,ltmp); */
   lnew.x += ltmp.x;
   lnew.y += ltmp.y;
   return(longbailout());
#endif
}

PopcornFractal()
{
   extern int row;
   tmp = old;
   tmp.x *= 3.0;
   tmp.y *= 3.0;
   FPUsincos(&tmp.x,&sinx,&cosx);
   FPUsincos(&tmp.y,&siny,&cosy);
   tmp.x = sinx/cosx + old.x;
   tmp.y = siny/cosy + old.y;
   FPUsincos(&tmp.x,&sinx,&cosx);
   FPUsincos(&tmp.y,&siny,&cosy);
   new.x = old.x - parm.x*siny;
   new.y = old.y - parm.x*sinx;
   /*
   new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
   new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
   */
   if(plot == noplot)
   {
      plot_orbit(new.x,new.y,1+row%colors);
      old = new;
   }
   else
   /* FLOATBAILOUT(); */
   /* PB The above line was weird, not what it seems to be!  But, bracketing
	 it or always doing it (either of which seem more likely to be what
	 was intended) changes the image for the worse, so I'm not touching it.
	 Same applies to int form in next routine. */
   /* PB later: recoded inline, still leaving it weird */
      tempsqrx = sqr(new.x);
   tempsqry = sqr(new.y);
   if((magnitude = tempsqrx + tempsqry) >= rqlim) return(1);
   old = new;
   return(0);
}

LPopcornFractal()
{
#ifndef XFRACT
   extern int row;
   ltmp = lold;
   ltmp.x *= 3L;
   ltmp.y *= 3L;
   LTRIGARG(ltmp.x);
   LTRIGARG(ltmp.y);
   SinCos086(ltmp.x,&lsinx,&lcosx);
   SinCos086(ltmp.y,&lsiny,&lcosy);
   ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
   ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
   LTRIGARG(ltmp.x);
   LTRIGARG(ltmp.y);
   SinCos086(ltmp.x,&lsinx,&lcosx);
   SinCos086(ltmp.y,&lsiny,&lcosy);
   lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
   lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
   if(plot == noplot)
   {
      iplot_orbit(lnew.x,lnew.y,1+row%colors);
      lold = lnew;
   }
   else
      LONGBAILOUT();
   /* PB above still the old way, is weird, see notes in FP popcorn case */
   return(0);
#endif
}

int MarksCplxMand(void)
{
   tmp.x = tempsqrx - tempsqry;
   tmp.y = 2*old.x*old.y;
   FPUcplxmul(&tmp, &Coefficient, &new);
   new.x += floatparm->x;
   new.y += floatparm->y;
   return(floatbailout());
}

int SpiderfpFractal(void)
{
   /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
   new.x = tempsqrx - tempsqry + tmp.x;
   new.y = 2 * old.x * old.y + tmp.y;
   tmp.x = tmp.x/2 + new.x;
   tmp.y = tmp.y/2 + new.y;
   return(floatbailout());
}

SpiderFractal(void)
{
#ifndef XFRACT
   /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
   lnew.x  = ltempsqrx - ltempsqry + ltmp.x;
   lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y;
   ltmp.x = (ltmp.x >> 1) + lnew.x;
   ltmp.y = (ltmp.y >> 1) + lnew.y;
   return(longbailout());
#endif
}

TetratefpFractal()
{
   /* Tetrate(XAXIS) { c=z=pixel: z=c^z, |z|<=(P1+3) } */
   new = ComplexPower(*floatparm,old);
   return(floatbailout());
}

ZXTrigPlusZFractal()
{
#ifndef XFRACT
   /* z = (p1*z*trig(z))+p2*z */
   LCMPLXtrig0(lold,ltmp);	    /* ltmp  = trig(old)	     */
   LCMPLXmult(lparm,ltmp,ltmp);      /* ltmp  = p1*trig(old)	      */
   LCMPLXmult(lold,ltmp,ltmp2);      /* ltmp2 = p1*old*trig(old)      */
   LCMPLXmult(lparm2,lold,ltmp);     /* ltmp  = p2*old		      */
   LCMPLXadd(ltmp2,ltmp,lnew);	     /* lnew  = p1*trig(old) + p2*old */
   return(longbailout());
#endif
}

ScottZXTrigPlusZFractal()
{
#ifndef XFRACT
   /* z = (z*trig(z))+z */
   LCMPLXtrig0(lold,ltmp);	    /* ltmp  = trig(old)       */
   LCMPLXmult(lold,ltmp,lnew);	     /* lnew  = old*trig(old)	*/
   LCMPLXadd(lnew,lold,lnew);	     /* lnew  = trig(old) + old */
   return(longbailout());
#endif
}

SkinnerZXTrigSubZFractal()
{
#ifndef XFRACT
   /* z = (z*trig(z))-z */
   LCMPLXtrig0(lold,ltmp);	    /* ltmp  = trig(old)       */
   LCMPLXmult(lold,ltmp,lnew);	     /* lnew  = old*trig(old)	*/
   LCMPLXsub(lnew,lold,lnew);	     /* lnew  = trig(old) - old */
   return(longbailout());
#endif
}

ZXTrigPlusZfpFractal()
{
   /* z = (p1*z*trig(z))+p2*z */
   CMPLXtrig0(old,tmp); 	 /* tmp  = trig(old)		 */
   CMPLXmult(parm,tmp,tmp);	 /* tmp  = p1*trig(old) 	 */
   CMPLXmult(old,tmp,tmp2);	 /* tmp2 = p1*old*trig(old)	 */
   CMPLXmult(parm2,old,tmp);	 /* tmp  = p2*old		 */
   CMPLXadd(tmp2,tmp,new);	 /* new  = p1*trig(old) + p2*old */
   return(floatbailout());
}

ScottZXTrigPlusZfpFractal()
{
   /* z = (z*trig(z))+z */
   CMPLXtrig0(old,tmp); 	/* tmp	= trig(old)	  */
   CMPLXmult(old,tmp,new);	 /* new  = old*trig(old)   */
   CMPLXadd(new,old,new);	 /* new  = trig(old) + old */
   return(floatbailout());
}

SkinnerZXTrigSubZfpFractal()
{
   /* z = (z*trig(z))-z */
   CMPLXtrig0(old,tmp); 	/* tmp	= trig(old)	  */
   CMPLXmult(old,tmp,new);	 /* new  = old*trig(old)   */
   CMPLXsub(new,old,new);	 /* new  = trig(old) - old */
   return(floatbailout());
}

Sqr1overTrigFractal()
{
#ifndef XFRACT
   /* z = sqr(1/trig(z)) */
   LCMPLXtrig0(lold,lold);
   LCMPLXrecip(lold,lold);
   LCMPLXsqr(lold,lnew);
   return(longbailout());
#endif
}

Sqr1overTrigfpFractal()
{
   /* z = sqr(1/trig(z)) */
   CMPLXtrig0(old,old);
   CMPLXrecip(old,old);
   CMPLXsqr(old,new);
   return(floatbailout());
}

TrigPlusTrigFractal()
{
#ifndef XFRACT
   /* z = trig(0,z)*p1+trig1(z)*p2 */
   LCMPLXtrig0(lold,ltmp);
   LCMPLXmult(lparm,ltmp,ltmp);
   LCMPLXtrig1(lold,ltmp2);
   LCMPLXmult(lparm2,ltmp2,lold);
   LCMPLXadd(ltmp,lold,lnew);
   return(longbailout());
#endif
}

TrigPlusTrigfpFractal()
{
   /* z = trig0(z)*p1+trig1(z)*p2 */
   CMPLXtrig0(old,tmp);
   CMPLXmult(parm,tmp,tmp);
   CMPLXtrig1(old,old);
   CMPLXmult(parm2,old,old);
   CMPLXadd(tmp,old,new);
   return(floatbailout());
}

/* The following four fractals are based on the idea of parallel
   or alternate calculations.  The shift is made when the mod
   reaches a given value.  JCO  5/6/92 */

LambdaTrigOrTrigFractal()
{
#ifndef XFRACT
   /* z = trig0(z)*p1 if mod(old) < p2.x and
          trig1(z)*p1 if mod(old) >= p2.x */
   if ((LCMPLXmod(lold)) < lparm2.x){
     LCMPLXtrig0(lold,ltmp);
     LCMPLXmult(*longparm,ltmp,lnew);}
   else{
     LCMPLXtrig1(lold,ltmp);
     LCMPLXmult(*longparm,ltmp,lnew);}
   return(longbailout());
#endif
}

LambdaTrigOrTrigfpFractal()
{
   /* z = trig0(z)*p1 if mod(old) < p2.x and
          trig1(z)*p1 if mod(old) >= p2.x */
   if (CMPLXmod(old) < parm2.x){
     CMPLXtrig0(old,old);
     FPUcplxmul(floatparm,&old,&new);}
   else{
     CMPLXtrig1(old,old);
     FPUcplxmul(floatparm,&old,&new);}
   return(floatbailout());
}

JuliaTrigOrTrigFractal()
{
#ifndef XFRACT
   /* z = trig0(z)+p1 if mod(old) < p2.x and
          trig1(z)+p1 if mod(old) >= p2.x */
   if (LCMPLXmod(lold) < lparm2.x){
     LCMPLXtrig0(lold,ltmp);
     LCMPLXadd(*longparm,ltmp,lnew);}
   else{
     LCMPLXtrig1(lold,ltmp);
     LCMPLXadd(*longparm,ltmp,lnew);}
   return(longbailout());
#endif
}

JuliaTrigOrTrigfpFractal()
{
   /* z = trig0(z)+p1 if mod(old) < p2.x and
          trig1(z)+p1 if mod(old) >= p2.x */
   if (CMPLXmod(old) < parm2.x){
     CMPLXtrig0(old,old);
     CMPLXadd(*floatparm,old,new);}
   else{
     CMPLXtrig1(old,old);
     CMPLXadd(*floatparm,old,new);}
   return(floatbailout());
}

static int AplusOne, Ap1deg;
static struct MP mpAplusOne, mpAp1deg, mptmpparmy;

int MPCHalleyFractal()
{
#ifndef XFRACT
   /*  X(X^a - 1) = 0, Halley Map */
   /*  a = parm.x,  relaxation coeff. = parm.y,  epsilon = parm2.x  */

int ihal;
struct MPC mpcXtoAlessOne, mpcXtoA;
struct MPC mpcXtoAplusOne; /* a-1, a, a+1 */
struct MPC mpcFX, mpcF1prime, mpcF2prime, mpcHalnumer1;
struct MPC mpcHalnumer2, mpcHaldenom, mpctmp;

   MPOverflow = 0;
   mpcXtoAlessOne.x = mpcold.x;
   mpcXtoAlessOne.y = mpcold.y;
   for(ihal=2; ihal<degree; ihal++) {
     mpctmp.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
     mpctmp.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
     mpcXtoAlessOne.x = mpctmp.x;
     mpcXtoAlessOne.y = mpctmp.y;
   }
   mpcXtoA.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
   mpcXtoA.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
   mpcXtoAplusOne.x = *pMPsub(*pMPmul(mpcXtoA.x,mpcold.x),*pMPmul(mpcXtoA.y,mpcold.y));
   mpcXtoAplusOne.y = *pMPadd(*pMPmul(mpcXtoA.x,mpcold.y),*pMPmul(mpcXtoA.y,mpcold.x));

   mpcFX.x = *pMPsub(mpcXtoAplusOne.x, mpcold.x);
   mpcFX.y = *pMPsub(mpcXtoAplusOne.y, mpcold.y); /* FX = X^(a+1) - X  = F */

   mpcF2prime.x = *pMPmul(mpAp1deg, mpcXtoAlessOne.x); /* mpAp1deg in setup */
   mpcF2prime.y = *pMPmul(mpAp1deg, mpcXtoAlessOne.y);        /* F" */

   mpcF1prime.x = *pMPsub(*pMPmul(mpAplusOne, mpcXtoA.x), mpone);
   mpcF1prime.y = *pMPmul(mpAplusOne, mpcXtoA.y);                   /*  F'  */

   mpctmp.x = *pMPsub(*pMPmul(mpcF2prime.x,mpcFX.x),*pMPmul(mpcF2prime.y,mpcFX.y));