/* mandelbrot tile                                                                  
mt -zone y0 x0 y1 x1 -size ny nx -ofile file ...                                    
minimum memory varsion. Allocate only one line                                      
*/                                                                                  
                                                                                    
#include <stdio.h>                                                                  
#include <math.h>                                                                   
                                                                                    
#define ERROR   -1                                                                  
                                                                                    
double  zw[4]   = {-2,-2,2,2};                                                      
int     size[2] = {24,80};                                                          
int     kf      = 0;            /* default to mandelbrot */                         
double  eps     =  1.0e-10;     /* -eps epsilon */                                  
double  lambda  = 0.35;         /* used by z <- lambda*cos(z) + c etc. */           
double  a       = -1.0;         /* a,c used by cubic and ecc */                     
double  c       = 0.0;                                                              
int     bailout = 35;           /* maximum iteration */                             
char   *ofile   = NULL;                                                             
char   *lut     = ".123456789ABCDEF";                                               
int     nlut;                                                                       
int     fcmax   = 8;                                                                
                                                                                    
int j_orbit(), j3_orbit(), j4_orbit(), j5_orbit(), j6_orbit(), j_wis(),             
    j_bes(), j_cubic(), j_ecc();                                                    
                                                                                    
int (*func[])() = {                                                                 
   j_orbit, j3_orbit, j4_orbit, j5_orbit, j6_orbit, j_wis, j_bes,                   
   j_cubic, j_ecc                                                                   
   };                                                                               
                                                                                    
main(argc, argv)                                                                    
int argc;                                                                           
char **argv;                                                                        
{                                                                                   
extern  double  atof();                                                             
extern  char   *malloc();                                                           
FILE   *ofp = NULL;                                                                 
char   *s;                                                                          
char   *rup, *t;                                                                    
int     i, j, k;                                                                    
double  x, y, dx, dy;                                                               
                                                                                    
for(argv++, argc--; 0 < argc--; argv++) {                                           
    s = *argv;                                                                      
    if (s[0] == '-') switch(s[1]) {                                                 
        case 'a': a = atof(*++argv);                                                
                  argc--;                                                           
        break;                                                                      
        case 'c': c = atof(*++argv);                                                
                  argc--;                                                           
        break;                                                                      
        case 'z':   /* zone */                                                      
                zw[0] = atof(*++argv);                                              
                zw[1] = atof(*++argv);                                              
                zw[2] = atof(*++argv);                                              
                zw[3] = atof(*++argv);                                              
                argc -= 4;                                                          
        break;                                                                      
        case 's':   /* size */                                                      
                size[0] = atoi(*++argv);                                            
                size[1] = atoi(*++argv);                                            
                argc -=2;                                                           
        break;                                                                      
        case 'l':       /* lut */                                                   
                if (s[2] == 'a')                                                    
                        lambda = atof(*++argv);                                     
                else    lut = *++argv;                                              
                argc--;                                                             
        break;                                                                      
        case 'o':       /* ofile */                                                 
                ofile = *++argv;                                                    
                argc--;                                                             
        break;                                                                      
        case 'b':       /* -bailout */                                              
                bailout = atoi(*++argv);                                            
                argc--;                                                             
        break;                                                                      
        case 'e':       /* -eps */                                                  
                eps = atof(*++argv);                                                
                argc--;                                                             
        break;                                                                      
        case 'f':       /* -function */                                             
                kf = atoi(*++argv);                                                 
                if (kf > fcmax || kf < 0 ) kf = 0;                                  
                argc--;                                                             
        break;                                                                      
        default: break;                                                             
        }                                                                           
    }                                                                               
nlut = strlen(lut);                                                                 
                                                                                    
if (ofile) {                                                                        
    ofp = fopen(ofile, "w");                                                        
    if (ofp == NULL) return ERROR;                                                  
    }                                                                               
else    ofp = stdout;                                                               
                                                                                    
fprintf(ofp,                                                                        
    "-fc %d -size %d %d -zone %f %f %f %f -lambda %f -eps %f -bailout %d -lut %s\n",
    kf, size[0], size[1], zw[0], zw[1], zw[2], zw[3], lambda, eps, bailout, lut);   
                                                                                    
if (NULL == (rup = malloc(1+size[1]))) return ERROR;                                
dy=(zw[2]-zw[0])/(size[0]-1);                                                       
dx=(zw[3]-zw[1])/(size[1]-1);                                                       
for (i=0, y=zw[0]; i < size[0]; i++, y +=dy) {                                      
    for (t=rup, j=0, x=zw[1]; j < size[1]; j++, x += dx) {                          
        k = (func[kf])(y,x);                                                        
       *t++ = lut[k % nlut];                                                        
        }                                                                           
    for (t = rup, j=0; j < size[1]; j++) fputc(*t++, ofp);                          
    fputc('\n', ofp);                                                               
    }                                                                               
if (ofile) fclose (ofp);                                                            
return 0;                                                                           
}                                                                                   
                                                                                    
/*      dynamic system z<-z^2+c with global parameters eps, bailout */              
                                                                                    
int     j_orbit(cy, cx)                                                             
double  cy, cx;                                                                     
{                                                                                   
double  x = 0.0, y = 0.0;                                                           
double  x1, y1, x2, y2;                                                             
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    x2= x*x; y2=y*y;                                                                
    if (4.1 < x2+y2) return j;                                                      
    x1 = cx + x2-y2;                                                                
    y1 = cy + 2*x*y;                                                                
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*      dynamic system z<-z^3+c  */                                                 
                                                                                    
int     j3_orbit(cy, cx)                                                            
double  cy, cx;                                                                     
{                                                                                   
double  x = 0.0, y = 0.0;                                                           
double  x1, y1, x2, y2;                                                             
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    x2= x*x; y2=y*y;                                                                
    if (4.1 < x2+y2) return j;                                                      
    x1 = cx + x*(x2-3*y2);                                                          
    y1 = cy + y*(3*x2-y2);                                                          
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*      dynamic system z<-z^4+c  */                                                 
                                                                                    
int     j4_orbit(cy, cx)                                                            
double  cy, cx;                                                                     
{                                                                                   
double  x = 0.0, y = 0.0;                                                           
double  x1, y1, x2, y2;                                                             
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    x2= x*x; y2=y*y;                                                                
    if (4.1 < x2+y2) return j;                                                      
    x1 = cx + x2*x2+y2*y2-6*x2*y2;                                                  
    y1 = cy + 4*x*y*(x2-y2);                                                        
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*      dynamic system z<-z^5+c  */                                                 
                                                                                    
int     j5_orbit(cy, cx)                                                            
double  cy, cx;                                                                     
{                                                                                   
double  x = 0.0, y = 0.0;                                                           
double  x1, y1, x2, y2, x2y2, x4, y4;                                               
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    x2= x*x; y2=y*y; x2y2=x2*y2; x4=x2*x2; y4=y2*y2;                                
    if (4.1 < x2+y2) return j;                                                      
    x1 = cx + x*(x2*x2+5*y2*y2-10*x2y2);                                            
    y1 = cy + y*(5*x4-10*x2y2+y4);                                                  
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*      dynamic system z<-z^6+c  */                                                 
                                                                                    
int     j6_orbit(cy, cx)                                                            
double  cy, cx;                                                                     
{                                                                                   
double  x = 0.0, y = 0.0;                                                           
double  x1, y1, x2, y2, x2y2, x4, y4;                                               
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    x2= x*x; y2=y*y; x2y2=x2*y2; x4=x2*x2; y4=y2*y2;                                
    if (4.1 < x2+y2) return j;                                                      
    x1 = cx + x2*x2*(x2-16*y2)+y2*y2*(15*x2-y2);                                    
    y1 = cy + x*y*(6*(x4+y4)-20*x2y2);                                              
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*      j_wis newton raphson for w = lambda*z+exp(z)            */                  
/*      iterate z <- z-f(z)/f'(z) until |f(z)/f'(z)|<eps        */                  
/*      f(z) = lambda*z+exp(z)-w                                */                  
/*      f'(z) = lambda+exp(z)                                   */                  
/*      here we set z = x+iy so exp(z)=exp(x)*cos(y)+i sin(y)   */                  
/*      here f(0)=1-w and f'(0)=1+lambda                        */                  
/*      so start the iteration with cy, cx                      */                  
                                                                                    
                                                                                    
int     j_wis(cy, cx)                                                               
double  cy, cx;                                                                     
{                                                                                   
double  exp(), cos(), sin();                                                        
double  x = cx, y = cy;                                                             
double  a, b, c, d, delta, ex, ezr, ezi, dx, dy;                                    
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    if (10000 < x*x+y*y) return 0;      /* infinity */                              
    ex = exp(x);  ezr = ex*cos(y); ezi = ex*sin(y);                                 
    delta = lambda*lambda+ex*ex+2*lambda*ezr;                                       
    a = lambda*x-cx+ezr;                                                            
    b = ezi-cy;                                                                     
    c = lambda + ezr;                                                               
/*  d = ezi          */                                                             
    dx = (a*c + b*ezi)/delta;                                                       
    dy = (b*c - a*ezi)/delta;                                                       
    if ((dx*dx+dy*dy) < eps*eps) return j;                                          
    x = x - dx;                                                                     
    y = y - dy;                                                                     
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/* solve w=j(z) by z[n+1]=z[n]-(j(z[n])-w)/j'(z[n])     */                          
/* compute j(z), j'(z) by series                        */                          
/* j(z) = sum {(z^2n)/4^n*(n!)^2 : n=0, infinity}       */                          
/* j(lambda, z) =                                       */                          
/* sum {z^2n/ 4^n*n! * product(lambda+i:i=1..n): n=0..infinity */                   
                                                                                    
int     j_bes(cy, cx)                                                               
double  cy, cx;                                                                     
{                                                                                   
int     j_bbox();                                                                   
double  x = cx, y = cy;                                                             
double  a,b,c,d, delta, dx, dy;                                                     
int     j, k;                                                                       
for(j = 0; j < bailout; j++) {                                                      
    if (10000 < x*x+y*y) return 0;              /* infinity */                      
    k = j_bbox(&a, &b, &c, &d, x, y, eps);      /* compute j(z), j'(z) */           
    a = a - cx;                                                                     
    b = b - cy;                                                                     
    delta = c*c + d*d;                                                              
    dx = (a*c + b*d)/delta;                                                         
    dy = (b*c - a*d)/delta;                                                         
    if ((dx*dx+dy*dy) < eps*eps) return j;                                          
    x = x - dx;                                                                     
    y = y - dy;                                                                     
    }                                                                               
return  0;                                                                          
}                                                                                   
                                                                                    
int     j_bbox(jx, jy, djx, djy, x, y, eps)                                         
double *jx, *jy, *djx, *djy;                                                        
double  x, y, eps;                                                                  
{                                                                                   
double  tx, ty, ux, uy, delta;                                                      
int     i;                                                                          
*jx = tx = 1.0; *jy = *djx = *djy = ty = 0.0;                                       
for (i=1; i<bailout; i++) {                                                         
    delta = 4*i*(i+lambda);                                                         
    ux = (tx*x-ty*y)/delta;                                                         
    uy = (tx*y+ty*x)/delta;                                                         
    *jx += (tx = (ux*x-uy*y));                                                      
    *jy += (ty = (ux*y+uy*x));                                                      
    *djx += 2*i*ux;                                                                 
    *djy += 2*i*uy;                                                                 
    if (ux*ux+uy*uy < eps*eps) break;                                               
    }                                                                               
return i;                                                                           
}                                                                                   
                                                                                    
                                                                                    
/*  solve z^3=az+b                                                                  
*   F(z) = 2*z^3+b/(3*x^2-a)                                                        
*   iterate z <- z-f(z)/f'(z) until |f(z)/f'(z)|<eps                                
*   here we set z = x+iy                                                            
*   start the iteration with cy, cx                                                 
*/                                                                                  
                                                                                    
int     j_cubic(cy, cx)                                                             
double  cy, cx;                                                                     
{                                                                                   
double  x = cx, y = cy;                                                             
double  x1, x2, y1, y2, nx, ny, dx, dy, delta;                                      
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    if (1e20 < x*x+y*y) return 0;      /* infinity */                               
    x2 = x*x; y2 = y*y;                                                             
    dx = 3*(x2-y2) - a;                                                             
    dy = 6*x*y;                                                                     
    delta =dx*dx+dy*dy;                                                             
    if (delta < eps) return 0;                                                      
    nx = c + 2*x*(x2-3*y2);                                                         
    ny = 2*y*(3*x2-y2);                                                             
    x1 = (nx*dx+ny*dy)/delta;                                                       
    y1 = (ny*dx-nx*dy)/delta;                                                       
    dx = x-x1;                                                                      
    dy = y-y1;                                                                      
    delta = dx*dx + dy*dy;                                                          
    if (delta < eps) return j;                                                      
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}                                                                                   
                                                                                    
/*  solve  Y^2 = X^3+aX+b                                                           
*   F(z) = 2*z^3-b+Y^Y/(3*z^2+a)                                                    
*   iterate z <- z-f(z)/f'(z) until |f(z)/f'(z)|<eps                                
*   here we set z = x+iy                                                            
*   start the iteration with cy, cx                                                 
*/                                                                                  
                                                                                    
int     j_ecc(cy, cx)                                                               
double  cy, cx;                                                                     
{                                                                                   
double  x = cx, y = cy;                                                             
double  reysq = cx*cx - cy*cy;                                                      
double  imysq = 2*cx*cy;                                                            
double  x1, x2, y1, y2, nx, ny, dx, dy, delta;                                      
int     j;                                                                          
for(j = 0; j < bailout; j++) {                                                      
    if (1e10 < x*x+y*y) return 0;      /* infinity */                               
    x2 = x*x; y2 = y*y;                                                             
    dx = 3*(x2-y2) + a;                                                             
    dy = 6*x*y;                                                                     
    delta =dx*dx+dy*dy;                                                             
    if (delta < eps) return 0;                                                      
    nx = c + 2*x*(x2-3*y2) - reysq;                                                 
    ny = 2*y*(3*x2-y2) - imysq;                                                     
    x1 = (nx*dx+ny*dy)/delta;                                                       
    y1 = (ny*dx-nx*dy)/delta;                                                       
    dx = x-x1;                                                                      
    dy = y-y1;                                                                      
    delta = dx*dx + dy*dy;                                                          
    if (delta < eps) return j;                                                      
    x = x1;                                                                         
    y = y1;                                                                         
    }                                                                               
return 0;                                                                           
}