/*

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 
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
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 -> 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 
xroundb (float x)
{
  float t = xround (x);

  if (fabsf (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
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);
}

// 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));
}

bool
xnegative_sign (double x)
{
  return __lo_ieee_signbit (x);
}

bool
xnegative_sign (float x)
{
  return __lo_ieee_float_signbit (x);
}

// Convert X to the nearest integer value.  Should not pass NaN to
// this function.

// Sometimes you need a large integer, but not always.

octave_idx_type
NINTbig (double x)
{
  if (x > std::numeric_limits<octave_idx_type>::max ())
    return std::numeric_limits<octave_idx_type>::max ();
  else if (x < std::numeric_limits<octave_idx_type>::min ())
    return std::numeric_limits<octave_idx_type>::min ();
  else
    return static_cast<octave_idx_type> ((x > 0) ? (x + 0.5) : (x - 0.5));
}

octave_idx_type
NINTbig (float x)
{
  if (x > std::numeric_limits<octave_idx_type>::max ())
    return std::numeric_limits<octave_idx_type>::max ();
  else if (x < std::numeric_limits<octave_idx_type>::min ())
    return std::numeric_limits<octave_idx_type>::min ();
  else
    return static_cast<octave_idx_type> ((x > 0) ? (x + 0.5) : (x - 0.5));
}

int
NINT (double x)
{
  if (x > std::numeric_limits<int>::max ())
    return std::numeric_limits<int>::max ();
  else if (x < std::numeric_limits<int>::min ())
    return std::numeric_limits<int>::min ();
  else
    return static_cast<int> ((x > 0) ? (x + 0.5) : (x - 0.5));
}

int
NINT (float x)
{
  if (x > std::numeric_limits<int>::max ())
    return std::numeric_limits<int>::max ();
  else if (x < std::numeric_limits<int>::min ())
    return std::numeric_limits<int>::min ();
  else
    return static_cast<int> ((x > 0) ? (x + 0.5) : (x - 0.5));
}