📄 fractals.c
字号:
for(i=0;i<degree;i++)
/* color in alternating shades with iteration according to
which root of 1 it converged to */
if(distance(roots[i],old) < threshold)
{
if (basin==2) {
tmpcolor = 1+(i&7)+((color&1)<<3);
} else {
tmpcolor = 1+i;
}
break;
}
if(tmpcolor == -1)
color = maxcolor;
else
color = tmpcolor;
}
return(1);
}
new.x = d1overd * new.x + roverd;
new.y *= d1overd;
/* Watch for divide underflow */
if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN)
return(1);
else
{
t2 = 1.0 / t2;
old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
}
return(0);
}
#endif /* newton code only used by xfractint */
complex_mult(arg1,arg2,pz)
_CMPLX arg1,arg2,*pz;
{
pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
return(0);
}
complex_div(numerator,denominator,pout)
_CMPLX numerator,denominator,*pout;
{
double mod;
if((mod = modulus(denominator)) < FLT_MIN)
return(1);
conjugate(&denominator);
complex_mult(numerator,denominator,pout);
pout->x = pout->x/mod;
pout->y = pout->y/mod;
return(0);
}
#ifndef XFRACT
struct MP mproverd, mpd1overd, mpthreshold,sqrmpthreshold;
struct MP mpt2;
struct MP mpone;
extern struct MPC MPCone;
extern int MPOverflow;
#endif
int MPCNewtonFractal()
{
#ifndef XFRACT
MPOverflow = 0;
mpctmp = MPCpow(mpcold,degree-1);
mpcnew.x = *pMPsub(*pMPmul(mpctmp.x,mpcold.x),*pMPmul(mpctmp.y,mpcold.y));
mpcnew.y = *pMPadd(*pMPmul(mpctmp.x,mpcold.y),*pMPmul(mpctmp.y,mpcold.x));
mpctmp1.x = *pMPsub(mpcnew.x, MPCone.x);
mpctmp1.y = *pMPsub(mpcnew.y, MPCone.y);
if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
{
if(fractype==MPNEWTBASIN)
{
int tmpcolor;
int i;
tmpcolor = -1;
for(i=0;i<degree;i++)
if(pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
{
if(basin==2)
tmpcolor = 1+(i&7) + ((color&1)<<3);
else
tmpcolor = 1+i;
break;
}
if(tmpcolor == -1)
color = maxcolor;
else
color = tmpcolor;
}
return(1);
}
mpcnew.x = *pMPadd(*pMPmul(mpd1overd,mpcnew.x),mproverd);
mpcnew.y = *pMPmul(mpcnew.y,mpd1overd);
mpt2 = MPCmod(mpctmp);
mpt2 = *pMPdiv(mpone,mpt2);
mpcold.x = *pMPmul(mpt2,(*pMPadd(*pMPmul(mpcnew.x,mpctmp.x),*pMPmul(mpcnew.y,mpctmp.y))));
mpcold.y = *pMPmul(mpt2,(*pMPsub(*pMPmul(mpcnew.y,mpctmp.x),*pMPmul(mpcnew.x,mpctmp.y))));
new.x = *pMP2d(mpcold.x);
new.y = *pMP2d(mpcold.y);
return(MPOverflow);
#endif
}
Barnsley1Fractal()
{
#ifndef XFRACT
/* Barnsley's Mandelbrot type M1 from "Fractals
Everywhere" by Michael Barnsley, p. 322 */
/* calculate intermediate products */
oldxinitx = multiply(lold.x, longparm->x, bitshift);
oldyinity = multiply(lold.y, longparm->y, bitshift);
oldxinity = multiply(lold.x, longparm->y, bitshift);
oldyinitx = multiply(lold.y, longparm->x, bitshift);
/* orbit calculation */
if(lold.x >= 0)
{
lnew.x = (oldxinitx - longparm->x - oldyinity);
lnew.y = (oldyinitx - longparm->y + oldxinity);
}
else
{
lnew.x = (oldxinitx + longparm->x - oldyinity);
lnew.y = (oldyinitx + longparm->y + oldxinity);
}
return(longbailout());
#endif
}
Barnsley1FPFractal()
{
/* Barnsley's Mandelbrot type M1 from "Fractals
Everywhere" by Michael Barnsley, p. 322 */
/* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
/* calculate intermediate products */
foldxinitx = old.x * floatparm->x;
foldyinity = old.y * floatparm->y;
foldxinity = old.x * floatparm->y;
foldyinitx = old.y * floatparm->x;
/* orbit calculation */
if(old.x >= 0)
{
new.x = (foldxinitx - floatparm->x - foldyinity);
new.y = (foldyinitx - floatparm->y + foldxinity);
}
else
{
new.x = (foldxinitx + floatparm->x - foldyinity);
new.y = (foldyinitx + floatparm->y + foldxinity);
}
return(floatbailout());
}
Barnsley2Fractal()
{
#ifndef XFRACT
/* An unnamed Mandelbrot/Julia function from "Fractals
Everywhere" by Michael Barnsley, p. 331, example 4.2 */
/* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
/* calculate intermediate products */
oldxinitx = multiply(lold.x, longparm->x, bitshift);
oldyinity = multiply(lold.y, longparm->y, bitshift);
oldxinity = multiply(lold.x, longparm->y, bitshift);
oldyinitx = multiply(lold.y, longparm->x, bitshift);
/* orbit calculation */
if(oldxinity + oldyinitx >= 0)
{
lnew.x = oldxinitx - longparm->x - oldyinity;
lnew.y = oldyinitx - longparm->y + oldxinity;
}
else
{
lnew.x = oldxinitx + longparm->x - oldyinity;
lnew.y = oldyinitx + longparm->y + oldxinity;
}
return(longbailout());
#endif
}
Barnsley2FPFractal()
{
/* An unnamed Mandelbrot/Julia function from "Fractals
Everywhere" by Michael Barnsley, p. 331, example 4.2 */
/* calculate intermediate products */
foldxinitx = old.x * floatparm->x;
foldyinity = old.y * floatparm->y;
foldxinity = old.x * floatparm->y;
foldyinitx = old.y * floatparm->x;
/* orbit calculation */
if(foldxinity + foldyinitx >= 0)
{
new.x = foldxinitx - floatparm->x - foldyinity;
new.y = foldyinitx - floatparm->y + foldxinity;
}
else
{
new.x = foldxinitx + floatparm->x - foldyinity;
new.y = foldyinitx + floatparm->y + foldxinity;
}
return(floatbailout());
}
JuliaFractal()
{
#ifndef XFRACT
/* used for C prototype of fast integer math routines for classic
Mandelbrot and Julia */
lnew.x = ltempsqrx - ltempsqry + longparm->x;
lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
return(longbailout());
#elif !defined(__386BSD__)
fprintf(stderr,"JuliaFractal called\n");
exit(-1);
#endif
}
JuliafpFractal()
{
/* floating point version of classical Mandelbrot/Julia */
/* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
new.x = tempsqrx - tempsqry + floatparm->x;
new.y = 2.0 * old.x * old.y + floatparm->y;
return(floatbailout());
}
static double f(double x, double y)
{
return(x - x*y);
}
static double g(double x, double y)
{
return(-y + x*y);
}
LambdaFPFractal()
{
/* variation of classical Mandelbrot/Julia */
/* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
tempsqrx = old.x - tempsqrx + tempsqry;
tempsqry = -(old.y * old.x);
tempsqry += tempsqry + old.y;
new.x = floatparm->x * tempsqrx - floatparm->y * tempsqry;
new.y = floatparm->x * tempsqry + floatparm->y * tempsqrx;
return(floatbailout());
}
LambdaFractal()
{
#ifndef XFRACT
/* variation of classical Mandelbrot/Julia */
/* in complex math) temp = Z * (1-Z) */
ltempsqrx = lold.x - ltempsqrx + ltempsqry;
ltempsqry = lold.y
- multiply(lold.y, lold.x, bitshiftless1);
/* (in complex math) Z = Lambda * Z */
lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
- multiply(longparm->y, ltempsqry, bitshift);
lnew.y = multiply(longparm->x, ltempsqry, bitshift)
+ multiply(longparm->y, ltempsqrx, bitshift);
return(longbailout());
#endif
}
SierpinskiFractal()
{
#ifndef XFRACT
/* following code translated from basic - see "Fractals
Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
lnew.x = (lold.x << 1); /* new.x = 2 * old.x */
lnew.y = (lold.y << 1); /* new.y = 2 * old.y */
if(lold.y > ltmp.y) /* if old.y > .5 */
lnew.y = lnew.y - ltmp.x; /* new.y = 2 * old.y - 1 */
else if(lold.x > ltmp.y) /* if old.x > .5 */
lnew.x = lnew.x - ltmp.x; /* new.x = 2 * old.x - 1 */
/* end barnsley code */
return(longbailout());
#endif
}
SierpinskiFPFractal()
{
/* following code translated from basic - see "Fractals
Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
new.x = old.x + old.x;
new.y = old.y + old.y;
if(old.y > .5)
new.y = new.y - 1;
else if (old.x > .5)
new.x = new.x - 1;
/* end barnsley code */
return(floatbailout());
}
LambdaexponentFractal()
{
/* found this in "Science of Fractal Images" */
FLOATEXPBAILOUT();
FPUsincos (&old.y,&siny,&cosy);
if (old.x >= rqlim && cosy >= 0.0) return(1);
tmpexp = exp(old.x);
tmp.x = tmpexp*cosy;
tmp.y = tmpexp*siny;
/*multiply by lamda */
new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
old = new;
return(0);
}
LongLambdaexponentFractal()
{
#ifndef XFRACT
/* found this in "Science of Fractal Images" */
LONGEXPBAILOUT();
SinCos086 (lold.y, &lsiny, &lcosy);
if (lold.x >= llimit && lcosy >= 0L) return(1);
longtmp = Exp086(lold.x);
ltmp.x = multiply(longtmp, lcosy, bitshift);
ltmp.y = multiply(longtmp, lsiny, bitshift);
lnew.x = multiply(longparm->x, ltmp.x, bitshift)
- multiply(longparm->y, ltmp.y, bitshift);
lnew.y = multiply(longparm->x, ltmp.y, bitshift)
+ multiply(longparm->y, ltmp.x, bitshift);
lold = lnew;
return(0);
#endif
}
FloatTrigPlusExponentFractal()
{
/* another Scientific American biomorph type */
/* z(n+1) = e**z(n) + trig(z(n)) + C */
if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
tmpexp = exp(old.x);
FPUsincos (&old.y,&siny,&cosy);
CMPLXtrig0(old,new);
/*new = trig(old) + e**old + C */
new.x += tmpexp*cosy + floatparm->x;
new.y += tmpexp*siny + floatparm->y;
return(floatbailout());
}
LongTrigPlusExponentFractal()
{
#ifndef XFRACT
/* calculate exp(z) */
/* domain check for fast transcendental functions */
TRIG16CHECK(lold.x);
TRIG16CHECK(lold.y);
longtmp = Exp086(lold.x);
SinCos086 (lold.y, &lsiny, &lcosy);
LCMPLXtrig0(lold,lnew);
lnew.x += multiply(longtmp, lcosy, bitshift) + longparm->x;
lnew.y += multiply(longtmp, lsiny, bitshift) + longparm->y;
return(longbailout());
#endif
}
MarksLambdaFractal()
{
/* Mark Peterson's variation of "lambda" function */
/* Z1 = (C^(exp-1) * Z**2) + C */
ltmp.x = ltempsqrx - ltempsqry;
ltmp.y = multiply(lold.x ,lold.y ,bitshiftless1);
lnew.x = multiply(lcoefficient.x, ltmp.x, bitshift)
- multiply(lcoefficient.y, ltmp.y, bitshift) + longparm->x;
lnew.y = multiply(lcoefficient.x, ltmp.y, bitshift)
+ multiply(lcoefficient.y, ltmp.x, bitshift) + longparm->y;
return(longbailout());
}
MarksLambdafpFractal()
{
/* Mark Peterson's variation of "lambda" function */
/* Z1 = (C^(exp-1) * Z**2) + C */
tmp.x = tempsqrx - tempsqry;
tmp.y = old.x * old.y *2;
new.x = coefficient.x * tmp.x - coefficient.y * tmp.y + floatparm->x;
new.y = coefficient.x * tmp.y + coefficient.y * tmp.x + floatparm->y;
return(floatbailout());
}
long XXOne, FgOne, FgTwo;
UnityFractal()
{
#ifndef XFRACT
/* brought to you by Mark Peterson - you won't find this in any fractal
books unless they saw it here first - Mark invented it! */
XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delmin))
return(1);
lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
lnew=lold; /* TW added this line */
return(0);
#endif
}
#define XXOne new.x
UnityfpFractal()
{
/* brought to you by Mark Peterson - you won't find this in any fractal
books unless they saw it here first - Mark invented it! */
XXOne = sqr(old.x) + sqr(old.y);
if((XXOne > 2.0) || (fabs(XXOne - 1.0) < ddelmin))
return(1);
old.y = (2.0 - XXOne)* old.x;
old.x = (2.0 - XXOne)* old.y;
new=old; /* TW added this line */
return(0);
}
#undef XXOne
Mandel4Fractal()
{
/* By writing this code, Bert has left behind the excuse "don't
know what a fractal is, just know how to make'em go fast".
Bert is hereby declared a bonafide fractal expert! Supposedly
this routine calculates the Mandelbrot/Julia set based on the
polynomial z**4 + lambda, but I wouldn't know -- can't follow
all that integer math speedup stuff - Tim */
/* first, compute (x + iy)**2 */
lnew.x = ltempsqrx - ltempsqry;
lnew.y = multiply(lold.x, lold.y, bitshiftless1);
if (longbailout()) return(1);
/* then, compute ((x + iy)**2)**2 + lambda */
lnew.x = ltempsqrx - ltempsqry + longparm->x;
lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
return(longbailout());
}
Mandel4fpFractal()
{
/* first, compute (x + iy)**2 */
new.x = tempsqrx - tempsqry;
new.y = old.x*old.y*2;
if (floatbailout()) return(1);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -