/* -*- mode: C -*-  */
/* 
   IGraph library.
   Copyright (C) 2007-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA
   
   This program 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 2 of the License, or
   (at your option) any later version.
   
   This program 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 this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
   02110-1301 USA

*/

#include <math.h>
#include <float.h>
#include <stdarg.h>
#include "config.h"
#include "igraph_math.h"
#include "igraph_types.h"

#ifdef _MSC_VER
#  define isinf(x) (!_finite(x) && !_isnan(x))
#endif

int igraph_finite(double x)
{
#ifdef isfinite 
    return isfinite(x);
#elif HAVE_ISFINITE == 1
    return isfinite(x);
#elif HAVE_FINITE == 1
    return finite(x);
#else
/* neither finite nor isfinite work. Do we really need the AIX exception? */
# ifdef _AIX
#  include <fp.h>
     return FINITE(x);
# else
    return (!isnan(x) & (x != IGRAPH_POSINFINITY) & (x != IGRAPH_NEGINFINITY));
# endif
#endif
}

double igraph_log2(const double a) {
  return log(a)/log(2.0);
}

int igraph_chebyshev_init(const double *dos, int nos, double eta)
{
    int i, ii;
    double err;

    if (nos < 1)
	return 0;

    err = 0.0;
    i = 0;			/* just to avoid compiler warnings */
    for (ii=1; ii<=nos; ii++) {
	i = nos - ii;
	err += fabs(dos[i]);
	if (err > eta) {
	    return i;
	}
    }
    return i;
}

double igraph_chebyshev_eval(double x, const double *a, const int n)
{
    double b0, b1, b2, twox;
    int i;

    if (n < 1 || n > 1000) IGRAPH_NAN;

    if (x < -1.1 || x > 1.1) IGRAPH_NAN;

    twox = x * 2;
    b2 = b1 = 0;
    b0 = 0;
    for (i = 1; i <= n; i++) {
	b2 = b1;
	b1 = b0;
	b0 = twox * b1 - b2 + a[n - i];
    }
    return (b0 - b2) * 0.5;
}

double igraph_log1p(double x)
{
    /* series for log1p on the interval -.375 to .375
     *				     with weighted error   6.35e-32
     *				      log weighted error  31.20
     *			    significant figures required  30.93
     *				 decimal places required  32.01
     */
    static const double alnrcs[43] = {
	+.10378693562743769800686267719098e+1,
	-.13364301504908918098766041553133e+0,
	+.19408249135520563357926199374750e-1,
	-.30107551127535777690376537776592e-2,
	+.48694614797154850090456366509137e-3,
	-.81054881893175356066809943008622e-4,
	+.13778847799559524782938251496059e-4,
	-.23802210894358970251369992914935e-5,
	+.41640416213865183476391859901989e-6,
	-.73595828378075994984266837031998e-7,
	+.13117611876241674949152294345011e-7,
	-.23546709317742425136696092330175e-8,
	+.42522773276034997775638052962567e-9,
	-.77190894134840796826108107493300e-10,
	+.14075746481359069909215356472191e-10,
	-.25769072058024680627537078627584e-11,
	+.47342406666294421849154395005938e-12,
	-.87249012674742641745301263292675e-13,
	+.16124614902740551465739833119115e-13,
	-.29875652015665773006710792416815e-14,
	+.55480701209082887983041321697279e-15,
	-.10324619158271569595141333961932e-15,
	+.19250239203049851177878503244868e-16,
	-.35955073465265150011189707844266e-17,
	+.67264542537876857892194574226773e-18,
	-.12602624168735219252082425637546e-18,
	+.23644884408606210044916158955519e-19,
	-.44419377050807936898878389179733e-20,
	+.83546594464034259016241293994666e-21,
	-.15731559416479562574899253521066e-21,
	+.29653128740247422686154369706666e-22,
	-.55949583481815947292156013226666e-23,
	+.10566354268835681048187284138666e-23,
	-.19972483680670204548314999466666e-24,
	+.37782977818839361421049855999999e-25,
	-.71531586889081740345038165333333e-26,
	+.13552488463674213646502024533333e-26,
	-.25694673048487567430079829333333e-27,
	+.48747756066216949076459519999999e-28,
	-.92542112530849715321132373333333e-29,
	+.17578597841760239233269760000000e-29,
	-.33410026677731010351377066666666e-30,
	+.63533936180236187354180266666666e-31,
    };

    static IGRAPH_THREAD_LOCAL int nlnrel = 0;
    static IGRAPH_THREAD_LOCAL double xmin = 0.0;

    if (xmin == 0.0) xmin = -1 + sqrt(DBL_EPSILON);/*was sqrt(d1mach(4)); */
    if (nlnrel == 0) /* initialize chebychev coefficients */
	nlnrel = igraph_chebyshev_init(alnrcs, 43, DBL_EPSILON/20);/*was .1*d1mach(3)*/

    if (x == 0.) return 0.;/* speed */
    if (x == -1) return(IGRAPH_NEGINFINITY);
    if (x  < -1) return(IGRAPH_NAN);

    if (fabs(x) <= .375) {
        /* Improve on speed (only);
	   again give result accurate to IEEE double precision: */
	if(fabs(x) < .5 * DBL_EPSILON)
	    return x;

	if( (0 < x && x < 1e-8) || (-1e-9 < x && x < 0))
	    return x * (1 - .5 * x);
	/* else */
	return x * (1 - x * igraph_chebyshev_eval(x / .375, alnrcs, nlnrel));
    }
    /* else */
/*     if (x < xmin) { */
/* 	/\* answer less than half precision because x too near -1 *\/ */
/*         ML_ERROR(ME_PRECISION, "log1p"); */
/*     } */
    return log(1 + x);
}

long double igraph_fabsl(long double a) { 
  if (a<0) {
    return -a;
  } else {
    return a;
  }
}

double igraph_fmin(double a, double b) { 
  if (b<a) {
    return b;
  } else {
    return a;
  }
}

double igraph_i_round(double X) {
  
  /* NaN */
  if (X != X) { return X; }

  if (X < 0.0) { 
    return floor(X);
  }

  return ceil(X);
}

#ifdef _MSC_VER
/**
 * Internal function, replacement for snprintf
 * Used only in case of the Microsoft Visual C compiler which does not
 * provide a proper sprintf implementation.
 * 
 * This implementation differs from the standard in the value returned
 * when the number of characters needed by the output, excluding the
 * terminating '\0' is larger than count
 */
int igraph_i_snprintf(char *buffer, size_t count, const char *format, ...) {
    int n;
    va_list args;
    if (count > 0) {
	va_start(args, format);
        n = _vsnprintf(buffer, count, format, args);
        buffer[count-1] = 0;
        va_end(args);
    } else n=0;
    return n;
}

#endif

int igraph_is_nan(double x) {
  return isnan(x);
}

int igraph_is_inf(double x) {
  return isinf(x) != 0;
}

int igraph_is_posinf(double x) {
  return isinf(x) == 1;
}

int igraph_is_neginf(double x) {
  return isinf(x) == -1;
}