/* 
* Copyright (C) 2001-2009, CompHEP Collaboration
* Copyright (C) 1997, Dmitry Kovalenko 
* ------------------------------------------------------
*/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>

#include "service2/include/chep_limits.h"
#include "chep_crt/include/chep_crt.h"

#include "kinaux.h"
#include "regul.h"
#include "regfunal.h"

static double 
sing1 (int ityp, double pole, double width, int deg, double x)
{
  double dx = x - pole;
  double res = 0.;

  if (ityp == 1)
    {
      if (deg == 2) res = 1. / (dx * dx);
      else res = 1. / fabs (dx);
    }
  else
    {
      if (width == 0.) {
        if (deg == 2) res = 1. / (dx * dx);
        else res = 1. / fabs (dx);
      }
      return 1 / (dx * dx + width * width);
    }
  return res;
}

static double 
sing2 (int ityp, double pole, double width, int deg, double x)
{
  double res = 0.;
  double dx = x - pole;

 if (ityp != 1 && ityp != 2) {
    printf ("comphep (error): sing2, invalid ityp = %d. ABORT\n", ityp);
    exit (99);
  }

  if (ityp == 1) {
    if (deg == 2) {
      res = -1 / dx;
    } else {
      if (dx > 0)
        res = log (dx);
      else
        res = -log (-dx);
    }
  } else {
    if (width == 0.) {
        if (deg == 2) {
          res = -1. / dx;
        } else {
          if (dx > 0)
            res = log (dx);
          else
            res = -log (-dx);
        }
    } else {
      res = atan (dx / width) / width;
    }
  }

  return res;
}

static double 
sing3 (int ityp, double pole, double width, int deg, double x, double xmin)
{
  double res = pole;

  if (ityp == 1) {
      if (deg == 2) {
        res -= 1. / x;
      } else {
        if (pole < xmin)
          res += exp (x);
        else
          res -= exp (-x);
      }
  } else {
      if (width == 0.) {
        if (deg == 2) {
          res -= 1. / x;
        } else {
          if (pole < xmin)
            res += exp (x);
          else
            res -= exp (-x);
        }
      } else {
        res += width * tan (x * width);
      }
  }

  return res;
}


static double
regfun_0_ (int factOnly, int itype, int nsing,
	   sing_struct * singar, double xmin, double xmax, double xx, double *xout)
{
  int tint0;
  int i, j, ityp;
  double delt0 = 1e-13;
  double ai[100];
  double bi[100];
  double rnorm_0;
  double delt;
  double tintxx;
  double factor;

  if (nsing == 0) {
    if (!factOnly) {
      *xout = xmin * (1 - xx) + xx * (xmax);
      if (*xout > xmax)
        *xout = xmax;
      if (*xout < xmin)
        *xout = xmin;
    }
    return xmax - xmin;
  }

  /* Parameter adjustments */

  if (itype > 0)   /* constant is included */
    {
      ityp = itype;
      tint0 = 1;
      rnorm_0 = 1 / (xmax - xmin);
    }
  else
    {
      ityp = -itype;
      tint0 = 0;
      rnorm_0 = 0.;
    }

  for (i = 0; i < nsing; ++i) {
    if (ityp == 1 || (ityp == 2 && singar[i].width == 0.))
      {
        if (singar[i].pos >= xmin && singar[i].pos <= xmax)
          {
            delt = delt0 * (fabs (xmin) + fabs (xmax));
            if (singar[i].pos - delt <= xmin)
              {
        	singar[i].pos = xmin - delt;
              }
            else if (singar[i].pos + delt >= xmax)
              {
        	singar[i].pos = xmax + delt;
              }
            else
              {
        	printf ("%f < x < %f\n", xmin, xmax);
        	printf ("Type= %d  Bad reg= %d\n", itype, i + 1);

        	for (j = 0; j < nsing; ++j) {
        	  printf ("%d %f %f %d\n", j + 1, singar[j].pos, singar[j].width, singar[j].power);
        	}
        	printf ("Error in %i-th regularization \n", i + 1);
        	exit (0);
              }
          }
      }
  }

  for (i = 0; i < nsing; ++i) {
    ai[i] = sing2 (ityp, singar[i].pos, singar[i].width, singar[i].power, xmin);
    bi[i] = sing2 (ityp, singar[i].pos, singar[i].width, singar[i].power, xmax);
  }

  if (!factOnly) {
    tintxx = (tint0 + nsing) * xx;
    if (tint0 > tintxx) {
      *xout = xmin + tintxx / rnorm_0;
    } else {
      tintxx -= tint0;
      j = tintxx;
      tintxx -= j;
      *xout = sing3 (ityp, singar[j].pos, singar[j].width, singar[j].power, tintxx * (bi[j] - ai[j]) + ai[j], xmin);
    }

    if (*xout > xmax)
      *xout = xmax;
    if (*xout < xmin)
      *xout = xmin;
  }

  factor = rnorm_0;
  for (i = 0; i < nsing; ++i) {
    factor += sing1 (ityp, singar[i].pos, singar[i].width, singar[i].power, *xout) / (bi[i] - ai[i]);
  }

  return (tint0 + nsing) / factor;
}


double 
regfun (int itype, int nsing, sing_struct * singar, double xmin, double xmax, double xx, double *xout)
{
  double fct = regfun_0_ (0, itype, nsing, singar, xmin, xmax, xx, xout);
  return fct;
}


double 
regfct (int itype, int nsing, sing_struct * singar, double xmin, double xmax, double xout)
{
  double fct = regfun_0_ (1, itype, nsing, singar, xmin, xmax, 0., &xout);
  return fct;
}