#define WANT_MATH
#define WANT_STREAM

#include "include.h"
#include "array1.h"
#include "cx.h"

//  Consider the sequence of complex numbers
//     z[0] =0; z[n+1] = z[n] ** 2 + c
//  where ** denotes exponentiation (i.e. square in this case).
//
//  The Mandelbrot set, M, is the set of values of c for which this sequence
//  remains bounded as n tends to infinity.
//
//  This program calculates the set, Z, of values of c for which z[n] = 0
//  for some n > 0.
//
//  Of course, Z is a subset of M. It follows from Montel's theorem in
//  complex analysis that the boundary of M is contained in the closure
//  of Z.
//
//  So if we calculate the points of Z we get an outline of M.
//
//  We can search for the values of Z using contour integration. Think of
//  z[n] as function of c and contour integrate d[n] = z[n]' / z[n]
//  around a square. Here ' denotes derivative wrt c.
//  If there are any zeros of z[n] in the square then the imaginary part of
//  integral will have a positive value depending on the number of zeros;
//  otherwise it will be zero.
//
//  The program uses gaussian integration to evaluate the integral numerically.
//  It uses an iterative scheme to work through values of n and a recursive
//  scheme to search for the zeros.
//
//  Start with a square covering the area we are interested in.
//  If the contour integration program finds it contains a zero divide
//  the square into 4 squares by dividing the horizontal and vertical
//  sides in half. Search each of these squares for zeros. Continue the
//  process until the zeros have been located with the desired accuracy.




#ifdef use_namespace
using namespace RBD_COMMON;             // needed for VC++
using namespace RBD_COMPLEX;
#endif

class SearchClass
{
   int maxit;                    // maximum number of iterations
   Real limit;                   // minimum for deciding a point is a "zero"
   const array1<Real>& gauss;    // integration points
   const array1<Real>& weight;   // integration weights
   int gs;                       // number of integration points
   int g4;                       // gs * 4
   array1<CX> C;                 // for integrating around square
   array1<CX> W;
   array1<CX> Z;
   array1<CX> D;
   array1<bool> skip;
public:
   SearchClass(int mi, Real li, const array1<Real>& g, const array1<Real>& w);
   void IntegrationPoints(CX centre, Real delta);
   int IntegrateSquare(Real delta);
   int Search(CX centre, Real delta, int depth);
};


void GaussianIntegration(array1<Real>& G, array1<Real>& W);


int my_main()
{
   int n = 16;                  // order of integration (16, 8, 4 or 2)
   CX centre(0,0);              // centre of square to be examined
   Real delta = 2.0;            // distance of side of square from centre
   int depth = 9;               // number of binary divisions of side of square
   int maxit = 1000;            // maximum number of iterations

//   int n = 16;
//   CX centre(-0.74916,0.06648);
//   Real delta = 0.00004;
//   int depth = 8;
//   int maxit = 10000;

   array1<Real> G(n);           // these will contain the gaussian integration
   array1<Real> W(n);           // parameters - n will set the number of values.

   GaussianIntegration(G, W);   // get Gaussian integration parameters and load
                                // them into G and W.

   SearchClass sc(maxit, 0.2, G, W); // a class that organises the search

   int s = sc.Search(centre, delta, depth); // search our initial square

   cout << "Number of points = " << s << endl;

   return 0;
}

SearchClass::SearchClass(int mi, Real li, const array1<Real>& g,
   const array1<Real>& w)
      : maxit(mi), limit(li), gauss(g), weight(w)
{
   gs = gauss.size(); g4 = 4 * gs;
   C.resize(g4); W.resize(g4); Z.resize(g4); D.resize(g4); skip.resize(g4);
}


// Setup the integration points and their weights for the square defined
// by centre and delta

void SearchClass::IntegrationPoints(CX centre, Real delta)
{
   int i; int j = 0;
   Imag i_delta = _I_ * delta; CX mid = centre + delta;
   for (i = 0; i < gs; i++)
      { W(j) = weight(i) * i_delta;   C(j) = mid + gauss(i) * i_delta; ++j; }
   mid = centre + i_delta;
   for (i = 0; i < gs; i++)
      { W(j) = - weight(i) * delta;   C(j) = mid - gauss(i) * delta;   ++j; }
   mid = centre - delta;
   for (i = 0; i < gs; i++)
      { W(j) = - weight(i) * i_delta; C(j) = mid - gauss(i) * i_delta; ++j; }
   mid = centre - i_delta;
   for (i = 0; i < gs; i++)
      { W(j) = weight(i) * delta;     C(j) = mid + gauss(i) * delta;   ++j; }
}

// do the integration around the square for each n

int SearchClass::IntegrateSquare(Real delta)
{
   Z = CX(0.0); D = CX(0.0); skip = false;    // set all elements of an array

   for (int n = 1; n <= maxit; n++)
   {
      Real sum = 0.0; bool converge = true;
      CX* z = Z.data(); CX* d = D.data(); CX* c = C.data(); CX* w = W.data();
      bool* s = skip.data();
      for (int j = 0; j < g4; j++)
      {
         if (!*s)
         {
            CX z_next = square(*z) + *c;  *z = z_next;
            if (z_next == 0) return n;  // we have hit an exact zero
            *d *= 2.0;
            CX e = (1.0 - *d * *c) / z_next;
            *d += e;                    // the new value of d
            sum += imag(*w * e);        // integration of imaginary part
                                        // we need do only e since we have
                                        // already shown that the integral of
                                        // the previous d is zero
            Real nz = norm1(z_next);
            if (nz < 4 || delta * ( 1.0 + norm1(*d) ) > 0.0001 * nz)
               converge = false;
            else *s = true;             // this term can't contribute
                                        // significantly in the future
         }
         ++z; ++d; ++c; ++w; ++s;
      }
      if (sum / pi_times_2 > limit) return n;  // have found a zero
      if (converge) return 0;           // no term can contribute in the future
   }

   return -1;                           // exit on iteration limit
                                        // usually means that our square
                                        // is completely in the interior
}

// the recursive search

int SearchClass::Search(CX centre, Real delta, int depth)
{
   IntegrationPoints(centre, delta); int n = IntegrateSquare(delta);
   if (n > 0)
   {
      if (depth > 0)
      {
         Real d2 = delta / 2.0; Imag d2i = d2 * _I_; int k = 0;
         k += Search(centre-d2-d2i, d2, depth-1);
         k += Search(centre-d2+d2i, d2, depth-1);
         k += Search(centre+d2-d2i, d2, depth-1);
         k += Search(centre+d2+d2i, d2, depth-1);
         return k;
      }
      else
      {
         // print coordinates of centre of square and iteration number
         cout << setw(20) << setprecision(15) << centre.real() << " "
              << setw(20) << setprecision(15) << centre.imag() << " "
              << setw(10) << n << endl;
         return 1;
      }
   }
   else return 0;
}

// load the integration parameters into G and W

void GaussianIntegration(array1<Real>& G, array1<Real>& W)
{
   switch (G.size())
   {
   case 16:

      G(8)  = 0.095012509837637;  G(7) = -G(8);
      G(9)  = 0.281603550779259;  G(6) = -G(9);
      G(10) = 0.458016777657227;  G(5) = -G(10);
      G(11) = 0.617876244402644;  G(4) = -G(11);
      G(12) = 0.755404408355003;  G(3) = -G(12);
      G(13) = 0.865631202387832;  G(2) = -G(13);
      G(14) = 0.944575023073233;  G(1) = -G(14);
      G(15) = 0.989400934991650;  G(0) = -G(15);

      W(8)  = 0.189450610455068;  W(7) = W(8);
      W(9)  = 0.182603415044924;  W(6) = W(9);
      W(10) = 0.169156519395003;  W(5) = W(10);
      W(11) = 0.149595988816577;  W(4) = W(11);
      W(12) = 0.124628971255534;  W(3) = W(12);
      W(13) = 0.095158511682493;  W(2) = W(13);
      W(14) = 0.062253523938648;  W(1) = W(14);
      W(15) = 0.027152459411754;  W(0) = W(15);

      break;

   case 8:

      G(4) = 0.183434642495650;  G(3) = -G(4);
      G(5) = 0.525532409916329;  G(2) = -G(5);
      G(6) = 0.796666477413627;  G(1) = -G(6);
      G(7) = 0.960289856497536;  G(0) = -G(7);

      W(4) = 0.362683783378362;  W(3) = W(4);
      W(5) = 0.313706645877887;  W(2) = W(5);
      W(6) = 0.222381034453374;  W(1) = W(6);
      W(7) = 0.101228536290376;  W(0) = W(7);

      break;

   case 4:

      G(2) = 0.339981043584856;  G(1) = -G(2);
      G(3) = 0.861136311594053;  G(0) = -G(3);

      W(2) = 0.652145154862546;  W(1) = W(2);
      W(3) = 0.347854845137454;  W(0) = W(3);

      break;

   case 2:

      G(1) = 0.577350269189626;  G(0) = -G(1);

      W(1) = 1.0; W(0) = 1.0;

      break;

   default: cout << "invalid parameter" << endl; exit(1);

   }
}

// call my_main() - use this to catch exceptions
// use macros for exception names for compatibility with simulated exceptions
int main()
{
   Tracer tr("main ");
   Try  { return my_main(); }
   Catch(BaseException) { cout << BaseException::what() << "\n"; }
   CatchAll { cout << "\nProgram fails - exception generated\n\n"; }
   return 0;
}