/////////////////////////////////////
// ftloo01.cpp
// 
// stephen c. ferguson
// created 10/18/98
// http://www.eclectasy.com/Iterations-et-Flarium24/
// modified 11-18-00
// modified 08-11-02

/////////////////////////////////////

#define STRICT
#define WIN32_LEAN_AND_MEAN
#include <windows.h>

#include "fll01.h"

void _formula(CDATA *v);

void Formula_00(CDATA *v);
void Formula_01(CDATA *v);
void Formula_02(CDATA *v);
void Formula_03(CDATA *v);
void Formula_04(CDATA *v);
void Formula_05(CDATA *v);
void Formula_06(CDATA *v);
void Formula_07(CDATA *v);
void Formula_08(CDATA *v);
void Formula_09(CDATA *v);
void Formula_10(CDATA *v);
void Formula_11(CDATA *v);
void Formula_12(CDATA *v);
void Formula_13(CDATA *v);
void Formula_14(CDATA *v);
void Formula_15(CDATA *v);
void Formula_16(CDATA *v);
void Formula_17(CDATA *v);
void Formula_18(CDATA *v);
void Formula_19(CDATA *v);
void Formula_20(CDATA *v);
void Formula_21(CDATA *v);
void Formula_22(CDATA *v);
void Formula_23(CDATA *v);
void Formula_24(CDATA *v);
void Formula_25(CDATA *v);
void Formula_26(CDATA *v);
void Formula_27(CDATA *v);
void Formula_28(CDATA *v);
void Formula_29(CDATA *v);
void Formula_30(CDATA *v);
void Formula_31(CDATA *v);
void Formula_32(CDATA *v);
void Formula_33(CDATA *v);
void Formula_34(CDATA *v);
void Formula_35(CDATA *v);
void Formula_36(CDATA *v);
void Formula_37(CDATA *v);
void Formula_38(CDATA *v);
void Formula_39(CDATA *v);
void Formula_40(CDATA *v);
void Formula_41(CDATA *v);
void Formula_42(CDATA *v);
void Formula_43(CDATA *v);
void Formula_44(CDATA *v);
void Formula_45(CDATA *v);
void Formula_46(CDATA *v);
void Formula_47(CDATA *v);
void Formula_48(CDATA *v);


#define FORMULA_TOTAL 49

typedef struct
{   // export pointer to address DLL
	int count;    // number of formulas in DLL
} TOTAL;

TOTAL Total;

////////////////////////////////////////////

typedef struct
{    // info structure
  void (*action)			(CDATA *v);		// action iteration
  char *title;											// for selection purposes
  char *description;								// identification string(allocated by dll)
}INFO;

////////////////////////////////////////////

INFO Info[FORMULA_TOTAL] = 
{
	Formula_00,
		"z*z+c "		,
			"Mandelbrot Method ",
	Formula_01,
		"z*z+c "		,
			"Newton Method ",
	Formula_02,
		"v->z = z-((z4-1)/(4*z3))+c; ",
			"4th Order Nova ",
	Formula_03,
		"z = z-(z3-1)/(3*z2)+c;",
			"3rd Order Nova ",
	Formula_04,
		"z = z-(z4-1)/(4*z3)+c;",
			"04 ",
	Formula_05,
		"z = z-c2*((z*z2-1)/(3*z2))+v->c;",
			"Twister Weed ",
	Formula_06,
		"z = c/ccos_z(); z=z*z;",
			"06 ",
	Formula_07,
		"Phoenix"	,
			"i=0.568",
	Formula_08,
		"z = z*z*z-aa3*z+b; ",
			"Cubic CBAP (by Holger Jaenisch)",
	Formula_09,
		"z = z-(z4-z)/(4*z3-z)+v->c;",
			"Nova Variation",
	Formula_10,
		"z = z*z*z-aa3*z+b;",
			"Cubic (CCAP by Holger Jaenisch",
	Formula_11,
		"z = z-(z4-1)/(4*z3-z)+v->c;",
			"Nova Variation ",
	Formula_12,"z = z = z-(z4-z2*csin(z))/(4*z3-c2*ccos(z))+v->c;"	,"12 ",
	Formula_13,"z = z-(z1*z1 - z1*z*c - z*z - zsin)/(4*z*z1 - 3*z1*z - 2*z - zcos) + c;"	,"13 ",
	Formula_14,"z = z = z-(z*z1*c + z1*c - zsin + z1)/(3*z1*c - 2*z*c - zcos) + c;"	,"14 ",
	Formula_15,"z = z-(c/z1*z1-z1*z*c+z1+c-zsin)/(c/4*z1*z*c-3*z1*c-2*z+1-zcos);"	,"15 ",
	Formula_16,"z = z-(z1*z1*c + z1*z + z + c)/(4*z1*z*c + 3*z1 + 1);"	,"16 ",
	Formula_17,"z = z-(z1*z - z1*c + z)/(3*z1 - 2*z*c + 1);"	,"17 ",
	Formula_18,"z = z-(z1*z1 + z1*z*c + z + c)/(4*z1*z + 3*z*z*c + 1);"	,"18 ",
	Formula_19,"z = z-(z*z4 + z2*z2*c + c)/(5*z4 + 4*z2*z*c); "	,"19 ",
	Formula_20,"z = z-(z1*z1*c + z1*z - z - c)/(4*z1*z*c + 3*z1 - 1); "	,"20 ",
	Formula_21,"z = z-(z1*z1*c-z1*z-z-c-1)/(4*z1*z*c-3*z1-1); "	,"21 ",
	Formula_22,"z = ztan*z-(zsin+z1*z1*c+z1*z)/(zcos+4*z1*z*c+3*z1);"	,"22 ",
	Formula_23,"z = z-(z1*z1*c+z1*z+z+c)/(4*z1*z*c+3*z1+1);"	,"23 ",
	Formula_24,"z = zcos*z/zsin-(z1*z1*c+z1*z+z+c)/(4*z1*z*c+3*z1+1);"	,"24 ",
	Formula_25,"z = z-(z1*z1*c+z1*z+z+c)/(4*z1*z*c+3*z1+1)+c;"	,"25 ",
	Formula_26,"z = z-ztan-(z1*z1*c-z1*z-z-c)/(4*z1*z*c-3*z1-1)+c;"	,"26 ",
	Formula_27,"z = z-(z1*z1*c-z1*z-z-zsin)/(4*z1*z*c-3*z1-1-zcos)*cexp_z();"	,"27 ",
	Formula_28,"z = z-(z2*z1*z*c+z1*c+z)/(3*z1*c+2*z*c+1);"	,"28 ",
	Formula_29,"z = z-cexp_z()/((z4-1-zsin)/(4*z3)-zcos)+c"	,"29 ",
	Formula_30,"z = z-((z4+1)/(4*z3))+c;"	,"30 ",
	Formula_31,"z = z-((z4-z)/(4*z3-1))+cexp_z()-c;"	,"31 ",
	Formula_32,"z = z-((z4-z )/(4*z3-1))+c+cexp_c();"	,"32 ",
	Formula_33,"z = z-((z3-z)/(3*z2-1))+c;"	,"33 ",
	Formula_34,"z = z-((z2-z)/(2*z-1))+c;"	,"34 ",
	Formula_35,"z = z-((z4-z3-z2-z)/(4*z3-3*z2-2*z-1))+c;"	,"35 ",
	Formula_36,"z = z-((z4*c-z)/(4*z3*c-1))+c;"	,"36 ",
	Formula_37,"z = z-((z4-z)/(4*z3-1))+c;"	,"37 ",
	Formula_38,"z = z-((z4*zsin-z)/(4*z3*zcos-1))+c;"	,"38 ",
	Formula_39,"z = z-((z3-z-zsin)/(3*z2-zcos))+c;"	,"39 ",
	Formula_40,"z = z-((z4*c*zsin-z-zsin-z)/(4*z3*c*zcos-1-zcos-1))+c;"	,"40 ",
	Formula_41,"z = z-((z4*c*zsin+z)/(4*z3*c+zcos))+c;"	,"41 ",
	Formula_42,"z = z-((z2*zsin+z)/(2*z*zsin+1))+c;"	,"42 ",
	Formula_43,"z = z-((z2*z3+z)/(2*z*z3+1))+c;"	,"43 ",
	Formula_44,"z = z-((z2*z3-z)/(2*z*z3-1))-c);"	,"44 ",
	Formula_45,"z = z-((z2*z3-z-zsin)/(2*z*z3-1-zcos))-c;"	,"45 ",
	Formula_46,"z = z-((z3*z4-z2-z)/(3*z2*z4-2*z-1))-c;"	,"46 ",
	Formula_47,"z = z-((z3*z4-z2-z)/(3*z2*z4-2*z-1))-c;"	,"47 ",		
	Formula_48,"z = z-((z4*z5-z3-z2-z-zsin)/(4*z3*z5-3*z2-2*z-1-zcos)) - c;"	,"48 "
} ;

////////////////////////////////////////////

int APIENTRY LibMain(HANDLE hInst, DWORD fdwReason, LPVOID lpReserved)
{	
  if (fdwReason == DLL_PROCESS_ATTACH)
  {
    wsprintf (cstr, "LibMain(): DLL loaded");

    Total.count = FORMULA_TOTAL;
  }
  else if (fdwReason == DLL_PROCESS_DETACH)
    wsprintf (cstr, "LibMain(): DLL unloaded");

  return TRUE;
}

////////////////////////////////////////////

int _total() { return FORMULA_TOTAL; }

//MessageBox(NULL, "step 1, formula z=z*z+c", "bug", MB_OK);

////////////////////////////////////
int _info(CDATA *v)
{
  v->StrTitle				= Info[v->formula].title;
  v->StrDescription = Info[v->formula].description;
  //MessageBox(NULL, v->StrTitle , "bug", MB_OK);
  return v->formula; 
}

////////////////////////////////////
void _formula(CDATA *v)
{
  //sprintf(cstr,"step 1, formula = %d", v->formula);
  //MessageBox(NULL, cstr, "bug", MB_OK);

  v->Iteration = -1;

  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
  {
    if (v->inverse)
    {
	//c = 1/c;
	//z = 1/z;
	
	temp        = __real__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__imag__ c  = __imag__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__real__ c  = temp;

	temp        = __real__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__imag__ z  = __imag__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__real__ z  = temp;
	
    }
  }
  else if (v->jul_save == 0)
  {
    if (v->inverse) 
    {
	//c = 1/c;
	//z = 1/z; 
	
	temp        = __real__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__imag__ c  = __imag__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__real__ c  = temp;

	temp        = __real__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__imag__ z  = __imag__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__real__ z  = temp;
	
    }
  }
  else
  {
    if (v->inverse)
    {
      //c = 1/c;
	//z = 1/z;

	temp        = __real__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__imag__ c  = __imag__ c / ((__real__ c * __real__ c) + (__imag__ c * __imag__ c));
	__real__ c  = temp;

	temp        = __real__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__imag__ z  = __imag__ z / ((__real__ z * __real__ z) + (__imag__ z * __imag__ z));
	__real__ z  = temp;
	
    }
  }
		
  v->z = z;
  v->c = c;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

  switch (v->formula)
	{
		case 0:
    	Formula_00(v);
			break;

		case 1:
    	Formula_01(v);
			break;

		case 2:
    	Formula_02(v);
			break;

		case 3:
    	Formula_03(v);
			break;

		case 4:
    	Formula_04(v);
			break;

		case 5:
    	Formula_05(v);
			break;

		case 6:
    	Formula_06(v);
			break;

		case 7:
    	Formula_07(v);
			break;

		case 8:
    	Formula_08(v);
			break;

		case 9:
    	Formula_09(v);
			break;

		case 10:
    	Formula_10(v);
			break;

		case 11:
    	Formula_11(v);
			break;

		case 12:
    	Formula_12(v);
			break;

		case 13:
    	Formula_13(v);
			break;
		
		case 14:
    	Formula_14(v);
			break;
		
		case 15:
    	Formula_15(v);
			break;

		case 16:
    	Formula_16(v);
			break;

		case 17:
    	Formula_17(v);
			break;

		case 18:
    	Formula_18(v);
			break;

		case 19:
    	Formula_19(v);
			break;

		case 20:
    	Formula_20(v);
			break;

		case 21:
    	Formula_21(v);
			break;

		case 22:
    	Formula_22(v);
			break;

		case 23:
    	Formula_23(v);
			break;

		case 24:
    	Formula_24(v);
			break;

		case 25:
    	Formula_25(v);
			break;

		case 26:
    	Formula_26(v);
			break;

		case 27:
    	Formula_27(v);
			break;

		case 28:
    	Formula_28(v);
			break;

		case 29:
    	Formula_29(v);
			break;

		case 30:
    	Formula_30(v);
			break;

		case 31:
    	Formula_31(v);
			break;

		case 32:
    	Formula_32(v);
			break;

		case 33:
    	Formula_33(v);
			break;

		case 34:
    	Formula_34(v);
			break;

		case 35:
    	Formula_35(v);
			break;

		case 36:
    	Formula_36(v);
			break;

		case 37:
    	Formula_37(v);
			break;

		case 38:
    	Formula_38(v);
			break;

		case 39:
    	Formula_39(v);
			break;

		case 40:
    	Formula_40(v);
			break;

		case 41:
    	Formula_41(v);
			break;

		case 42:
    	Formula_42(v);
			break;

		case 43:
    	Formula_43(v);
			break;

		case 44:
    	Formula_44(v);
			break;

		case 45:
    	Formula_45(v);
			break;

		case 46:
    	Formula_46(v);
			break;

		case 47:
    	Formula_47(v);
			break;

		case 48:
    	Formula_48(v);
			break;

		case 49:
    	Formula_49(v);
			break;

	}	
}

//MessageBox(NULL, "step 2, formula z=z*z+c", "bug", MB_OK);
void Formula_00(CDATA *v)
{
  //MessageBox(NULL, "step 1, formula z=z*z+c", "bug", MB_OK);

  // z=z*z+c, Mandelbrot method
  while ( v->Iteration++ < v->Maxit && 
		( mx = __real__ v->z * __real__ v->z ) +
		//( my = __imag__ v->z * __imag__ v->z ) < 1e24 )
		( my = __imag__ v->z * __imag__ v->z ) < v->Bailout )
  {
    v->pXTemp[v->Iteration] = __real__ v->z;
	  v->pYTemp[v->Iteration] = __imag__ v->z;
      
		zz = __real__ v->z * __imag__ v->z;
         __imag__ v->z = zz + zz + __imag__ v->c;
         __real__ v->z = mx - my + __real__ v->c;

			//if (filter) Do_Filter();
    //MessageBox(NULL, "step 1, formula z=z*z+c", "bug", MB_OK);
  }

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_01(CDATA *v)
{
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    //temp = __real__ v->z;
    //__real__ v->z = __imag__ v->z;
    //__imag__ v->z = temp;    
	
    z1 = v->z;

    v->z = v->z*v->z+v->c;    

    v->pXTemp[v->Iteration] = __real__ v->z;
    v->pYTemp[v->Iteration] = __imag__ v->z;
	
    zd = v->z-z1;

    temp = (sum_sqrs_zd());
  }

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
	
  //MessageBox(NULL, "step 3, formula z=z*z+c", "bug", MB_OK);
}

void Formula_02(CDATA *v)
{
  //MessageBox(NULL, "formula 02", "bug", MB_OK);
  // 4th order Nova
  if (v->jul == 0)
    v->z += 1;
  temp = 1;
  c = v->c;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z = z1 = v->z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;

    v->z = z-((z4-1)/(4*z3))+c;

    v->pXTemp[v->Iteration] = __real__ v->z;
    v->pYTemp[v->Iteration] = __imag__ v->z;

    zd = z1-v->z;
    temp = sum_sqrs_zd();
  }
}

void Formula_03(CDATA *v)
{
  // 3rd order Nova
  if (v->jul == 0)
		v->z = v->z+1;
	c = v->c;
	z = v->z;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;

    z = z-(z3-1)/(3*z2)+c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_04(CDATA *v)
{
  // Newton method
  if (v->jul == 0)
    v->z = 1;
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z = z1 = v->z;

    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    z = z-(z4-1)/(4*z3)+c;
    v->z = z*z+c;
    v->pXTemp[v->Iteration] = __real__ v->z;
    v->pYTemp[v->Iteration] = __imag__ v->z;

    zd = z1-v->z;
    temp = sum_sqrs_zd();
  }
}

void Formula_05(CDATA *v)
{
  // Twister weed method
  z = v->z;
  c2 = cexp_z();
  temp = 1;

  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;

    z = z-c2*((z*z2-1)/(3*z2))+v->c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_06(CDATA *v)
{
  temp = 1;
  c = v->c;
  z = v->z;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;

    z = c/ccos_z(); z=z*z;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z-z1;

    temp = sum_sqrs_zd();
  }
}

void Formula_07(CDATA *v)  // Phoenix
{
  __real__ c1 = 0;
  __imag__ c1 = 1;

  c = v->c;
  z = v->z;

  // Rotate it so its pointing upward
  z*=c1;
  c*=c1;

  __imag__ c = 0;

  __real__ c = 0.568;
  //__real__ c = 0.563;   // Paul Carlson's point
  //__real__ c = 0.56667; // Ushiki's / Pickover's point

  z3=0;
  temp = 1;

  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;

    z  = z*z - .5*z  + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z*z - .5*z3 + c;
    z3 = z;
    z  = z2;

    zd = z1-z3;
    temp = sum_sqrs_zd();
  }
}

void Formula_08(CDATA *v)
{
  // CBAP
	z = v->z;
  c = v->c;

  t=c;
  t3=3*t;
  t2=t*t;
  a=(t2+1)/t3;
  b=2*a*a*a+(t2-2)/t3;
  aa3=a*a*3;

  z=-a;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp < dMAX)
  {
    z1 = z;

    z = z*z*z-aa3*z+b;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_09(CDATA *v)
{
  // Nova variation method
  z = 1;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    z = z-(z4-z)/(4*z3-z)+v->c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_10(CDATA *v)
{
  // CCAP
	z = v->z;
  c = v->c;

  t=c;
  a=t;
  b=t+2*t*t*t;
  aa3=a*a*3;

  //if (jul == 0)
  z=-a;

	temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;

    z = z*z*z-aa3*z+b;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_11(CDATA *v)
{
  // Nova variation method
  z = v->z;
  z += 1;

  c2 = csin(z);
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    z = z-(z4-1)/(4*z3-c2*z)+v->c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_12(CDATA *v)
{
  // Nova variation method
  z = v->z;
  z += 1;

  c2 = csin(1/z);
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    z = z-(z4-z2*csin(z))/(4*z3-c2*ccos(z))+v->c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_13(CDATA *v)
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z - z*z*z*c - z - c;  // the function
  z'  = 4*z*z*z - 3*z*z*c - 1;      // 1st derivative
  z'' = 12*z*z - 6*z*c;             // 2nd derivative

  12*z*z = 6*z*c;
  z   = (6*z*c)/(12*z);             // solve for z
  z   = c/2;

  */

  z = v->z;
  c = csin(v->c);

  z = z/2;
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    sine_cosine();
    
    z = z-(z1*z1 - z1*z*c - z*z - zsin)/(4*z*z1 - 3*z1*z - 2*z - zcos) + c;
    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_14(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*c + z*z + z + c;  // the function
  z'  = 3*z*z*c + 2*z + 1;      // 1st derivative
  z'' = 6*z*c   + 2;            // 2nd derivative

  6*z*c = -2;
  z     = -2/(6*c);             // solve for z
  z     = -1/(3*c)

  */

  z = v->z;
  c = v->c;

  z *= 1/(3*z);

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;
    
    sine_cosine();
    
    z = z-(z*z1*c + z1*c - zsin + z1)/(3*z1*c - 2*z*c - zcos) + c;
    
    zd = z - z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_15(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z*c + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z*c + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z*c          ;  // 2nd derivative

  6*z*c = -12*z*z*c;
  z     = -(12*z*z*c)/(6*c);
  z     = -2*z*z;

  12*z*z*c = -6*z*c;

  z     = -(6*z*c)/(12*z*c);

  */

  //z     = -2*z*z;
  //z     = -(6*z*c)/(12*z*c);

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;
    z  = z2;

    //ztan = v->dStrands/(1e-15+ctan_z());
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - (c/z1*z1 - z1*z*c + z1 + c - zsin)/(c/4*z1*z*c - 3*z1*c - 2*z + 1 - zcos);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_16(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z-(z1*z1*c + z1*z + z + c)/(4*z1*z*c + 3*z1 + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_17(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z + z*z*c + c;  // the function
  z'  = 3*z*z + 2*z*c    ;  // 1st derivative
  z'' = 6*z + 2*c        ;  // 2nd derivative

  6*z = -2*c;
  z     = -2*c/6;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = z/3;
  else
    z = c/3;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z  = z - (z1*z - z1*c + z)/(3*z1 - 2*z*c + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }
  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;
  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_18(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z + z*z*z*c + c;  // the function
  z'  = 4*z*z*z + 3*z*z*c    ;  // 1st derivative
  z'' = 12*z*z  + 6*z*c      ;  // 2nd derivative

  12*z*z = -6*z*c;
  z      = -(6*z*c)/(12*z);
  z      = -c/2;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -z/2;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z  = z-(z1*z1 + z1*z*c + z + c)/(4*z1*z + 3*z*z*c + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_19(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*z + z*z*z*z*c + c;  // the function
  z'  = 5*z*z*z*z + 4*z*z*z*c    ;  // 1st derivative
  z'' = 20*z*z*z  + 12*z*z*c      ;  // 2nd derivative

	20*z*z*z = -12*z*z*c;
  z        = -(12*z*z*c)/(20*z*z);
  z        = -(6*c)/10;
  z        = -(3*c)/5;


  */


  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -(3*z)/5;

	temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z1 = z;
    z2 = z*z;
    z4 = z2*z2;			

    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z-(z*z4 + z2*z2*c + c)/(5*z4 + 4*z2*z*c);

    zd = z - z1;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_20(CDATA *v)  // testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z-(z1*z1*c + z1*z - z - c)/(4*z1*z*c + 3*z1 - 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_21(CDATA *v)  // 9th order Newton Mset
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = .5/z;

  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z-(z1*z1*c - z1*z - z - c - 1)/(4*z1*z*c - 3*z1 - 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_22(CDATA *v)  // 
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = .5/z;

  //z = cexp_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    ztan = ctan_z();
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = ztan * z - (zsin + z1*z1*c + z1*z)/(zcos + 4*z1*z*c + 3*z1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_23(CDATA *v)  // Testing again
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  //z = cexp_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - (z1*z1*c + z1*z + z + c)/(4*z1*z*c + 3*z1 + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_24(CDATA *v)
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  //z = cexp_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = zcos*z/zsin - (z1*z1*c + z1*z + z + c)/(4*z1*z*c + 3*z1 + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_25(CDATA *v)
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  c = cexp_c();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - (z1*z1*c + z1*z + z + c)/(4*z1*z*c + 3*z1 + 1) + c;

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_26(CDATA *v) 
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -.5/z;

  c = cexp_c();
  ztan = ctan_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - ztan - (z1*z1*c - z1*z - z - c)/(4*z1*z*c - 3*z1 - 1) + c;

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;

}

void Formula_27(CDATA *v)
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*z*c + z*z*z + z + c;  // the function
  z'  = 4*z*z*z*c + 3*z*z + 1    ;  // 1st derivative
  z'' = 12*z*z*c  + 6*z          ;  // 2nd derivative

  12*z*z*c = -6*z;

  z     = -(6*z)/(12*z*c);
  z     = =z/(2*z*c);
  z     = -.5/c;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = .5/z;
    //z = -3;

  //c = cexp_c();
  //ztan = ctan_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - (z1*z1*c - z1*z - z - zsin)/(4*z1*z*c - 3*z1 - 1 - zcos) * cexp_z();

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_28(CDATA *v)  // Testing
{
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z + z*z*c + z + c   ;  // the function
  z'  = 3*z*z + 2*z*c + 1         ;  // 1st derivative
  z'' = 6*z  + 2*c                ;  // 2nd derivative

  2*c = -6*z;

  z     = - (2*c)/6*z;

  */

  z = v->z;
  c = v->c;

  if (v->jul == 0)
    z = -z/3;

  //c = cexp_c();
  //ztan = ctan_z();
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    z2 = z;
    z1 = z*z;

    //ztan = ctan_z();
    //sine_cosine();
    //zcos = (1e-15+zcos)/v->dStrands;
    //zsin = (1e-15+zsin)/v->dStrands;
    
    z = z - (z2*z1*z*c + z1*c + z)/(3*z1*c + 2*z*c + 1);

    zd = z-z2;
    temp = sum_sqrs_zd();
  }

  __real__ v->z = __real__ z;
  __imag__ v->z = __imag__ z;

  v->pXTemp[v->Iteration] = __real__ v->z;
  v->pYTemp[v->Iteration] = __imag__ v->z;
}

void Formula_29(CDATA *v)
{
  //MessageBox(NULL, "formula 02", "bug", MB_OK);
  // 4th order Nova
  if (v->jul == 0)
    v->z += 1;
  temp = 1;
  c = v->c;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z = z1 = v->z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    v->z = z-cexp_z()/((z4-1-zsin)/(4*z3)-zcos)+c;

    v->pXTemp[v->Iteration] = __real__ v->z;
    v->pYTemp[v->Iteration] = __imag__ v->z;

    zd = z1-v->z;
    temp = sum_sqrs_zd();
  }
}

void Formula_30(CDATA *v)
{
  //MessageBox(NULL, "formula 02", "bug", MB_OK);
  // 4th order Nova
	
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z -= 1;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    //z = z1 = v->z;
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    //sine_cosine();

    z = z - ((z4+1)/(4*z3)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_31(CDATA *v)
{
  //MessageBox(NULL, "formula 02", "bug", MB_OK);
  // 4th order Nova
	
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z -= 1;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    //z = z1 = v->z;
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    //sine_cosine();

    z = z - ((z4 - z )/(4*z3-1)) + cexp_z() - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_32(CDATA *v)
{
  //MessageBox(NULL, "formula 02", "bug", MB_OK);
  // 4th order Nova
	
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    //z = z1 = v->z;
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z4 - z )/(4*z3-1)) + c + cexp_c();

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_33(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    //z = z1 = v->z;
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z3 - z )/(3*z2-1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_34(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    //z3 = z2*z;
    //z4 = z2*z2;
    //sine_cosine();

    z = z - ((z2 - z)/(2*z - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_35(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z4 - z3 - z2 - z)/(4*z3 - 3*z2 - 2*z - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_36(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    //sine_cosine();

    z = z - ((z4*c - z)/(4*z3*c - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_37(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z += 1;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    //sine_cosine();

    z = z - ((z4 - z)/(4*z3 - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_38(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z4*zsin - z)/(4*z3*zcos - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_39(CDATA *v)
{
	
  /*
  Solve the roots of the function.
  Initialize z to one of the roots.
  Apply Newton's method for solving roots.
  f(z) = z - z/z'

  z   = z*z*z*c - z         ;  // the function
  z'  = 3*z*z*c - 1           ;  // 1st derivative
  z'' = 6*z  - 1                  ;  // 2nd derivative

  6*z   = 1;

  */	
  
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {	
    z1 = z;	
    z2 = z*z;
    z3 = z2*z;

    //z4 = z2*z2;
    sine_cosine();
		
    z = z - ((z3 - z - zsin)/(3*z2 - zcos)) + c;
	   	
    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_40(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z4*c*zsin - z - zsin - z)/(4*z3*c*zcos - 1 - zcos - 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_41(CDATA *v)
{
  z = v->z;
  c = v->c;
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z = z - ((z4*c*zsin + z)/(4*z3*c + zcos)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_42(CDATA *v)
{
  z = v->z;
  //c = v->c;
  c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    //z3 = z2*z;
    //z4 = z2*z2;
    //sine_cosine();

    zsin = z*z + c;

    z = z - ((z2*zsin + z)/(2*z*zsin + 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_43(CDATA *v)
{
  z = v->z;
  //c = v->c;
  c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    //z3 = z2*z;
    //z4 = z2*z2;
    sine_cosine();

    z3 = z*z*z + c;

    z = z - ((z2*z3 + z)/(2*z*z3 + 1)) + c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_44(CDATA *v)
{
  z = v->z;
  c = v->c;
  //c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    //z3 = z2*z;
    //z4 = z2*z2;
    sine_cosine();

    z3 = z*z+c;

    z = z - ((z2*z3 - z)/(2*z*z3 - 1)) - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_45(CDATA *v)
{
  z = v->z;
  c = v->c;
  //c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    //z3 = z2*z;
    //z4 = z2*z2;
    sine_cosine();

    z3 = z*z+c;

    z = z - ((z2*z3 - z - zsin)/(2*z*z3 - 1 - zcos)) - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_46(CDATA *v)
{
  z = v->z;
  c = v->c;
  //c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    //z4 = z2*z2;
    //sine_cosine();

    z4 = z*z*z*z+c;

    z = z - ((z3*z4 - z2 - z)/(3*z2*z4 - 2*z - 1)) - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_47(CDATA *v)
{
  z = v->z;
  c = v->c;
  //c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    //z4 = z2*z2;
    //sine_cosine();

    z4 = 1-z*z+c;

    z = z - ((z3*z4 - z2 - z)/(3*z2*z4 - 2*z - 1)) - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}

void Formula_48(CDATA *v)
{
  z = v->z;
  c = v->c;
  //c = cexp_c();
	
  if (v->jul == 0)
    z = 0;
	
  temp = 1;
  while( v->Iteration++ < v->Maxit && temp > zerotol && temp < dMAX)
  {
    z1 = z;
    z2 = z*z;
    z3 = z2*z;
    z4 = z2*z2;
    sine_cosine();

    z5 = 1+z*z-c;

    z = z - ((z4*z5 - z3 - z2 - z - zsin)/(4*z3*z5 - 3*z2 - 2*z - 1 - zcos)) - c;

    v->pXTemp[v->Iteration] = __real__ z;
    v->pYTemp[v->Iteration] = __imag__ z;

    zd = z1-z;
    temp = sum_sqrs_zd();
  }
}