/****************************************************************************/
/*                                                                          */
/*        COUNTIT - Iteration for convergence test                          */
/*                                                                          */
/****************************************************************************/

#define STRICT
#include <windows.h>
#include <commctrl.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>
#include <new.h>
#include <float.h>

#include "const.h"
#include "classes.h"
#include "extern.h"
#include "funcdefs.h"
#include "julia.rh"

// These define the loop test region; they are not critical and may be changed
#define LOOP_TEST_START 50
#define LOOP_TEST_END 2000

/*
   For Julia sets, c defines the particular set and z_start is the
   location on the plane being evaluated.
   For Mandelbrot, c is the location in the plane being evaluated,
   and the starting Z for iteration is 0.  But when this is called
   for Mandelbrot, we might as well start the iteration at c
   because if we start at 0 the next Z will be c.
*/

unsigned int countit (
   COMPLEX *c ,        // Constant in iteration Z(k+1) = Z(k)^2 + c
   COMPLEX *z_start ,  // Starting point of iteration; evaluated point in plane
   COMPLEX *end ,      // ending point; needed only for BinaryDecomp
   double *loop_real , // Work area DwellLimit long
   double *loop_imag   // Ditto
   )
{
   unsigned int count ;
   double z_real, z_imag, t_real, t_imag, rsq, isq, *lp_real, *lp_imag, escape ;

   z_real = z_start->real ;
   z_imag = z_start->imag ;

   escape = EscapeConstant * EscapeConstant ;

   count = 0 ;
   lp_real = loop_real ;
   lp_imag = loop_imag ;

   for (;;) {
      rsq = z_real * z_real ;
      isq = z_imag * z_imag ;

      if ((rsq + isq) > escape) {
#if MAX_THREADS == 1
         ++out_by_escape ;
#endif
         break ;
         }

      if (++count == DwellLimit) { // Breaks one too soon for Julia; no big deal
#if MAX_THREADS == 1
         ++out_by_max ;
#endif
         break ;
         }

      t_real = rsq - isq  +  c->real ;
      t_imag = 2. * z_real * z_imag  +  c->imag ;

      //  Every 50 iterations check for convergence to fixed point

      if ((count % 50 == 0)
       && (fabs (z_real - t_real) * 1.e14 < fabs ( z_real ))
       && (fabs (z_imag - t_imag) * 1.e14 < fabs ( z_imag ))) {
#if MAX_THREADS == 1
         ++out_by_conv ;
#endif
         count = DwellLimit ;
         break ;
         }

      //  If count gets to LOOP starting limit, commence loop test

      if (count >= LOOP_TEST_START  &&  count <= LOOP_TEST_END) {
         if ((count > LOOP_TEST_START)  // Ensure first element is defined
          && (fabs ( t_real - loop_real[(count-LOOP_TEST_START)/2] ) * 1.e14 <= fabs ( t_real ))
          && (fabs ( t_imag - loop_imag[(count-LOOP_TEST_START)/2] ) * 1.e14 <= fabs ( t_imag ))) {
            count = DwellLimit ;
#if MAX_THREADS == 1
            ++out_by_loop ;
#endif
            break ;
            }
         *lp_real++ = t_real ;
         *lp_imag++ = t_imag ;
         }

      z_real = t_real ;
      z_imag = t_imag ;
      } // Endless iteration loop

   end -> real = z_real ;
   end -> imag = z_imag ;
   return count ;
}
