📄 fractals.c
字号:
mpctmp.y = *pMPadd(*pMPmul(mpcF2prime.x,mpcFX.y),*pMPmul(mpcF2prime.y,mpcFX.x));
/* F * F" */
mpcHaldenom.x = *pMPadd(mpcF1prime.x, mpcF1prime.x);
mpcHaldenom.y = *pMPadd(mpcF1prime.y, mpcF1prime.y); /* 2 * F' */
mpcHalnumer1 = MPCdiv(mpctmp, mpcHaldenom); /* F"F/2F' */
mpctmp.x = *pMPsub(mpcF1prime.x, mpcHalnumer1.x);
mpctmp.y = *pMPsub(mpcF1prime.y, mpcHalnumer1.y); /* F' - F"F/2F' */
mpcHalnumer2 = MPCdiv(mpcFX, mpctmp);
mpctmp.x = *pMPmul(mptmpparmy,mpcHalnumer2.x); /* mptmpparmy is */
mpctmp.y = *pMPmul(mptmpparmy,mpcHalnumer2.y); /* relaxation coef. */
mpcnew.x = *pMPsub(mpcold.x, mpctmp.x);
mpcnew.y = *pMPsub(mpcold.y, mpctmp.y);
new.x = *pMP2d(mpcnew.x);
new.y = *pMP2d(mpcnew.y);
return(MPCHalleybailout()||MPOverflow);
#endif
}
HalleyFractal()
{
/* X(X^a - 1) = 0, Halley Map */
/* a = parm.x = degree, relaxation coeff. = parm.y, epsilon = parm2.x */
int ihal;
_CMPLX XtoAlessOne, XtoA, XtoAplusOne; /* a-1, a, a+1 */
_CMPLX FX, F1prime, F2prime, Halnumer1, Halnumer2, Haldenom;
XtoAlessOne = old;
for(ihal=2; ihal<degree; ihal++) {
FPUcplxmul(&old, &XtoAlessOne, &XtoAlessOne);
}
FPUcplxmul(&old, &XtoAlessOne, &XtoA);
FPUcplxmul(&old, &XtoA, &XtoAplusOne);
CMPLXsub(XtoAplusOne, old, FX); /* FX = X^(a+1) - X = F */
F2prime.x = Ap1deg * XtoAlessOne.x; /* Ap1deg in setup */
F2prime.y = Ap1deg * XtoAlessOne.y; /* F" */
F1prime.x = AplusOne * XtoA.x - 1.0;
F1prime.y = AplusOne * XtoA.y; /* F' */
FPUcplxmul(&F2prime, &FX, &Halnumer1); /* F * F" */
Haldenom.x = F1prime.x + F1prime.x;
Haldenom.y = F1prime.y + F1prime.y; /* 2 * F' */
FPUcplxdiv(&Halnumer1, &Haldenom, &Halnumer1); /* F"F/2F' */
CMPLXsub(F1prime, Halnumer1, Halnumer2); /* F' - F"F/2F' */
FPUcplxdiv(&FX, &Halnumer2, &Halnumer2);
new.x = old.x - (parm.y * Halnumer2.x); /* parm.y is relaxation coef. */
new.y = old.y - (parm.y * Halnumer2.y);
return(Halleybailout());
}
LongPhoenixFractal()
{
#ifndef XFRACT
/* z(n+1) = z(n)^2 + p + qy(n), y(n+1) = z(n) */
ltmp.x = multiply(lold.x, lold.y, bitshift);
lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(longparm->y,ltmp2.x,bitshift);
lnew.y = (ltmp.x + ltmp.x) + multiply(longparm->y,ltmp2.y,bitshift);
ltmp2 = lold; /* set ltmp2 to Y value */
return(longbailout());
#endif
}
PhoenixFractal()
{
/* z(n+1) = z(n)^2 + p + qy(n), y(n+1) = z(n) */
tmp.x = old.x * old.y;
new.x = tempsqrx - tempsqry + floatparm->x + (floatparm->y * tmp2.x);
new.y = (tmp.x + tmp.x) + (floatparm->y * tmp2.y);
tmp2 = old; /* set tmp2 to Y value */
return(floatbailout());
}
LongPhoenixPlusFractal()
{
#ifndef XFRACT
/* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
int i;
_LCMPLX loldplus, lnewminus;
loldplus = lold;
ltmp = lold;
for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
}
loldplus.x += longparm->x;
LCMPLXmult(ltmp, loldplus, lnewminus);
lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
ltmp2 = lold; /* set ltmp2 to Y value */
return(longbailout());
#endif
}
PhoenixPlusFractal()
{
/* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n), y(n+1) = z(n) */
int i;
_CMPLX oldplus, newminus;
oldplus = old;
tmp = old;
for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
}
oldplus.x += floatparm->x;
FPUcplxmul(&tmp, &oldplus, &newminus);
new.x = newminus.x + (floatparm->y * tmp2.x);
new.y = newminus.y + (floatparm->y * tmp2.y);
tmp2 = old; /* set tmp2 to Y value */
return(floatbailout());
}
LongPhoenixMinusFractal()
{
#ifndef XFRACT
/* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
int i;
_LCMPLX loldsqr, lnewminus;
LCMPLXmult(lold,lold,loldsqr);
ltmp = lold;
for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
}
loldsqr.x += longparm->x;
LCMPLXmult(ltmp, loldsqr, lnewminus);
lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
ltmp2 = lold; /* set ltmp2 to Y value */
return(longbailout());
#endif
}
PhoenixMinusFractal()
{
/* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n), y(n+1) = z(n) */
int i;
_CMPLX oldsqr, newminus;
FPUcplxmul(&old, &old, &oldsqr);
tmp = old;
for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
}
oldsqr.x += floatparm->x;
FPUcplxmul(&tmp, &oldsqr, &newminus);
new.x = newminus.x + (floatparm->y * tmp2.x);
new.y = newminus.y + (floatparm->y * tmp2.y);
tmp2 = old; /* set tmp2 to Y value */
return(floatbailout());
}
ScottTrigPlusTrigFractal()
{
#ifndef XFRACT
/* z = trig0(z)+trig1(z) */
LCMPLXtrig0(lold,ltmp);
LCMPLXtrig1(lold,lold);
LCMPLXadd(ltmp,lold,lnew);
return(longbailout());
#endif
}
ScottTrigPlusTrigfpFractal()
{
/* z = trig0(z)+trig1(z) */
CMPLXtrig0(old,tmp);
CMPLXtrig1(old,tmp2);
CMPLXadd(tmp,tmp2,new);
return(floatbailout());
}
SkinnerTrigSubTrigFractal()
{
#ifndef XFRACT
/* z = trig(0,z)-trig1(z) */
LCMPLXtrig0(lold,ltmp);
LCMPLXtrig1(lold,ltmp2);
LCMPLXsub(ltmp,ltmp2,lnew);
return(longbailout());
#endif
}
SkinnerTrigSubTrigfpFractal()
{
/* z = trig0(z)-trig1(z) */
CMPLXtrig0(old,tmp);
CMPLXtrig1(old,tmp2);
CMPLXsub(tmp,tmp2,new);
return(floatbailout());
}
TrigXTrigfpFractal()
{
/* z = trig0(z)*trig1(z) */
CMPLXtrig0(old,tmp);
CMPLXtrig1(old,old);
CMPLXmult(tmp,old,new);
return(floatbailout());
}
TrigXTrigFractal()
{
#ifndef XFRACT
_LCMPLX ltmp2;
/* z = trig0(z)*trig1(z) */
LCMPLXtrig0(lold,ltmp);
LCMPLXtrig1(lold,ltmp2);
LCMPLXmult(ltmp,ltmp2,lnew);
if(overflow)
TryFloatFractal(TrigXTrigfpFractal);
return(longbailout());
#endif
}
/* call float version of fractal if integer math overflow */
TryFloatFractal(int (*fpFractal)())
{
overflow=0;
/* lold had better not be changed! */
old.x = lold.x; old.x /= fudge;
old.y = lold.y; old.y /= fudge;
tempsqrx = sqr(old.x);
tempsqry = sqr(old.y);
fpFractal();
lnew.x = new.x/fudge;
lnew.y = new.y/fudge;
return(0);
}
/********************************************************************/
/* Next six orbit functions are one type - extra functions are */
/* special cases written for speed. */
/********************************************************************/
TrigPlusSqrFractal() /* generalization of Scott and Skinner types */
{
#ifndef XFRACT
/* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
LCMPLXmult(lparm,ltmp,lnew); /* lnew = lparm*trig(lold) */
LCMPLXsqr_old(ltmp); /* ltmp = sqr(lold) */
LCMPLXmult(lparm2,ltmp,ltmp);/* ltmp = lparm2*sqr(lold) */
LCMPLXadd(lnew,ltmp,lnew); /* lnew = lparm*trig(lold)+lparm2*sqr(lold) */
return(longbailout());
#endif
}
TrigPlusSqrfpFractal() /* generalization of Scott and Skinner types */
{
/* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
CMPLXtrig0(old,tmp); /* tmp = trig(old) */
CMPLXmult(parm,tmp,new); /* new = parm*trig(old) */
CMPLXsqr_old(tmp); /* tmp = sqr(old) */
CMPLXmult(parm2,tmp,tmp2); /* tmp = parm2*sqr(old) */
CMPLXadd(new,tmp2,new); /* new = parm*trig(old)+parm2*sqr(old) */
return(floatbailout());
}
ScottTrigPlusSqrFractal()
{
#ifndef XFRACT
/* { z=pixel: z=trig(z)+sqr(z), |z|<BAILOUT } */
LCMPLXtrig0(lold,lnew); /* lnew = trig(lold) */
LCMPLXsqr_old(ltmp); /* lold = sqr(lold) */
LCMPLXadd(ltmp,lnew,lnew); /* lnew = trig(lold)+sqr(lold) */
return(longbailout());
#endif
}
ScottTrigPlusSqrfpFractal() /* float version */
{
/* { z=pixel: z=sin(z)+sqr(z), |z|<BAILOUT } */
CMPLXtrig0(old,new); /* new = trig(old) */
CMPLXsqr_old(tmp); /* tmp = sqr(old) */
CMPLXadd(new,tmp,new); /* new = trig(old)+sqr(old) */
return(floatbailout());
}
SkinnerTrigSubSqrFractal()
{
#ifndef XFRACT
/* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
LCMPLXtrig0(lold,lnew); /* lnew = trig(lold) */
LCMPLXsqr_old(ltmp); /* lold = sqr(lold) */
LCMPLXsub(lnew,ltmp,lnew); /* lnew = trig(lold)-sqr(lold) */
return(longbailout());
#endif
}
SkinnerTrigSubSqrfpFractal()
{
/* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
CMPLXtrig0(old,new); /* new = trig(old) */
CMPLXsqr_old(tmp); /* old = sqr(old) */
CMPLXsub(new,tmp,new); /* new = trig(old)-sqr(old) */
return(floatbailout());
}
TrigZsqrdfpFractal()
{
/* { z=pixel: z=trig(z*z), |z|<TEST } */
CMPLXsqr_old(tmp);
CMPLXtrig0(tmp,new);
return(floatbailout());
}
TrigZsqrdFractal() /* this doesn't work very well */
{
#ifndef XFRACT
/* { z=pixel: z=trig(z*z), |z|<TEST } */
LCMPLXsqr_old(ltmp);
LCMPLXtrig0(ltmp,lnew);
if(overflow)
TryFloatFractal(TrigZsqrdfpFractal);
return(longbailout());
#endif
}
SqrTrigFractal()
{
#ifndef XFRACT
/* { z=pixel: z=sqr(trig(z)), |z|<TEST} */
LCMPLXtrig0(lold,ltmp);
LCMPLXsqr(ltmp,lnew);
return(longbailout());
#endif
}
SqrTrigfpFractal()
{
/* SZSB(XYAXIS) { z=pixel, TEST=(p1+3): z=sin(z)*sin(z), |z|<TEST} */
CMPLXtrig0(old,tmp);
CMPLXsqr(tmp,new);
return(floatbailout());
}
Magnet1Fractal() /* Z = ((Z**2 + C - 1)/(2Z + C - 2))**2 */
{ /* In "Beauty of Fractals", code by Kev Allen. */
_CMPLX top, bot, tmp;
double div;
top.x = tempsqrx - tempsqry + floatparm->x - 1; /* top = Z**2+C-1 */
top.y = old.x * old.y;
top.y = top.y + top.y + floatparm->y;
bot.x = old.x + old.x + floatparm->x - 2; /* bot = 2*Z+C-2 */
bot.y = old.y + old.y + floatparm->y;
div = bot.x*bot.x + bot.y*bot.y; /* tmp = top/bot */
if (div < FLT_MIN) return(1);
tmp.x = (top.x*bot.x + top.y*bot.y)/div;
tmp.y = (top.y*bot.x - top.x*bot.y)/div;
new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y); /* Z = tmp**2 */
new.y = tmp.x * tmp.y;
new.y += new.y;
return(floatbailout());
}
Magnet2Fractal() /* Z = ((Z**3 + 3(C-1)Z + (C-1)(C-2) ) / */
/* (3Z**2 + 3(C-2)Z + (C-1)(C-2)+1) )**2 */
{ /* In "Beauty of Fractals", code by Kev Allen. */
_CMPLX top, bot, tmp;
double div;
top.x = old.x * (tempsqrx-tempsqry-tempsqry-tempsqry + T_Cm1.x)
- old.y * T_Cm1.y + T_Cm1Cm2.x;
top.y = old.y * (tempsqrx+tempsqrx+tempsqrx-tempsqry + T_Cm1.x)
+ old.x * T_Cm1.y + T_Cm1Cm2.y;
bot.x = tempsqrx - tempsqry;
bot.x = bot.x + bot.x + bot.x
+ old.x * T_Cm2.x - old.y * T_Cm2.y
+ T_Cm1Cm2.x + 1.0;
bot.y = old.x * old.y;
bot.y += bot.y;
bot.y = bot.y + bot.y + bot.y
+ old.x * T_Cm2.y + old.y * T_Cm2.x
+ T_Cm1Cm2.y;
div = bot.x*bot.x + bot.y*bot.y; /* tmp = top/bot */
if (div < FLT_MIN) return(1);
tmp.x = (top.x*bot.x + top.y*bot.y)/div;
tmp.y = (top.y*bot.x - top.x*bot.y)/div;
new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y); /* Z = tmp**2 */
new.y = tmp.x * tmp.y;
new.y += new.y;
return(floatbailout());
}
LambdaTrigFractal()
{
#ifndef XFRACT
LONGXYTRIGBAILOUT();
LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
lold = lnew;
return(0);
#endif
}
LambdaTrigfpFractal()
{
FLOATXYTRIGBAILOUT();
CMPLXtrig0(old,tmp); /* tmp = trig(old) */
CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
old = new;
return(0);
}
/* bailouts are different for different trig functions */
LambdaTrigFractal1()
{
#ifndef XFRACT
LONGTRIGBAILOUT(); /* sin,cos */
LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
lold = lnew;
return(0);
#endif
}
LambdaTrigfpFractal1()
{
FLOATTRIGBAILOUT(); /* sin,cos */
CMPLXtrig0(old,tmp); /* tmp = trig(old) */
CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
old = new;
return(0);
}
LambdaTrigFractal2()
{
#ifndef XFRACT
LONGHTRIGBAILOUT(); /* sinh,cosh */
LCMPLXtrig0(lold,ltmp); /* ltmp = trig(lold) */
LCMPLXmult(*longparm,ltmp,lnew); /* lnew = longparm*trig(lold) */
lold = lnew;
return(0);
#endif
}
LambdaTrigfpFractal2()
{
#ifndef XFRACT
FLOATHTRIGBAILOUT(); /* sinh,cosh */
CMPLXtrig0(old,tmp); /* tmp = trig(old) */
CMPLXmult(*floatparm,tmp,new); /* new = longparm*trig(old) */
old = new;
return(0);
#endif
}
ManOWarFractal()
{
#ifndef XFRACT
/* From Art Matrix via Lee Skinner */
lnew.x = ltempsqrx - ltempsqry + ltmp.x + longparm->x;
lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y + longparm->y;
ltmp = lold;
return(longbailout());
#endif
}
ManOWarfpFractal()
{
/* From Art Matrix via Lee Skinner */
/* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
new.x = tempsqrx - tempsqry + tmp.x + floatparm->x;
new.y = 2.0 * old.x * old.y + tmp.y + floatparm->y;
tmp = old;
return(floatbailout());
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -