/*

Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
              2005, 2006, 2007, 2008, 2009 John W. Eaton
Copyright (C) 2010 VZLU Prague

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
<http://www.gnu.org/licenses/>.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cfloat>

#include "lo-error.h"
#include "lo-ieee.h"
#include "lo-mappers.h"
#include "lo-math.h"
#include "lo-specfun.h"
#include "lo-utils.h"
#include "oct-cmplx.h"

#include "f77-fcn.h"

// double -> double mappers.

double
arg (double x)
{
  return atan2 (0.0, x);
}

double
conj (double x)
{
  return x;
}

double
fix (double x)
{
  return gnulib::trunc (x);
}

double
imag (double)
{
  return 0.0;
}

double
real (double x)
{
  return x;
}

double
xround (double x)
{
  return gnulib::round (x);
}

double
xtrunc (double x)
{
  return gnulib::trunc (x);
}

double 
xroundb (double x)
{
  double t = xround (x);

  if (fabs (x - t) == 0.5)
    t = 2 * xtrunc (0.5 * t);

  return t;
}

double
signum (double x)
{
  double tmp = 0.0;

  if (x < 0.0)
    tmp = -1.0;
  else if (x > 0.0)
    tmp = 1.0;

  return xisnan (x) ? octave_NaN : tmp;
}

double
mod (double x, double y)
{
  if (y == 0)
    return x;

  double r = fmod (x, y);
  return ((r < 0) != (y < 0)) ? y+r : r;
}

double
xlog2 (double x)
{
#if defined (HAVE_LOG2)
  return log2 (x);
#else
#if defined (M_LN2)
  static double ln2 = M_LN2;
#else
  static double ln2 = log (2);
#endif

  return log (x) / ln2;
#endif
}

Complex
xlog2 (const Complex& x)
{
#if defined (M_LN2)
  static double ln2 = M_LN2;
#else
  static double ln2 = log (2);
#endif

  return std::log (x) / ln2;
}

double
xexp2 (double x)
{
#if defined (HAVE_EXP2)
  return exp2 (x);
#else
#if defined (M_LN2)
  static double ln2 = M_LN2;
#else
  static double ln2 = log (2);
#endif

  return exp (x * ln2);
#endif
}

double
xlog2 (double x, int& exp)
{
  return frexp (x, &exp);
}

Complex
xlog2 (const Complex& x, int& exp)
{
  double ax = std::abs (x);
  double lax = xlog2 (ax, exp);
  return (ax != lax) ? (x / ax) * lax : x;
}

// double -> bool mappers.

#if ! defined(HAVE_CMATH_ISNAN)
bool
xisnan (double x)
{
  return lo_ieee_isnan (x);
}
#endif

#if ! defined(HAVE_CMATH_ISFINITE)
bool
xfinite (double x)
{
  return lo_ieee_finite (x);
}
#endif

#if ! defined(HAVE_CMATH_ISINF)
bool
xisinf (double x)
{
  return lo_ieee_isinf (x);
}
#endif

bool
octave_is_NA (double x)
{
  return lo_ieee_is_NA (x);
}

bool
octave_is_NaN_or_NA (double x)
{
  return lo_ieee_isnan (x);
}

// (double, double) -> double mappers.

// complex -> complex mappers.

Complex
acos (const Complex& x)
{
  static Complex i (0, 1);

  return -i * (log (x + i * (sqrt (1.0 - x*x))));
}

Complex
acosh (const Complex& x)
{
  return log (x + sqrt (x*x - 1.0));
}

Complex
asin (const Complex& x)
{
  static Complex i (0, 1);

  return -i * log (i*x + sqrt (1.0 - x*x));
}

Complex
asinh (const Complex& x)
{
  return log (x + sqrt (x*x + 1.0));
}

Complex
atan (const Complex& x)
{
  static Complex i (0, 1);

  return i * log ((i + x) / (i - x)) / 2.0;
}

Complex
atanh (const Complex& x)
{
  return log ((1.0 + x) / (1.0 - x)) / 2.0;
}

Complex
ceil (const Complex& x)
{
  return Complex (ceil (real (x)), ceil (imag (x)));
}

Complex
fix (const Complex& x)
{
  return Complex (fix (real (x)), fix (imag (x)));
}

Complex
floor (const Complex& x)
{
  return Complex (floor (real (x)), floor (imag (x)));
}

Complex
xround (const Complex& x)
{
  return Complex (xround (real (x)), xround (imag (x)));
}

Complex
xroundb (const Complex& x)
{
  return Complex (xroundb (real (x)), xroundb (imag (x)));
}

Complex
signum (const Complex& x)
{
  double tmp = abs (x);

  return tmp == 0 ? 0.0 : x / tmp;
}

// complex -> bool mappers.

bool
octave_is_NA (const Complex& x)
{
  return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
}

bool
octave_is_NaN_or_NA (const Complex& x)
{
  return (xisnan (real (x)) || xisnan (imag (x)));
}

// (complex, complex) -> complex mappers.

// FIXME -- need to handle NA too?

Complex
xmin (const Complex& x, const Complex& y)
{
  return abs (x) <= abs (y) ? x : (xisnan (x) ? x : y);
}

Complex
xmax (const Complex& x, const Complex& y)
{
  return abs (x) >= abs (y) ? x : (xisnan (x) ? x : y);
}


// float -> float mappers.

float
arg (float x)
{
  return atan2 (0.0f, x);
}

float
conj (float x)
{
  return x;
}

float
fix (float x)
{
  return gnulib::truncf (x);
}

float
imag (float)
{
  return 0.0;
}

float
real (float x)
{
  return x;
}

float
xround (float x)
{
  return gnulib::round (x);
}

float
xtrunc (float x)
{
  return gnulib::truncf (x);
}

float 
xroundb (float x)
{
  float t = xround (x);

  if (fabs (x - t) == 0.5)
    t = 2 * xtrunc (0.5 * t);

  return t;
}

float
signum (float x)
{
  float tmp = 0.0;

  if (x < 0.0)
    tmp = -1.0;
  else if (x > 0.0)
    tmp = 1.0;

  return xisnan (x) ? octave_Float_NaN : tmp;
}

float
mod (float x, float y)
{
  if (y == 0)
    return x;

  float r = fmodf (x, y);
  return ((r < 0) != (y < 0)) ? y+r : r;
}

float
xlog2 (float x)
{
#if defined (HAVE_LOG2)
  return log2 (x);
#else
#if defined (M_LN2)
  static float ln2 = M_LN2;
#else
  static float ln2 = log2 (2);
#endif

  return log (x) / ln2;
#endif
}

FloatComplex
xlog2 (const FloatComplex& x)
{
#if defined (M_LN2)
  static float ln2 = M_LN2;
#else
  static float ln2 = log (2);
#endif

  return std::log (x) / ln2;
}

float
xexp2 (float x)
{
#if defined (HAVE_EXP2)
  return exp2 (x);
#else
#if defined (M_LN2)
  static float ln2 = M_LN2;
#else
  static float ln2 = log2 (2);
#endif

  return exp (x * ln2);
#endif
}

float
xlog2 (float x, int& exp)
{
  return frexpf (x, &exp);
}

FloatComplex
xlog2 (const FloatComplex& x, int& exp)
{
  float ax = std::abs (x);
  float lax = xlog2 (ax, exp);
  return (ax != lax) ? (x / ax) * lax : x;
}

// float -> bool mappers.

#if ! defined(HAVE_CMATH_ISNANF)
bool
xisnan (float x)
{
  return lo_ieee_isnan (x);
}
#endif

#if ! defined(HAVE_CMATH_ISFINITEF)
bool
xfinite (float x)
{
  return lo_ieee_finite (x);
}
#endif

#if ! defined(HAVE_CMATH_ISINFF)
bool
xisinf (float x)
{
  return lo_ieee_isinf (x);
}
#endif

bool
octave_is_NA (float x)
{
  return lo_ieee_is_NA (x);
}

bool
octave_is_NaN_or_NA (float x)
{
  return lo_ieee_isnan (x);
}

// (float, float) -> float mappers.

// complex -> complex mappers.

FloatComplex
acos (const FloatComplex& x)
{
  static FloatComplex i (0, 1);

  return -i * (log (x + i * (sqrt (static_cast<float>(1.0) - x*x))));
}

FloatComplex
acosh (const FloatComplex& x)
{
  return log (x + sqrt (x*x - static_cast<float>(1.0)));
}

FloatComplex
asin (const FloatComplex& x)
{
  static FloatComplex i (0, 1);

  return -i * log (i*x + sqrt (static_cast<float>(1.0) - x*x));
}

FloatComplex
asinh (const FloatComplex& x)
{
  return log (x + sqrt (x*x + static_cast<float>(1.0)));
}

FloatComplex
atan (const FloatComplex& x)
{
  static FloatComplex i (0, 1);

  return i * log ((i + x) / (i - x)) / static_cast<float>(2.0);
}

FloatComplex
atanh (const FloatComplex& x)
{
  return log ((static_cast<float>(1.0) + x) / (static_cast<float>(1.0) - x)) / static_cast<float>(2.0);
}

FloatComplex
ceil (const FloatComplex& x)
{
  return FloatComplex (ceil (real (x)), ceil (imag (x)));
}

FloatComplex
fix (const FloatComplex& x)
{
  return FloatComplex (fix (real (x)), fix (imag (x)));
}

FloatComplex
floor (const FloatComplex& x)
{
  return FloatComplex (floor (real (x)), floor (imag (x)));
}

FloatComplex
xround (const FloatComplex& x)
{
  return FloatComplex (xround (real (x)), xround (imag (x)));
}

FloatComplex
xroundb (const FloatComplex& x)
{
  return FloatComplex (xroundb (real (x)), xroundb (imag (x)));
}

FloatComplex
signum (const FloatComplex& x)
{
  float tmp = abs (x);

  return tmp == 0 ? 0.0 : x / tmp;
}

// complex -> bool mappers.

bool
octave_is_NA (const FloatComplex& x)
{
  return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
}

bool
octave_is_NaN_or_NA (const FloatComplex& x)
{
  return (xisnan (real (x)) || xisnan (imag (x)));
}

// (complex, complex) -> complex mappers.

// FIXME -- need to handle NA too?

FloatComplex
xmin (const FloatComplex& x, const FloatComplex& y)
{
  return abs (x) <= abs (y) ? x : (xisnan (x) ? x : y);
}

FloatComplex
xmax (const FloatComplex& x, const FloatComplex& y)
{
  return abs (x) >= abs (y) ? x : (xisnan (x) ? x : y);
}

Complex
rc_acos (double x)
{
  return fabs (x) > 1.0 ? acos (Complex (x)) : Complex (acos (x));
}

FloatComplex
rc_acos (float x)
{
  return fabsf (x) > 1.0f ? acos (FloatComplex (x)) : FloatComplex (acosf (x));
}

Complex rc_acosh (double x)
{
  return x < 1.0 ? acosh (Complex (x)) : Complex (acosh (x));
}

FloatComplex rc_acosh (float x)
{
  return x < 1.0f ? acosh (FloatComplex (x)) : FloatComplex (acoshf (x));
}

Complex rc_asin (double x)
{
  return fabs (x) > 1.0 ? asin (Complex (x)) : Complex (asin (x));
}

FloatComplex rc_asin (float x)
{
  return fabsf (x) > 1.0f ? asin (FloatComplex (x)) : FloatComplex (asinf (x));
}

Complex rc_atanh (double x)
{
  return fabs (x) > 1.0 ? atanh (Complex (x)) : Complex (atanh (x));
}

FloatComplex rc_atanh (float x)
{
  return fabsf (x) > 1.0f ? atanh (FloatComplex (x)) : FloatComplex (atanhf (x));
}

Complex rc_log (double x)
{
  const double pi = 3.14159265358979323846;
  return x < 0.0 ? Complex (log (-x), pi) : Complex (log (x));
}

FloatComplex rc_log (float x)
{
  const float pi = 3.14159265358979323846f;
  return x < 0.0f ? FloatComplex (logf (-x), pi) : FloatComplex (logf (x));
}

Complex rc_log2 (double x)
{
  const double pil2 = 4.53236014182719380962; // = pi / log(2)
  return x < 0.0 ? Complex (xlog2 (-x), pil2) : Complex (xlog2 (x));
}

FloatComplex rc_log2 (float x)
{
  const float pil2 = 4.53236014182719380962f; // = pi / log(2)
  return x < 0.0f ? FloatComplex (xlog2 (-x), pil2) : FloatComplex (xlog2 (x));
}

Complex rc_log10 (double x)
{
  const double pil10 = 1.36437635384184134748; // = pi / log(10)
  return x < 0.0 ? Complex (log10 (-x), pil10) : Complex (log10 (x));
}

FloatComplex rc_log10 (float x)
{
  const float pil10 = 1.36437635384184134748f; // = pi / log(10)
  return x < 0.0f ? FloatComplex (log10 (-x), pil10) : FloatComplex (log10f (x));
}

Complex rc_sqrt (double x)
{
  return x < 0.0 ? Complex (0.0, sqrt (-x)) : Complex (sqrt (x));
}

FloatComplex rc_sqrt (float x)
{
  return x < 0.0f ? FloatComplex (0.0f, sqrtf (-x)) : FloatComplex (sqrtf (x));
}