/*   File: radius.cc

     A program that proves that
     $$
        \forall x\in X \exist y\in Y : |F(x+iy)| > 1,
     $$
     where $F$ is the modulus function described in the article. 
     Here, $X\times Y$ describes a rectangular domain in $\C$.

     Uses: The C-XSC Class Library, available from 
           http://www.uni-karlsruhe.de/~iam/html/language/cxsc/cxsc.html

     Author: Warwick Tucker <warwick@math.uu.se>
     Latest edit: Fri Mar 14 10:15:11 CET 2003
*/

#include <stdlib.h>
#include <iostream>
#include "list.h"              // My own linked template list.
#include "rmath.hpp"           // Include real standard functions.
#include "complex.hpp"         // Include complex package.
#include "interval.hpp"        // Include interval arithmetic package.
#include "imath.hpp"           // Include interval standard functions.
#include "cimath.h"            // Include complex interval standard functions.
#include "cinterval.hpp"       // Include complex interval arithmetic package.
#include "civector.hpp"        // Include complex interval vector package.

using namespace std;
using namespace cxsc;

#define CRECT cinterval // A rectangle in the complex plane.
#define BOX   ivector   // A box in R^n.

////////////////////////////////////////////////////////////////////////////
#undef  POSITIVE_BRANCH
#define POSITIVE_BRANCH
#undef  NEGATIVE_BRANCH
//#define NEGATIVE_BRANCH

////////////////////////////////////////////////////////////////////////////
// Constants and parameters.

const CRECT I(interval(0), interval(1)); // The imaginary unit.

const interval PI         = 4 * atan(interval(1));
const interval PI_HALF    = 2 * atan(interval(1));
const interval THETA      = PI_HALF;
const interval ALPHA      = PI_HALF;
const interval LAMBDA     = interval(101)/10;
const interval k          = interval(1)/10;
const interval R          = 1/k;
const interval a          = LAMBDA - sqr(k);
const interval b          = 1 - LAMBDA;
const interval kb         = k*b;
const interval kab        = k*a*b;
 
//////////////////////////////////////////////////////////////////////////////
// Random functions.

static void InitRandoms  ();
static real RandomReal   (real min, real max);

static void InitRandoms()
{ srand(time(NULL)); }

static real RandomReal(real min, real max)
{ return ((real) rand() / RAND_MAX) * (max - min) + min; }

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

void printConstants()
{
  cout << "Printing all user-defined constants:" << endl
       << "PI      " << PI << endl
       << "PI_HALF " << PI_HALF << endl
       << "THETA   " << THETA << endl
       << "ALPHA   " << ALPHA << endl
       << "LAMBDA  " <<  LAMBDA << endl
       << "k       " << k << endl
       << "R       " << R << endl
       << "a       " << a << endl
       << "b       " << b << endl
       << "kb      " << kb << endl
       << "kab     " << kab << endl;
}

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

CRECT pow(const interval &x, const CRECT &n)
{
  CRECT X = CRECT(x, interval(0));
  return pow(X, n);
}

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

void BCD(CRECT &B, CRECT &C, CRECT &D, const CRECT &z)
{
  CRECT denominator = 2*(a + kb)*cos(PI_HALF*z);
  CRECT localPow = pow(k, -2*z);

  B = (kb + a*localPow)/denominator;
  C = (1 - localPow)/(2*(a + kb)*sin(PI_HALF*z));
  D = (a + kb*localPow)/denominator;
}

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

void MN(CRECT &M, CRECT &N, const CRECT &z)
{
  CRECT Bneg, Cneg, Dneg;
  BCD(Bneg, Cneg, Dneg, -z);

  CRECT Bpos, Cpos, Dpos;
  BCD(Bpos, Cpos, Dpos, +z);

  M = -2*kab*Cneg*Cpos + Bneg*Bpos + Dneg*Dpos;
  N = (Dpos*Bpos + kab*sqr(Cpos)) * (Dneg*Bneg + kab*sqr(Cneg));
}

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

void G(CRECT &Gpos, CRECT &Gneg, const CRECT &z)
{
  CRECT M, N;
  MN(M, N, z); 
  
  Gpos = +(M + sqrt(sqr(M) - 4*N)) / 2; 
  Gneg = +(M - sqrt(sqr(M) - 4*N)) / 2; 
}

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

interval Func(const CRECT &z)
{
  CRECT Gpos, Gneg;
  G(Gpos, Gneg, z);
#ifdef POSITIVE_BRANCH
  return abs(Gpos);
#endif // POSITIVE_BRANCH
#ifdef NEGATIVE_BRANCH
  return abs(Gneg);
#endif // NEGATIVE_BRANCH
}

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

void splitAndStore(const CRECT &localX, List<CRECT> &domainList)
{
  CRECT loX = localX;
  CRECT hiX = localX;

  if ( diam(Re(localX)) > diam(Im(localX)) ) {
    SetRe(loX, interval(Inf(Re(localX)), mid(Re(localX))));
    SetRe(hiX, interval(mid(Re(localX)), Sup(Re(localX))));
  }
  else {
    SetIm(loX, interval(Inf(Im(localX)), mid(Im(localX))));
    SetIm(hiX, interval(mid(Im(localX)), Sup(Im(localX))));
  }
  domainList += loX;
  domainList += hiX;
}

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

void valGreaterThanOne(List<CRECT> &verifiedList,  List<CRECT> &refineList, 
		       List<CRECT> &unverifiedList,List<CRECT> &removedList, 
		       const CRECT &domain, const real &TOL)
{
  refineList += domain;
  while( !IsEmpty(refineList) ) {
    CRECT    localX = First(refineList);
    interval localY = Func(localX);  
    RemoveCurrent(refineList);   
    if ( (Inf(localY) <= 1) && (1 <= Sup(localY)) ) 
      if ( max(diam(Re(localX)), diam(Im(localX))) > TOL )
	splitAndStore(localX, refineList);
      else
	unverifiedList += localX;
    else if ( 1 < Inf(localY) )
      verifiedList += localX;
    else 
      removedList += localX;
  }
}

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

void printList(List<CRECT> &printList, const int &label)
{
  if ( IsEmpty(printList) ) {
    cout << "list # " << label << " is empty." << endl;
    return;
  }
  
  char * txtfile;
  if ( label == 1 )
    txtfile = "verified.txt";
  else if ( label == 2 )
    txtfile = "refine.txt";
  else if ( label == 3 )
    txtfile = "unverified.txt";
  else
    txtfile = "removed.txt";
  
  FILE * TextOutFile = fopen(txtfile,"w");

  cout << "Writing the list to '" << txtfile 
       << "'. (" << Length(printList) << " units)." << endl;
  interval X, Y;
  while( !IsEmpty(printList) ) {
    X = Re(First(printList));
    Y = Im(First(printList));
    RemoveCurrent(printList);
    fprintf(TextOutFile,"%f %f %f %f\n",
	    _double(Inf(X)), _double(Sup(X)), 
	    _double(Inf(Y)), _double(Sup(Y)));
  }
  fflush(TextOutFile);
  fclose(TextOutFile);
}

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

int main(int argc, char * argv[])
{
  CRECT domain = CRECT(interval(1, 4)/4, interval(0, 1));    // Global region.
  //domain = CRECT(interval(10,26)/40, interval(9,11)/100);  // Strip $s_1$.
  //domain = CRECT(interval(1,2)/2, interval(29, 31)/100);   // Strip $s_2$.

  printConstants();
  cout << "domain = " << domain << endl;
  List<CRECT> verifiedList;
  List<CRECT> refineList;
  List<CRECT> unverifiedList;
  List<CRECT> removedList;
  real tolerance = pow(2, -3);
  if ( argc > 1 )
    tolerance = _real(pow(2,atoi(argv[1])));
  cout << "Tolerance      = " << tolerance << endl;
  
#ifdef POSITIVE_BRANCH
  cout << "Computing on the POSITIVE_BRANCH." << endl;
#endif // POSITIVE_BRANCH
#ifdef NEGATIVE_BRANCH
  cout << "Computing on the NEGATIVE_BRANCH." << endl;
#endif // NEGATIVE_BRANCH

  valGreaterThanOne(verifiedList, refineList, unverifiedList, 
		    removedList, domain, tolerance);
  printList(verifiedList, 1);
  printList(refineList, 2);
  printList(unverifiedList, 3);
  printList(removedList, 4);
  
  return 0;
}

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