// Samples from -RANGE to +RANGE
#define REAL_RANGE 8.0
#define IMAG_RANGE 8.0

// Actually performs steps 2*STEP+1
#define REAL_STEPS 512
#define IMAG_STEPS 512

// Numer of times to evaluate the grid when measuring timing
#define TIMING_LOOPS 10000

#include <stdio.h>

#include <gmp.h>
#include <mpfr.h>

#include <complex.h>
#include <math.h>
#include <float.h>

#include <time.h>
#include <sys/times.h>
#include <unistd.h>

#define MP_PRECISION 500

mpfr_t a, b, den, tmp1, tmp2 ;
mpfr_t rp1, ip1, sinha, cosha, sinb, cosb ;
mpfr_t rp2, ip2, sinh2a, cosh2a, sin2b, cos2b ;

complex double
do_ctanh_1 ( complex double z )
{
  /* 
                   sinh(2*a) + i*sin(2*b)
     ctanh(a+bi) = ----------------------
                   cosh(2*a) + cos(2*b)
  */

  /* z = a + bi */
  mpfr_set_d (a, 2.0*creal(z), GMP_RNDN);
  mpfr_set_d (b, 2.0*cimag(z), GMP_RNDN);
	  
  mpfr_sinh_cosh (sinh2a, cosh2a, a, GMP_RNDN);
  mpfr_sin_cos (sin2b, cos2b, b, GMP_RNDN);
	  
  /* Denominator */
  mpfr_add (den, cosh2a, cos2b, GMP_RNDN);
	  
  /* Real component */
  mpfr_div (rp1, sinh2a, den, GMP_RNDN);
  
  /* Imaginary component */
  mpfr_div (ip1, sin2b, den, GMP_RNDN);

  return mpfr_get_d (rp1, GMP_RNDN) + mpfr_get_d (ip1, GMP_RNDN)*I ;
}

complex double
do_ctanh_2 ( complex double z )
{
  /* 
                   sinh(a)*cosh(a) + i*sin(b)*cos(b)
     ctanh(a+bi) = ---------------------------------
                          cos(b)^2 + sinh(a)^2
  */

  /* z = a + bi */
  mpfr_set_d (a, creal(z), GMP_RNDN);
  mpfr_set_d (b, cimag(z), GMP_RNDN);
	  
  mpfr_sinh_cosh (sinha, cosha, a, GMP_RNDN);
  mpfr_sin_cos (sinb, cosb, b, GMP_RNDN);
	  
  /* Denominator */
  mpfr_mul (tmp1, cosb, cosb, GMP_RNDN);
  mpfr_mul (tmp2, sinha, sinha, GMP_RNDN);
  mpfr_add (den, tmp1, tmp2, GMP_RNDN);
	  
  /* Real component */
  mpfr_mul (tmp1, sinha, cosha, GMP_RNDN);
  mpfr_div (rp2, tmp1, den, GMP_RNDN);
  
  /* Imaginary component */
  mpfr_mul (tmp1, sinb, cosb, GMP_RNDN);
  mpfr_div (ip2, tmp1, den, GMP_RNDN);

  return mpfr_get_d (rp2, GMP_RNDN) + mpfr_get_d (ip2, GMP_RNDN)*I ;
}

/* function x = compare_double(a,b,del) */
/*   A = min(a,b); */
/*   B = max(a,b); */
/*   A = max(A.*(1+del),A.*(1-del)); */
/*   B = min(B.*(1+del),B.*(1-del)); */
/*   x = (A <= B) ; */
/*   y = ((a==0) & (b==0)) ; */
/*   x(y) = 0 ; */
/* endfunction */

int
doubles_differ ( double a, double b, double tol )
{
  double A = fmin (a, b) ;
  double B = fmax (a, b) ;
  A = fmax (A*(1+tol), A*(1-tol)) ;
  B = fmin (B*(1+tol), B*(1-tol)) ;

  if (a==0.0 & b==0.0)
    return 0 ;
  else
    return (A <= B) ;
}

/* function [x,y] = compare_complex(a,b,del) */
/*   x = compare_double(real(a), real(b),del); */
/*   y = compare_double(imag(a), imag(b),del); */
/* endfunction */

int
complex_differ ( double complex a, double complex b, double tol )
{
  int retval = 0 ;
  if (doubles_differ (creal(a), creal(b), tol) > 0)
    retval += 1 ;
  if (doubles_differ (cimag(a), cimag(b), tol) > 0)
    retval += 2 ;
  // if (retval > 0)
  //   printf("Differ! %g %g | %g %g\n", creal(a), cimag(a), creal(b), cimag(b)) ;
  return retval ;
}

int main (void)
{
  mpfr_init2 (a, MP_PRECISION);
  mpfr_init2 (b, MP_PRECISION);
  mpfr_init2 (den, MP_PRECISION);
  mpfr_init2 (tmp1, MP_PRECISION);
  mpfr_init2 (tmp2, MP_PRECISION);
  mpfr_init2 (rp1, MP_PRECISION);
  mpfr_init2 (ip1, MP_PRECISION);
  mpfr_init2 (sinh2a, MP_PRECISION);
  mpfr_init2 (cosh2a, MP_PRECISION);
  mpfr_init2 (sin2b, MP_PRECISION);
  mpfr_init2 (cos2b, MP_PRECISION);
  mpfr_init2 (rp2, MP_PRECISION);
  mpfr_init2 (ip2, MP_PRECISION);
  mpfr_init2 (sinha, MP_PRECISION);
  mpfr_init2 (cosha, MP_PRECISION);
  mpfr_init2 (sinb, MP_PRECISION);
  mpfr_init2 (cosb, MP_PRECISION);

  int delta[10][4] = {0} ; 
  for ( int i=0; i<10; i++ )
    for ( int j=1; j<4; j++ )
      delta[i][j] = 0 ;

  int internal_errors = 0 ;

  for ( int i=-REAL_STEPS; i<=REAL_STEPS; i++)
    {
      for ( int j=-IMAG_STEPS; j<=IMAG_STEPS; j++)
	{
	  double di = REAL_RANGE/REAL_STEPS*i ;
	  double dj = IMAG_RANGE/IMAG_STEPS*j ;
	  
	  complex double z = di + dj*I;
	  complex double t = ctanh(z);
	  complex double r1 = do_ctanh_1 (z) ;
	  complex double r2 = do_ctanh_2 (z) ;

	  if ( complex_differ(r1,r2,DBL_EPSILON) > 0 )
	    {
	      internal_errors++ ;
	      fprintf(stderr, "\nInternal inconsistency: %g %g | %g %g | %g %g \n", 
		      creal(z), cimag(z), creal(r1), cimag(r1), creal(r2), cimag(r2)) ;
	      fflush(stderr) ;
	    }
	  for ( int k=0; k<10; k++ )
	    {
	      int d = complex_differ(r1,t,(k+1)*DBL_EPSILON) ;
	      delta[k][d]++ ;
	    }
	}
      printf(".", i) ;
      fflush(stdout) ;
    }

  printf("# Timing\n") ;
  struct tms start ;
  struct tms end ;
  times (&start);
  for ( int l=0 ; l<TIMING_LOOPS ; l++ )
    for ( int i=-REAL_STEPS; i<=REAL_STEPS; i++)
      {
	double r = REAL_RANGE/REAL_STEPS*i ;
	for ( int j=-IMAG_STEPS; j<=IMAG_STEPS; j++)
	  ctanh(r + IMAG_RANGE/IMAG_STEPS*j*I);
      }
  times (&end);
  printf("\n\n") ;

  printf(" Error magn  Real part  Imag part  Both\n") ;
  printf(" ----------- ---------- ---------- ---------\n") ;
  for ( int i=0; i<10; i++ )
    {
      printf("  >=%2d eps  ", i+1) ;
      for ( int j=1; j<4; j++ )
	printf("%10d ", delta[i][j]) ;
      printf("\n") ;
    }
  printf("\n Grid size:                %d x %d\n", 2*REAL_STEPS+1, 2*IMAG_STEPS+1) ;
  printf(" Top Left corner:           %+f + %+f i\n", -REAL_RANGE, IMAG_RANGE) ;
  printf(" Bottom Right corner:       %+f + %+f i\n",	REAL_RANGE, -IMAG_RANGE) ;  
  printf(" Time to evaluate grid:     %f (sec) [%d loops]\n",	((double)(end.tms_utime - start.tms_utime))/sysconf(_SC_CLK_TCK), TIMING_LOOPS ) ;
  printf(" Internal inconsistencies:  %d\n\n", internal_errors) ;

  return 0;
}