/*  This file is used for general iteration of rational
    functions.  It consists of the following routines:

    pick_map(name)
    pick_map(name,param1,...) - 
        Picks the map to be iterated.
	Name is a string specifying the map.
	Parameters are optional for some values of name.

    arg_map(argv,argc,i)
	scans argument list for map

    iter(z)
	iterates selected map on z;
	sets z=INFINITY for large values

    attractor(at,na)   - returns a list at[] of
	the attractor for each critical point of f.

    critical_points(cps,nc) - returns a list of
	the critical points of f.

Possible names:
With no parameters:
    P3_4 - deg 4 map, Pilgrim example I
    S3_4 - deg 4 map with S3 symmetry
    S3Z2_4 - Z/2 quotient of deg 4 map
    S3Z3_4 - Z/3 quotient of deg 4 map
    S3S3_4 - S3  quotient of deg 4 map
    A5_11 - deg 11 map with A5 symmetry
    A5_19 - deg 19 map with A5 symmetry
    A5_29 - deg 29 map with A5 symmetry
    V4_15 - deg 15 map with V4 symmetry
    NEW_3 - Newton's method failing for cubics
    PET_4  - Preparabolic petals

With one parameter:
    BAS_2  - deg 2 rational map, w/cp of per 2.
    POLY_2 - quadratic polynomial z^2+c
    Z3_3   - cubic polynomial z^3+c
    Z4_4   - quartic polynomial z^4+c
    FCP_3  - cubic polynomial z^3 + c z^2
    P3_3   - deg 3 map, Pilgrim example II 
    Q3_3   - deg 3 map, variant.
    NCUB_3 - Newton's method for general cubics 
    CAN_5  - deg 5 map, with CantorxS^1 Julia set
    EXP    - a exp z
    SIN    - a sin z
    ROT_3  - a z^2 (z-3)/(1-3z)
    NCUB_3 - Newton's method for cubics
    NAP_2  - Napolean mapping
    FIX_3  - cubic rational map with 3 fixed critical points

With two parameters:
    FP_3   - cubic polynomial a(z^3/3-(1/2)(t+1/t)z^2+z)
    PAR_2  - rational map with fixed point of mult. b.
    ESC_2  - exp(a) (z + 1/z + b)
    BIC_2  - bicritical b/(z+ 1/z + a) (HALF_2)
    BIC_3  - bicritical cubic
    BIC_4  - bicritical quartic

With six parameters:
    RAT_2 - arbitrary deg 2 rational map 
		(az^2+bz + c)/(dz^2+ez + f)
*/

/*  To add a new function:
	(1) Define fn, fncps, fnnc.
	(2) Modify pick_map.
	(3) Modify arg_map if fn takes params.
	(4) Add name to fnusage and list above.
*/

/* Global variables  - map, critical points, number
   of critical points, parameters on which the
   map may depend
*/

#include "cx.h"
#include <stdio.h>
#define INFINITY 1e+20

complex (*f)();
complex *cps;
int nc;
complex param1, param2, param3, param4, param5, param6;

/* Defined variables - point to be considered
   infinity; resolution of different attractors;
   number of iterates to seek an attractor.
*/

#define ATT_EPS 0.0001
#define ATT_ITER 10000
static double att_iter = ATT_ITER;

/* MAPPINGS AND THEIR CRITICAL POINTS */

/* z^2 + a */
complex poly_2(z)
	complex z;
{
	return(add(mult(z,z),param1));
}
static complex poly_2cps[2]={0.0,0.0,INFINITY,0.0};
static int poly_2nc = 2;

/* z^3 + a */
complex z3_3(z)
	complex z;
{
	return(add(mult(mult(z,z),z),param1));
}
static complex z3_3cps[2]={0.0,0.0,INFINITY,0.0};
static int z3_3nc = 2;

/* z^4 + a */
complex z4_4(z)
	complex z;
{
	return(add(mult(mult(z,z),mult(z,z)),param1));
}
static complex z4_4cps[2]={0.0,0.0,INFINITY,0.0};
static int z4_4nc = 2;

/* z^3 + a z^2 */
complex fcp_3(z)
	complex z;
{
	complex z2, z3;
	z2 = mult(z,z);
	z3 = mult(z2,z);
	return(add(z3,mult(param1,z2)));
}
static complex fcp_3cps[3]=
	{0.0,0.0,0.0,0.0,INFINITY,0.0};
static int fcp_3nc = 3;
fcp_3set()
{
	static complex m2over3 = {-0.66666666666666,0.0};
	fcp_3cps[0] = mult(m2over3,param1);
}

/* a(z^3/3-(1/2)(t+1/t)z^2 + z) */
/* Critical points at t, 1/t; multiplier at fp = a. */
/* t = param1, a = param2 */
complex fp_3(z)
	complex z;
{
	complex b, bz2, w, z2, z3;
	z2 = mult(z,z);
	z3 = mult(z2,z);
	b = add(param1,recip(param1));
	bz2 = mult(b,z2);
	w.x = z3.x/3 - bz2.x/2 + z.x; 
	w.y = z3.y/3 - bz2.y/2 + z.y; 
	w = mult(param2,w);
	return(w);
}
complex fp_3cps[3];
static int fp_3nc = 3;
fp_3set()
{
	fp_3cps[0]   = param1;
	fp_3cps[1]   = recip(param1);
	fp_3cps[2].x = INFINITY;
	fp_3cps[2].y = 0.0;
}

/*  (az^2+bz + c)/(dz^2+ez + f) */
complex rat_2(z)
	complex z;
{
	static complex one={1,0};
	static complex infinity={INFINITY,0};
	complex num, den, z2;
	z2 = mult(z,z);
	num = add(add(mult(z2,param1),mult(z,param2)),param3);
	den = add(add(mult(z2,param4),mult(z,param5)),param6);
	if(cx_abs(den) == 0.0) return(infinity);
	return(divide(num,den));
}

/* CPS: (ae-bd)z^2 + 2(af-cd)z + bf-ce = 0 */
complex rat_2cps[2];
static int rat_2nc = 2;
rat_2set()
{
	complex a,b,c,disc,e;
	static complex two = {2.0,0.0}, four = {4.0,0.0};
	a = sub(mult(param1,param5),mult(param2,param4)); 
	b = mult(two,sub(mult(param1,param6),mult(param3,param4)));
	c = sub(mult(param2,param6),mult(param3,param5));
	disc = sub(mult(b,b),mult(four,mult(a,c)));
	e = cx_sqrt(disc);
	b.x = -b.x; b.y = -b.y;
	if(a.x == 0 && a.y == 0)
	{	rat_2cps[0] = divide(c,b);
		rat_2cps[1].x = INFINITY;
		rat_2cps[1].y = 0;
		return;
	}
	rat_2cps[0] = divide(add(b,e),mult(two,a));
	rat_2cps[1] = divide(sub(b,e),mult(two,a));
}

/* (a*(z-1)^3/(1 - 3*z) */
/* a=-1 a semi-basilica */

complex q3_3(z)
	complex z;
{
	static complex one={1.0,0.0},
		three={3.0,0.0};
	complex a, numer, denom;

	a = sub(z,one);
	a = mult(a,mult(a,a));
	numer= mult(a,param1);
	denom= sub(one,mult(three,z));
	return(divide(numer,denom));
}
static complex q3_3cps[2]={0,0,INFINITY,0};
static int q3_3nc = 2;

/* (a*(-1 + z)^2*(2 + z))/(-2 + 3*z) */
/* a=2 is the semi-basilica */

complex p3_3(z)
	complex z;
{
	static complex one={1.0,0.0},
		two={2.0,0.0}, three={3.0,0.0};
	complex a, numer, denom;

	a = sub(z,one);
	a = mult(a,a);
	a = mult(a,add(z,two));
	numer= mult(a,param1);
	denom= sub(mult(three,z),two);
	return(divide(numer,denom));
}
static complex p3_3cps[2]={0,0,INFINITY,0};
static int p3_3nc = 2;

/* z^2 + a/z^3  */
/* J is a Cantor set x S^1 for a = 10^{-10} */
/* We suppress the finite critical points */

complex can_5(z)
	complex z;
{	
	return(add(mult(z,z),
		divide(param1,mult(mult(z,z),z))));
}
static complex can_5cps[1]={INFINITY,0.0};
static int can_5nc = 1;

/* (z^2+1)/(a z^2 + 1) */
/* Matings with basilica */
complex bas_2(z)
	complex z;
{
	static complex one={1.0,0.0};
	complex z2;
	z2 = mult(z,z);
	z = divide(sub(z2,one),
		add(mult(param1,z2),one)
		);
	return(z);
}

static complex bas_2cps[2]={0.0,0.0,INFINITY,0.0};
static int bas_2nc = 2;

/* 3 (z-1)^3 (z+3)/(6 z^2 - 8z + 3) */

complex p3_4(z)
	complex z;	
{
	static complex one={1.0,0}, 
		three={3.0,0}, six={6.0,0.0}, 
		eight={8.0,0};

	complex a, a3, b, c, denom, numer;

	a = sub(z,one);
	a3 = mult(a,mult(a,a));
	b = add(z,three);
	numer = mult(three,mult(a3,b));

	a = mult(six,mult(z,z));
	b = mult(eight,z);
	denom = add(sub(three,b),a);

	return(divide(numer,denom));
}

static complex p3_4cps[3]={0.0,0.0,1.0,0.0,INFINITY,0.0};
static int p3_4nc = 3;

/* s3_4(z) = (z^4 + 2z)/(2z^3 + 1) */

complex s3_4(z)
	complex z;
{
	static complex one={1,0}, two={2,0};
	complex z3;

	z3=mult(z,mult(z,z));
	return(divide(
		mult(z,add(z3,two)),
		add(one,mult(two,z3))
		));
}

static complex s3_4cps[3]={
	{1.0,0.0},
	{-0.5, 0.8660254037844386},
	{-0.5,-0.8660254037844386}
	};
static int s3_4nc = 3;

/* s3z2_4(z) = (z^4 + 4z - 6)/(2z^3 - 6z + 5) */

complex s3z2_4(z)
	complex z;
{
	complex z3, z4;
	static complex c_6={-6,0}, c4={4,0}, c6={6,0}, c2={2,0},
		c5={5,0};

	z3=mult(z,mult(z,z));
	z4=mult(z,z3);
	return(divide(
		add(z4,sub(mult(c4,z),c6)),
		add(mult(c2,z3),add(mult(z,c_6),c5))
		));
}

static complex s3z2_4cps[4]={
	{-1.0,0.0},
	{ 2.0,0.0},
	{ 0.732050807568877, 0.0},
	{-2.732050807568877, 0.0}
	};
static int s3z2_4nc = 4;

/* s3z3_4(z) = (z^4+6z^3+12z^2+8z)/(8z^3+12z^2+6z+1) */

complex s3z3_4(z)
	complex z;
{
	complex num, den, z2, z3, z4;
	static complex c1={1,0}, c6={6,0}, c8={8,0}, c12={12,0};

	z2=mult(z,z);
	z3=mult(z,z2);
	z4=mult(z,z3);
	num=add(add(z4,mult(c6,z3)),
		add(mult(c12,z2),mult(c8,z))
	);
	den=add(add(mult(c8,z3),mult(c12,z2)),
		add(mult(c6,z),c1)
	);
	return(divide(num,den));
}

static complex s3z3_4cps[2]={
	{-2.0,0.0},
	{ 1.0,0.0}
	};
static int s3z3_4nc = 2;

/* s3s3_4(z) = (z^4+12z^3+120z^2+316z+234)/
		(8z^3+60z^2+150z+125) */

complex s3s3_4(z)
	complex z;
{
	complex num, den, z2, z3, z4;
	static complex c12={12,0}, c120={120,0}, c316={316,0},
	c234={234,0}, c8={8,0}, c60={60,0}, c150={150,0},
	c125={125,0};

	z2=mult(z,z);
	z3=mult(z,z2);
	z4=mult(z,z3);
	num=add(c234,add(
		add(z4,mult(c12,z3)),
		add(mult(c120,z2),mult(c316,z))
		));
	den=add(add(mult(c8,z3),mult(c60,z2)),
		add(mult(c150,z),c125)
	);
	return(divide(num,den));
}

static complex s3s3_4cps[3]={
	{2.0,0.0},
	{-1.803847577293369, 0.0},
	{-12.19615242270663, 0.0}
	};
static int s3s3_4nc = 3;

/* a5_11(z) = (z^11 + 66z^6 -11z)/(-11z^10-66z^5+1) */
complex a5_11(z)
	complex z;
{
	complex z2, z4, z5, z6, z10, z11;
	complex numer, denom;
	complex add(), divide(), mult(); 

	static complex p66 = {66.0,0.0}, m66 = {-66.0,0.0},
		m11 = {-11.0, 0.0}, p1 = {1.0,0.0};

	z2 = mult (z, z);
	z4 = mult (z2, z2);
	z5 = mult (z, z4);
	z6 = mult (z2, z4);
	z10 = mult (z5, z5);
	z11 = mult (z10, z);
	numer = add(z11,add(mult (z6,p66),mult(z,m11)));
	denom = add(p1,add(mult(z5,m66),mult(z10,m11)));	
	return (divide (numer, denom));
}

static complex a5_11cps[20] ={
	{-2.391693058376407, -1.737666720845667}, 
	{-2.391693058376406, 1.737666720845668}, 
	{-0.978147600733806, -0.7106658305995174}, 
	{-0.978147600733806, 0.7106658305995177}, 
	{-0.827090915285201, 0}, 
	{-0.3382612127175641, 0}, 
	{-0.255585148716257, -0.786610204550513}, 
	{-0.2555851487162569, 0.786610204550513}, 
	{-0.1045284632676064, -0.3217055305649404}, 
	{-0.1045284632676063, 0.3217055305649404}, 
	{0.2736590696263884, -0.1988249522579513}, 
	{0.2736590696263884, 0.1988249522579513}, 
	{0.3736191374661519, -1.149881468553194}, 
	{0.3736191374661521, 1.149881468553194}, 
	{0.6691306063588573, -0.4861518423097244}, 
	{0.6691306063588574, 0.4861518423097242}, 
	{0.913545457642601, -2.811603815447866}, 
	{0.913545457642601, 2.811603815447866}, 
	{1.209056926535307, 0}, 
	{2.956295201467612, 0}
	};
static int a5_11nc = 20;

/* a5_19(z) = 
(-57z^15+247z^10+171z^5+1)/(-z^19+171z^14-247z^9-57z^4)
*/

complex a5_19(z)
	complex z;
{
	complex z2, z4, z5, z9, z10, z14, z15, z19;
	complex numer, denom;
	complex add(), divide(), mult(); 

	static complex m57 = {-57.0,0.0}, p247 = {247.0,0.0},
		p171 = {171.0, 0.0}, m247 = {-247.0,0.0},
		p1 = {1.0,0.0}, m1 = {-1.0,0.0};

	if(z.x == 0 && z.y == 0)
	{	z.x = INFINITY;
		z.y = 0.0;
		return(z);
	}
	if(z.x == INFINITY && z.y == 0)
	{	z.x = z.y = 0.0;
		return(z);
	}
	z2 = mult (z, z);
	z4 = mult (z2, z2);
	z5 = mult (z, z4);
	z9 = mult (z4, z5);
	z10 = mult (z5, z5);
	z14 = mult (z10, z4);
	z15 = mult (z10, z5);
	z19 = mult (z10, z9);
	numer = add(
		add(p1, mult (z5,p171)),
		add(mult(z10,p247), mult(z15,m57)));
	denom = add(
		add(mult(m57,z4),mult(m247,z9)),
		add(mult(p171,z14),mult(m1,z19)));
	return (divide (numer, denom));
}

static complex a5_19cps[12] = {
	{-1.618033988749895, 0}, 
	{-0.5000000000000002, -0.3632712640026805}, 
	{-0.5000000000000002, 0.3632712640026807},
	{-0.5, -1.538841768587627}, 
	{-0.4999999999999998, 1.538841768587627}, 
	{0.1909830056250525, -0.5877852522924735}, 
	{0.1909830056250527, 0.5877852522924734}, {0.6180339887498951, 0}, 
	{1.309016994374947, -0.951056516295154}, 
	{1.309016994374948, 0.951056516295153}, 
	{0, 0},
	{INFINITY,0}
	};
static int a5_19nc = 12;

/* a5_29(z) = (1-435z^5-6670z^10-3335z^20+87z^25)/
	    (87z^4+3335z^9+6670z^19-435z^24-z^29)
*/
complex a5_29(z)
complex z;
{
	complex z2, z4, z5, z9, z10, z19,z20,
		z24, z25, z29;
	complex numer, denom;
	complex add(), divide(), mult(), sub(); 
	static complex     p1 = {     1.0,0.0};
	static complex    p87 = {    87.0,0.0};
	static complex   m435 = {  -435.0,0.0};
	static complex  m3335 = { -3335.0,0.0};
	static complex  p3335 = {  3335.0,0.0};
	static complex  m6670 = { -6670.0,0.0};
	static complex  p6670 = {  6670.0,0.0};

	if(z.x == 0 && z.y == 0)
	{	z.x = INFINITY;
		z.y = 0.0;
		return(z);
	}
	if(z.x == INFINITY && z.y == 0)
	{	z.x = z.y = 0.0;
		return(z);
	}
	z2 = mult (z, z);
	z4 = mult (z2, z2);
	z5 = mult (z, z4);
	z9 = mult (z4, z5);
	z10 = mult (z5, z5);
	z19 = mult (z10, z9);
	z20 = mult (z10, z10);
	z24 = mult (z20, z4);
	z25 = mult (z20, z5);
	z29 = mult (z20, z9);
	numer = add(p1,add(
		add(mult(z5,m435), mult(z10,m6670)),
		add(mult(z20,m3335), mult(z25,p87))));
	denom = sub(
		add(add(mult(z24,m435), mult(z19,p6670)),
		add(mult(z9,p3335), mult(z4,p87))),
		z29);
	return (divide (numer, denom));
}

static complex a5_29cps[] = {
	{0, 0}, 
	{INFINITY,0},
	{-2.391693058376407, -1.737666720845667}, 
	{-2.391693058376406, 1.737666720845668}, 
	{-1.618033988749895, 0}, 
	{-0.978147600733806, -0.7106658305995174}, 
	{-0.978147600733806, 0.7106658305995177},
	{-0.827090915285201, 0}, 
	{-0.5000000000000002, -0.3632712640026805}, 
	{-0.5000000000000002, 0.3632712640026807},
	{-0.5, -1.538841768587627}, 
	{-0.4999999999999998, 1.538841768587627},
	{-0.3382612127175641, 0}, 
	{-0.255585148716257, -0.786610204550513}, 
	{-0.2555851487162569, 0.786610204550513}, 
	{-0.1045284632676064, -0.3217055305649404}, 
	{-0.1045284632676063, 0.3217055305649404}, 
	{0.1909830056250525, -0.5877852522924735}, 
	{0.1909830056250527, 0.5877852522924734}, 
	{0.2736590696263884, -0.1988249522579513}, 
	{0.2736590696263884, 0.1988249522579513}, 
	{0.3736191374661519, -1.149881468553194}, 
	{0.3736191374661521, 1.149881468553194},
	{0.6180339887498951, 0}, 
	{0.6691306063588573, -0.4861518423097244}, 
	{0.6691306063588574, 0.4861518423097242}, 
	{0.913545457642601, -2.811603815447866}, 
	{0.913545457642601, 2.811603815447866},
	{1.209056926535307, 0}, 
	{1.309016994374947, -0.951056516295154}, 
	{1.309016994374948, 0.951056516295153},
	{2.956295201467612, 0}
	};

static int a5_29nc=32;

complex v4_15(z)
	complex z;
{	
	static complex p1875={1875, 0};
	static complex m24375={-24375, 0};
	static complex m1125={-1125, 0};
	static complex m104775={-104775, 0};
	static complex p7425={7425, 0};
	static complex p22275={22275, 0};
	static complex p495={495, 0};
	static complex p891={891, 0};

	static complex m4125={-4125, 0};
	static complex m1375={-1375, 0};
	static complex m37125={-37125, 0};
	static complex p62865={62865, 0};
	static complex m405={-405, 0};
	static complex p5265={5265, 0};
	static complex p243={243, 0};

	complex z2, z3, z4, z5, z6, z7, z8;
	complex z9, z10, z11, z12, z13, z14, z15;
	complex numer, denom;

	z2=mult(z,z);
	z3=mult(z2,z);
	z4=mult(z3,z);
	z5=mult(z4,z);
	z6=mult(z5,z);
	z7=mult(z6,z);
	z8=mult(z7,z);
	z9=mult(z8,z);
	z10=mult(z9,z);
	z11=mult(z10,z);
	z12=mult(z11,z);
	z13=mult(z12,z);
	z14=mult(z13,z);
	z15=mult(z14,z);

	numer = 
	add( add( 
		add(p1875,
		    mult(m24375,z2)),
		add(mult(m1125,z4),
		    mult(m104775,z6)) 
		),
	     add(
		add(mult(p7425,z8),
		    mult(p22275,z10)),
		add(mult(p495,z12),
		    mult(p891,z14))
		)
	);

	denom = 
	add( add( 
		add(mult(m4125,z),
		    mult(m1375,z3)),
		add(mult(m37125,z5),
		    mult(p7425,z7)) 
		),
	     add(
		add(mult(p62865,z9),
		    mult(m405,z11)),
		add(mult(p5265,z13),
		    mult(p243,z15))
		)
	);

	return(divide(numer,denom));
	}

static complex v4_15cps[16] = {
	{-3.188099401672292, 0}, 
	{-0.803428418944652, -0.803428418944652},
	{-0.803428418944652, 0.803428418944652}, 
	{-0.7832270778663048, -2.050515110763247}, 
	{-0.7832270778663048, 2.050515110763247}, 
	{-0.5494338679968265, 0.2098650630044637}, 
	{-0.5494338679968256, -0.209865063004464},
	{0., -0.4049417179585518}, 
	{0., 0.4049417179585518},
	{0.5494338679968256, 0.209865063004464}, 
	{0.5494338679968265, -0.2098650630044637}, 
	{0.7832270778663048, -2.050515110763247}, 
	{0.7832270778663048, 2.050515110763247}, 
	{0.803428418944652, -0.803428418944652}, 
	{0.803428418944652, 0.803428418944652},
	{3.188099401672292, 0}
	};

static int v4_15nc = 16;

/* Newton's method for 4z^3 - 2z + 1:
   f[z] = (8z^3-1)/(16z^2-2)  */
complex new_3(z)
	complex z;
{
	complex z2, z3;
	static complex p8={8.0,0.0}, m1={-1.0,0.0},
		p16={16.0,0.0}, m2 = {-2.0,0.0};
	
	z2=mult(z,z);
	z3=mult(z2,z);
	return(
	divide(add(mult(p8,z3),m1),
		add(mult(p16,z2),m2)
	));
}

static complex new_3cps[4]={
	{0.0,0.0},
	{0.442323088559658, -0.2948714025111029},
	{0.442323088559658,  0.2948714025111029},
	{0.884646177119316, 0.0}
	};
static int new_3nc = 4;

/* a exp(z) */
complex aexp(z)
	complex z;
{
	return(mult(param1,cx_exp(z)));
}

static complex expcps[1]={{0.0,0.0}};
static int expnc=1;

/* a sin(z) */
complex a_sin(z)
	complex z;
{
	return(mult(param1,cx_sin(z)));
}

static complex sincps[2]={{PIVAL/2,0.0},{-PIVAL/2,0.0}};
static int sinnc=2;

/* a z^2 (z-3)/(1-3z) */
complex rot_3(z)
	complex z;
{
	complex w;
	static complex p1={1.0,0.0}, p3 = {3.0,0.0};

	w = divide(sub(z,p3),sub(p1,mult(p3,z)));
	w = mult(mult(z,z),w);
	w = mult(param1,w);
	return(w);
}
static int rot_3nc=3;
static complex rot_3cps[3]=
	{{0.0,0.0},{1.0,0.0},{INFINITY,0}};

/* z (1+z)^3 
   cp of mult 3 at -1 lands on parabolic;
   other cp at -1/4 in basin.
*/
complex pet_4(z)
	complex z;
{
	static complex p1={1.0,0.0};
	complex w, w2;

	w = add(z,p1);
	w2 = mult(w,w);
	w  = mult(z,mult(w2,w));
	return(w);
}
static int pet_4nc=2;
static complex pet_4cps[2]=
	{{-1.0,0.0},{INFINITY,0}};

/* Newton's method for p(z)=(z+1)(z-1)(z-a)
   f(z) = z - p(z)/p'(z).
   cps at zeros and at roots.
*/
complex ncub_3(z)
	complex z;
{
	static complex p1={1.0,0.0}, p2={2.0,0.0}, p3={3.0,0.0};
	complex p3z2, p2az, p, pprime, w;

	p = mult(add(z,p1),mult(sub(z,p1),sub(z,param1)));
	p3z2 = mult(p3,mult(z,z));
	p2az = mult(p2,mult(param1,z));
	pprime = sub(p3z2,add(p1,p2az));
	w = sub(z,divide(p,pprime));
	return(w);
}
static int ncub_3nc=4;
static complex ncub_3cps[4]=
	{{0.0,0.0},{0.0,0.0},{-1.0,0.0},{1.0,0.0}};
/* Cps are (0+1+a)/3 and a,-1,1. */
ncub_3set()
{
	ncub_3cps[1] = param1;
	ncub_3cps[0].x = param1.x / 3;
	ncub_3cps[0].y = param1.y / 3;
}

/* Napolean mapping -- stellating simplices */
/*      f(z) = z(z+a^2-2a)/(z-2az+a^2)      */
complex nap_2(z)
	complex z;
{
	static complex two = {2.0,0.0};
	complex denom, num, p2, w;

	p2 = mult(param1,param1);
	num = sub(add(z,p2),mult(two,param1));
	num = mult(z,num);
	denom = sub(add(z,p2),mult(two,mult(param1,z)));
	w = divide(num,denom);
	return(w);
}

static int nap_2nc=2;
complex nap_2cps[2];

/* Cps are roots of (1-2a)z^2 + (2a^2) z + a^4-2a^3 = 0 */
nap_2set()
{
	complex a, b, c, p2, p3, p4;
	complex two = {2.0,0.0};

	p2 = mult(param1,param1);
	p3 = mult(param1,p2);
	p4 = mult(param1,p3);

	a.x = 1 - 2*param1.x; 
	a.y = 0 - 2*param1.y; 
	b = mult(two,p2);
	c = sub(p4,mult(two,p3));
	quadsol(a,b,c,&nap_2cps[0],&nap_2cps[1]);
}

/* Quadratic Bicritical map f(z) = b/(z + 1/z + a) */
/* f(0)=0, f'(0)=b; cps at +/- 1 */
complex bic_2(z)
	complex z;
{
	complex w;

	w = add(z,recip(z));
	w = add(w,param1);
	w = divide(param2,w);
	return(w);
}

static int bic_2nc=2;
static complex bic_2cps[2] = {-1.0,0,1.0,0};

/* Cubic Bicritical map f(z) = b/(1/z + 8/(z+3/z) + a) */
/* f(0)=0, f'(0)=b; cps at +/- 1 */
complex bic_3(z)
	complex z;
{
	complex w;
	complex three = {3.0,0.0}, eight = {8.0,0.0};

	w = add(z,divide(three,z));
	w = divide(eight,w);
	w = add(w,recip(z));
	w = add(w,param1);
	w = divide(param2,w);
	return(w);
}

static int bic_3nc=2;
static complex bic_3cps[2] = {-1.0,0,1.0,0};

/* Quartic Bicritical map f(z) = b/((1+6z^2+z^4)/(z+z^3)+a)
/* f(0)=0, f'(0)=b; cps at +/- 1 */
complex bic_4(z)
	complex z;
{
	complex den, num, w, z2, z3, z4;
	complex one={1.0,0.0}, six = {6.0,0.0};

	z2 = mult(z,z);	
	z3 = mult(z2,z);
	z4 = mult(z3,z);
	num = add(z4,add(one,mult(six,z2)));
	den = add(z,z3);
	w = add(divide(num,den),param1);
	w = divide(param2,w);
	return(w);
}

static int bic_4nc=2;
static complex bic_4cps[2] = {-1.0,0,1.0,0};

/* Parabolic rational map f(z) = z/(1/b + z/a + z^2) */
/* f(0)=0, f'(0)=b */
complex par_2(z)
	complex z;
{
	static complex one = {1,0};
	complex den, num, z2, w;

	z2 = mult(z,z);
	den = divide(z,param1);
	den = add(den,recip(param2));
	den = add(den,z2);
	num = z;
	w = divide(num,den);
	return(w);
}

static int par_2nc=2;
complex par_2cps[2];

par_2set()
{
	complex z;
	static complex minusone = {-1,0};

	par_2cps[0] = cx_sqrt(recip(param2));
	par_2cps[1] = mult(minusone,par_2cps[0]);
}

/* Degree two rational map f(z) = exp(a)/(z+1/z+b) */
/* Preset param1 to exp(a) */
complex esc_2(z)
	complex z;
{
	static complex one = {1,0};
	complex w;

	w = divide(param1,add(z,add(recip(z),param2)));
	return(w);
}
static int esc_2nc=2;
static complex esc_2cps[2]={1,0,-1,0};

/* Degree three rational map f(z) = z^2 (2+3d - (2d+1) z)/(d+z) */
/* Fixed cps at 0,1,inf */
complex fix_3(z)
	complex z;
{
	complex den,num,w;
	static complex one = {1,0}, two = {2,0}, three = {3,0};

	num = add(two,mult(three,param1));
	num = sub(num,mult(z,add(mult(two,param1),one)));
	num = mult(num,mult(z,z));
	den = add(z,param1);
	w = divide(num,den);
	return(w);
}
static int fix_3nc=4;
static complex fix_3cps[4]={{0,0},{1,0},{INFINITY,0},{0,0}};
fix_3set()
{
	complex den, num, sq;

	sq = mult(param1,param1);
	
	num.x = -2*param1.x - 3*sq.x; 
	num.y = -2*param1.x - 3*sq.x; 
	den.x = 1 + 2*param1.x;
	den.y = 1 + 2*param1.y;
	fix_3cps[3] = divide(num,den);
}

/* ATTRACTORS AND CRITICAL POINTS */

critical_points(acps,anc)
	complex *acps[];
	int *anc;
{
	*anc = nc;
	*acps = cps;
}

complex iter(z)
	complex z;
{
	complex w;
	w = (*f)(z);
	if(cx_abs(w) > 1/ATT_EPS)
	if((*f) != aexp && (*f) != a_sin)  
	{	w.x = INFINITY;
		w.y=0.0;}
	return(w);
}

/* One attractor for each critical point */
/* Insure attractors agree if in same cycle */
attractor(at, na)
	complex at[];
	int *na;
{
	int i, j, k;
	complex z;

	for(i=0; i<nc; i++)
	{	z = cps[i];
		for(j=0; j<att_iter; j++) z = iter(z);
		at[i] = z;
		for(j=0; j<att_iter; j++) 
		{	z = iter(z);
			for(k=0; k<i; k++)
			if(cx_abs(sub(z,at[k])) < ATT_EPS) at[i] = at[k];
		}
	}
	*na = nc;
}

/* MAP SELECTOR */

/* Sets pointer to function and critical points */
/* Returns 1 if name not known */

pick_map(name,p1,p2,p3,p4,p5,p6)
	char *name;
	complex p1, p2, p3, p4, p5, p6;
{
	if(!strcmp(name,"POLY_2"))
	{	f = poly_2;
		param1 = p1;
		cps = poly_2cps;
		nc  = poly_2nc;
		return(0);
	}

	if(!strcmp(name,"RAT_2"))
	{	f = rat_2;
		param1 = p1;
		param2 = p2;
		param3 = p3;
		param4 = p4;
		param5 = p5;
		param6 = p6;
		rat_2set();
		cps = rat_2cps;
		nc  = rat_2nc;
		return(0);
	}

	if(!strcmp(name,"Z3_3"))
	{	f = z3_3;
		param1 = p1;
		cps = z3_3cps;
		nc  = z3_3nc;
		return(0);
	}

	if(!strcmp(name,"Z4_4"))
	{	f = z4_4;
		param1 = p1;
		cps = z4_4cps;
		nc  = z4_4nc;
		return(0);
	}

	if(!strcmp(name,"FCP_3"))
	{	f = fcp_3;
		param1 = p1;
		cps = fcp_3cps;
		fcp_3set();
		nc  = fcp_3nc;
		return(0);
	}

	if(!strcmp(name,"FP_3"))
	{	f = fp_3;
		param1 = p1;
		param2 = p2;
		cps = fp_3cps;
		fp_3set();
		nc  = fp_3nc;
		return(0);
	}

	if(!strcmp(name,"BAS_2"))
	{	f = bas_2;
		param1 = p1;
		cps = bas_2cps;
		nc  = bas_2nc;
		return(0);
	}

	if(!strcmp(name,"CAN_5"))
	{	f = can_5;
		param1 = p1;
		cps = can_5cps;
		nc  = can_5nc;
		return(0);
	}

	if(!strcmp(name,"P3_3"))
	{	f = p3_3;
		param1 = p1;
		cps = p3_3cps;
		nc  = p3_3nc;
		return(0);
	}

	if(!strcmp(name,"Q3_3"))
	{	f = q3_3;
		param1 = p1;
		cps = q3_3cps;
		nc  = q3_3nc;
		return(0);
	}

	if(!strcmp(name,"P3_4"))
	{	f = p3_4;
		cps = p3_4cps;
		nc  = p3_4nc;
		return(0);
	}

	if(!strcmp(name,"S3_4"))
	{	f = s3_4;
		cps = s3_4cps;
		nc  = s3_4nc;
		return(0);
	}

	if(!strcmp(name,"S3Z2_4"))
	{	f = s3z2_4;
		cps = s3z2_4cps;
		nc  = s3z2_4nc;
		return(0);
	}

	if(!strcmp(name,"S3Z3_4"))
	{	f = s3z3_4;
		cps = s3z3_4cps;
		nc  = s3z3_4nc;
		return(0);
	}

	if(!strcmp(name,"S3S3_4"))
	{	f = s3s3_4;
		cps = s3s3_4cps;
		nc  = s3s3_4nc;
		return(0);
	}

	if(!strcmp(name,"A5_11"))
	{	f = a5_11;
		cps = a5_11cps;
		nc  = a5_11nc;
		return(0);
	}

	if(!strcmp(name,"A5_19"))
	{	f = a5_19;
		cps = a5_19cps;
		nc  = a5_19nc;
		return(0);
	}

	if(!strcmp(name,"A5_29"))
	{	f = a5_29;
		cps = a5_29cps;
		nc  = a5_29nc;
		return(0);
	}

	if(!strcmp(name,"V4_15"))
	{	f = v4_15;
		cps = v4_15cps;
		nc  = v4_15nc;
		return(0);
	}

	if(!strcmp(name,"NEW_3"))
	{	f = new_3;
		cps = new_3cps;
		nc  = new_3nc;
		return(0);
	}

	if(!strcmp(name,"EXP"))
	{	f = aexp;
		cps = expcps;
		nc  = expnc;
		param1=p1;
		return(0);
	}

	if(!strcmp(name,"SIN"))
	{	f = a_sin;
		cps = sincps;
		nc  = sinnc;
		param1=p1;
		return(0);
	}

	if(!strcmp(name,"ROT_3"))
	{	f = rot_3;
		cps = rot_3cps;
		nc  = rot_3nc;
		param1=p1;
		return(0);
	}

	if(!strcmp(name,"PET_4"))
	{	f = pet_4;
		cps = pet_4cps;
		nc  = pet_4nc;
		return(0);
	}

	if(!strcmp(name,"NCUB_3"))
	{	f = ncub_3;
		param1=p1;
		cps = ncub_3cps;
		nc  = ncub_3nc;
		ncub_3set();
		return(0);
	}

	if(!strcmp(name,"NAP_2"))
	{	f = nap_2;
		param1=p1;
		cps = nap_2cps;
		nc  = nap_2nc;
		nap_2set();
		return(0);
	}

	if(!strcmp(name,"HALF_2"))
	{	f = bic_2;
		param1=p1;
		param2=p2;
		cps = bic_2cps;
		nc  = bic_2nc;
		return(0);
	}

	if(!strcmp(name,"BIC_2"))
	{	f = bic_2;
		param1=p1;
		param2=p2;
		cps = bic_2cps;
		nc  = bic_2nc;
		return(0);
	}

	if(!strcmp(name,"BIC_3"))
	{	f = bic_3;
		param1=p1;
		param2=p2;
		cps = bic_3cps;
		nc  = bic_3nc;
		return(0);
	}

	if(!strcmp(name,"BIC_4"))
	{	f = bic_4;
		param1=p1;
		param2=p2;
		cps = bic_4cps;
		nc  = bic_4nc;
		return(0);
	}

	if(!strcmp(name,"PAR_2"))
	{	f = par_2;
		param1=p1;
		param2=p2;
		cps = par_2cps;
		nc  = par_2nc;
		par_2set();
		return(0);
	}

	if(!strcmp(name,"ESC_2"))
	{	f = esc_2;
		param1=p1;
		param2=p2;
		param1=cx_exp(param1);
		cps = esc_2cps;
		nc  = esc_2nc;
		return(0);
	}

	if(!strcmp(name,"FIX_3"))
	{	f = fix_3;
		param1=p1;
		cps = fix_3cps;
		nc  = fix_3nc;
		fix_3set();
		return(0);
	}

	return(1);
}

/* Scans argument list for map and parameters */
arg_map(argc, argv, i)
	int argc;
	char *argv[];
	int *i;
{
	complex q1, q2, q3, q4, q5, q6;
	char name[20];
	if(argc <= (1+*i)) return(1);
	sscanf(argv[++*i],"%s",name);

/* These need to be scanned for one parameter */
	if(0==  strcmp(name,"POLY_2")*
		strcmp(name,"Z3_3")*
		strcmp(name,"Z4_4")*
		strcmp(name,"FCP_3")*
		strcmp(name,"BAS_2")*
		strcmp(name,"P3_3")*
		strcmp(name,"Q3_3")*
		strcmp(name,"CAN_5")*
		strcmp(name,"EXP")*
		strcmp(name,"SIN")*
		strcmp(name,"ROT_3")*
		strcmp(name,"NCUB_3")*
		strcmp(name,"NAP_2")*
		strcmp(name,"FIX_3")
	) { 	if(argc <= (2+*i)) return(1);
		sscanf(argv[++*i],"%lf",&q1.x);
		sscanf(argv[++*i],"%lf",&q1.y);
	}

/* These need to be scanned for two parameters */
	if(0==  strcmp(name,"FP_3")*
		strcmp(name,"PAR_2")*
		strcmp(name,"HALF_2")*
		strcmp(name,"BIC_2")*
		strcmp(name,"BIC_3")*
		strcmp(name,"BIC_4")*
		strcmp(name,"ESC_2")
	) { 	if(argc <= (4+*i)) return(1);
		sscanf(argv[++*i],"%lf",&q1.x);
		sscanf(argv[++*i],"%lf",&q1.y);
		sscanf(argv[++*i],"%lf",&q2.x);
		sscanf(argv[++*i],"%lf",&q2.y);
	}

/* These need to be scanned for six parameters */
	if(0==  strcmp(name,"RAT_2")
	) { 	if(argc <= (10+*i)) return(1);
		sscanf(argv[++*i],"%lf",&q1.x);
		sscanf(argv[++*i],"%lf",&q1.y);
		sscanf(argv[++*i],"%lf",&q2.x);
		sscanf(argv[++*i],"%lf",&q2.y);
		sscanf(argv[++*i],"%lf",&q3.x);
		sscanf(argv[++*i],"%lf",&q3.y);
		sscanf(argv[++*i],"%lf",&q4.x);
		sscanf(argv[++*i],"%lf",&q4.y);
		sscanf(argv[++*i],"%lf",&q5.x);
		sscanf(argv[++*i],"%lf",&q5.y);
		sscanf(argv[++*i],"%lf",&q6.x);
		sscanf(argv[++*i],"%lf",&q6.y);
	}

	return(pick_map(name,q1,q2,q3,q4,q5,q6));
}

/* Explain how to select a function */
fnusage()
{
        fprintf(stderr,"  -f [name][parameters] Select function\n");
	fprintf(stderr,"  Possible names:\n");
	fprintf(stderr,"  With six parameters:  RAT_2\n");
	fprintf(stderr,"  With two parameters: FP_3, PAR_2, BIC_2, BIC_3, BIC_4, ESC_2\n");
	fprintf(stderr,"  With one parameter:\n");
	fprintf(stderr,"  POLY_2 BAS_2 FCP_3 Z3_3 Z4_4 P3_3 Q3_3 CAN_5 EXP SIN ROT_3 NCUB_3 NAP_2 FIX_3\n");
	fprintf(stderr,"  With no parameters:\n");
	fprintf(stderr,"  P3_4 S3_4 S3Z2_4 S3Z3_4 S3S3_4 A5_11 A5_19 A5_29 V4_15\n");
	fprintf(stderr,"  NEW_3 PET_4\n");
}

set_att_iter(i)
	int i;
{
	att_iter = i;
}

/* Solve quadratic equation a r^2 + b r + c =  0 */
quadsol(a,b,c,r1,r2)
	complex a,b,c,*r1,*r2;
{
	complex disc;
	static complex two = {2.0,0.0}, four = {4.0,0.0};

	disc=sub(mult(b,b),mult(four,mult(a,c)));
	disc=cx_sqrt(disc);
	*r1 = divide(sub(disc,b),mult(two,a));
	disc.x = -disc.x;  disc.y = -disc.y;
	*r2 = divide(sub(disc,b),mult(two,a));
}
