/* ------------------------------------------------------------------------
   COMPLEX NUMBERS (2D-NUMBERS)
   ------------------------------------------------------------------------ */
class complexNumber
{
protected:

  double r;  // real part
  double i;  // imaginary part

public:

  // create 0 / equal to other / with real & imaginary part
  complexNumber (void) { r=0; i=0; }
  complexNumber (const class complexNumber &c) { r=c.r; i=c.i; }
  complexNumber (double realPart, double imagPart) { r=realPart; i=imagPart; }

  // set number (real & imaginary part)
  void set (double realPart, double imagPart) { r=realPart; i=imagPart; }

  // get real / imaginary part
  double real (void) const { return (r); }
  double imag (void) const { return (i); }

  // euclidean dist to origin / other complex Number
  double origin (void) const { return ( sqrt(r*r + i*i) ); }
  double dist (const class complexNumber &c) class { 
    return ( sqrt( (r-c.r)*(r-c.r) + (i-c.i)*(i-c.i)) ); }

  // make equal to / add to this / subtract from this / multiplicate by
  class complexNumber operator =  (const class complexNumber &c) { 
    r = c.r; 
    i = c.i; 
    return(*this);     }
  class complexNumber operator += (const class complexNumber &c) { 
    r += c.r; 
    i += c.i; 
    return(*this);     }
  class complexNumber operator -= (const class complexNumber &c) { 
    r -= c.r; 
    i -= c.i; 
    return(*this);     }
  class complexNumber operator *= (const class complexNumber &c) { 
    double R=r, I=i;
    r = (R*c.r)-(I*c.i); 
    i = (R*c.i)+(I*c.r);
    return (*this);     }

  // add / subtract / multiplicate two numbers
  class complexNumber operator + (const class complexNumber &c) const {
    class complexNumber s (*this); s += c; return (s); }
  
  class complexNumber operator - (const class complexNumber &c) const {
    class complexNumber s (*this); s -= c; return (s); }
  
  class complexNumber operator * (const class complexNumber &c) const {
    class complexNumber s (*this); s *= c; return (s); }

  // compare two numbers
  bool operator == (const class complexNumber &c) const { 
    return ((real()==c.real()) && (imag()==c.imag())); }
};


/* ------------------------------------------------------------------------
   QUATERNION NUMBERS (4D-NUMBERS) 
   ------------------------------------------------------------------------ */
class quaternion 
{
protected:

  double x;
  double y;
  double z;
  double w;

public:

  // create 0 / equal to other / with 4 parts a,b,c,d
  quaternion (void) { x=0; y=0; z=0; w=0; }
  quaternion (const class quaternion &q) { x=q.x; y=q.y; z=q.z; w=q.w; }
  quaternion (double a, double b, double c, double d) { x=a; y=b; z=c; w=d; }

  // set 4 parts 
  void set (double a, double b, double c, double d) { x=a; y=b; z=c; w=d; }

  // get 4 parts
  double first  (void) const { return (x); }
  double second (void) const { return (y); }
  double third  (void) const { return (z); }
  double fourth (void) const { return (w); }

  // euclidean distances to origin / other quaternion
  double origin (void) const { return (sqrt(x*x + y*y + z*z + w*w)); }
  double dist (const class quaternion &q) const { 
    return ( sqrt((x-q.x)*(x-q.x) + 
		  (y-q.y)*(y-q.y) +
		  (z-q.z)*(z-q.z) +
		  (w-q.w)*(w-q.w)) ); }

  // make equal to / add to this / subtract from this / multiplicate by
  quaternion operator = (const class quaternion &q) { 
    x = q.x; 
    y = q.y; 
    z = q.z; 
    w = q.w; 
    return (*this); }
  quaternion operator += (const class quaternion &q) { 
    x += q.x;
    y += q.y;
    z += q.z;
    w += q.w;
    return (*this); }
  quaternion operator -= (const class quaternion &q) { 
    x -= q.x;
    y -= q.y;
    z -= q.z;
    w -= q.w;
    return (*this); }
  quaternion operator *= (const class quaternion &q) { 
    double X=x, Y=y, Z=z, W=w;
    x = X*q.x - Y*q.y - Z*q.z - W*q.w;
    y = X*q.y + Y*q.x + Z*q.w - W*q.z;
    z = X*q.z + Z*q.x + W*q.y - Y*q.w;
    w = X*q.w + W*q.x + Y*q.z - Z*q.y;
    return (*this); }

  // add / subtract / multiplicate two numbers
  class quaternion operator + (const class quaternion &c) const {
    class quaternion s (*this); s += c; return (s); }
  class quaternion operator - (const class quaternion &c) const {
    class quaternion s (*this); s -= c; return (s); }
  class quaternion operator * (const class quaternion &c) const {
    class quaternion s (*this); s *= c; return (s); }

  // compare two numbers
  bool operator == (const class quaternion &c) const { 
    return ((first()==c.first()) && (second()==c.second()) && 
	    (third()==c.third()) && (fourth()==c.fourth())); }

/* ------------------------------------------------------------------------
   MANDELBROT FRACTAL 
   ------------------------------------------------------------------------ */
template <class numbertype>
class mandelbrot
{
protected:

  double bout;   // bailout - stop if dist of z to origin greater bout
  long   maxit;  // max iteration - stop if iterations >= maxit

public:

  // create default / with special bailout and max. iterations
  mandelbrot (void) { bout = 2.0; maxit = 50; }
  mandelbrot (double bailout, long maxiter) { bout = bailout; maxit = maxiter;}

  // set bailout / max. iterations 
  void set (double bailout, long maxiter) { bout = bailout; maxit = maxiter; }

  // get bailout / max. iterations 
  double bailout (void) const { return (bout);  }
  long   maxiter (void) const { return (maxit); }

  // eval point c (return number of steps)
  virtual long eval (const numbertype &c) const { 
    numbertype z; return (eval(c,z,0)); }
 
  // is point c part of mandelbrot set (eval >= max.iterations)
  virtual bool partof (const numbertype &c) const { 
    return ((eval(c)>=maxit)); }
    
protected:

  virtual long eval (const numbertype &c, const numbertype &z, long n) const {
    numbertype zi (z);
    for ( ; n<maxit; zi = (zi*zi)+c, n++) 
      if (zi.origin()>=bout) return (n);
    return (maxit); }

  // rucursive version
  // virtual long eval (const numbertype &c, numbertype z, long n) {
  //   if (n>=maxit || z.origin()>=bout) return(n);
  //   else return( eval (c, (z*z)+c, n+1) ); } 

};


/* ------------------------------------------------------------------------
   JULIA FRACTAL 
   ------------------------------------------------------------------------ */
template <class numbertype>
class julia : public mandelbrot<numbertype>
{
protected:
public:

  // create default / with special bailout and max. iterations
  julia (void) { bout = 2.0; maxit = 50; }
  julia (double bailout, long maxiter) { bout = bailout; maxit = maxiter; }

  // eval point z of mandelbrot point c (return number of steps)
  virtual long eval (const numbertype &c, const numbertype &z) const {
    return (mandelbrot<numbertype>::eval(c,z,0)); }

  // is point z part of mandelbrot set (eval >= max.iterations)
  virtual int partof (const numbertype &c, const numbertype &z) const { 
    return ((eval(c,z)>=maxit)); }
};