/*
Most of the code in this module has been taken from Guile source code.
This file is distributed under GNU Lesser General Public License. 
*/

#include <libguile.h>
#include <math.h>
#include <assert.h>
#include "gc/gc.h"

#define ASSERT(x) assert(x)

#ifndef MINI_GMP_LIMB_TYPE
#define MINI_GMP_LIMB_TYPE long
#endif

typedef unsigned MINI_GMP_LIMB_TYPE mp_limb_t;
typedef long mp_size_t;
typedef unsigned long mp_bitcnt_t;

typedef mp_limb_t *mp_ptr;
typedef const mp_limb_t *mp_srcptr;

struct scm_bignum
{
  scm_t_bits tag;
  /* FIXME: In Guile 3.2, replace this union with just a "size" member.
     Digits are always allocated inline.  */
  union {
    mpz_t mpz;
    struct {
      int zero;
      int size;
      mp_limb_t *limbs;
    } z;
  } u;
  mp_limb_t limbs[];
};

#define NLIMBS_MAX (SSIZE_MAX / sizeof(mp_limb_t))

static SCM _scm_i_normbig (SCM);

static struct scm_bignum *
allocate_bignum (size_t nlimbs)
{
  ASSERT (nlimbs <= (size_t)INT_MAX);
  ASSERT (nlimbs <= NLIMBS_MAX);

  size_t size = sizeof (struct scm_bignum) + nlimbs * sizeof(mp_limb_t);
  struct scm_bignum *z = scm_gc_malloc_pointerless (size, "bignum");

  z->tag = scm_tc16_big;

  z->u.z.zero = 0;
  z->u.z.size = nlimbs;
  z->u.z.limbs = z->limbs;

  // _mp_alloc == 0 means GMP will never try to free this memory.
  ASSERT (z->u.mpz[0]._mp_alloc == 0);
  // Our "size" field should alias the mpz's _mp_size field.
  ASSERT (z->u.mpz[0]._mp_size == nlimbs);
  // Limbs are always allocated inline.
  ASSERT (z->u.mpz[0]._mp_d == z->limbs);

  // z->limbs left uninitialized.
  return z;
}

static struct scm_bignum *
negate_bignum (struct scm_bignum *z)
{
  z->u.z.size = -z->u.z.size;
  return z;
}

static inline unsigned long
long_magnitude (long l)
{
  unsigned long mag = l;
  return l < 0 ? ~mag + 1 : mag;
}

static struct scm_bignum *
make_bignum_0 (void)
{
  return allocate_bignum (0);
}

static struct scm_bignum *
make_bignum_1 (int is_negative, mp_limb_t limb)
{
  struct scm_bignum *z = allocate_bignum (1);
  z->limbs[0] = limb;
  return is_negative ? negate_bignum(z) : z;
}


static inline SCM
scm_from_bignum (struct scm_bignum *x)
{
  return SCM_PACK (x);
}

static struct scm_bignum *
ulong_to_bignum (unsigned long u)
{
  return u == 0 ? make_bignum_0 () : make_bignum_1 (0, u);
};

static struct scm_bignum *
long_to_bignum (long i)
{
  if (i > 0)
    return ulong_to_bignum (i);

  return i == 0 ? make_bignum_0 () : make_bignum_1 (1, long_magnitude (i));
};

static SCM
long_to_scm (long i)
{
  if (SCM_FIXABLE (i))
    return SCM_I_MAKINUM (i);
  return scm_from_bignum (long_to_bignum (i));
}

static mp_limb_t*
bignum_limbs (struct scm_bignum *z)
{
  // FIXME: In the future we can just return z->limbs.
  return z->u.z.limbs;
}

static struct scm_bignum *
make_bignum_from_mpz (mpz_srcptr mpz)
{
  size_t nlimbs = mpz_size (mpz);
  struct scm_bignum *ret = allocate_bignum (nlimbs);
  mpn_copyi (bignum_limbs (ret), mpz_limbs_read (mpz), nlimbs);
  return mpz_sgn (mpz) < 0 ? negate_bignum (ret) : ret;
}

static SCM
take_mpz (mpz_ptr mpz)
{
  SCM ret;
  if (mpz_fits_slong_p (mpz))
    ret = long_to_scm (mpz_get_si (mpz));
  else
    ret = scm_from_bignum (make_bignum_from_mpz (mpz));
  mpz_clear (mpz);
  return ret;
}

SCM
scm_integer_from_double (double val)
{
  if (!isfinite (val))
    scm_out_of_range ("inexact->exact", scm_from_double (val));

  if (((double) INT64_MIN) <= val && val <= ((double) INT64_MAX))
    return scm_from_int64 (val);

  mpz_t result;
  mpz_init_set_d (result, val);
  return take_mpz (result);
}

static int convert_to_integer( SCM x )
  {
  if( SCM_LIKELY( SCM_I_INUMP( x ) ) )
    {
    return SCM_I_INUM( x );
    }
  else if( SCM_BIGP( x ) )
    {
    SCM xNormed = _scm_i_normbig( x );
    if( SCM_I_INUMP( xNormed ) )
      {
      return SCM_I_INUM( xNormed );
      }
    else
      {
      scm_wrong_type_arg_msg( "convert_to_integer", SCM_ARG1, x,
                              "inum");
      }
    }
  else
    {
    scm_wrong_type_arg_msg( "convert_to_integer", SCM_ARG1, x,
                            "inum");
    }
  }

static SCM
_scm_i_normbig (SCM b)
{
  /* convert a big back to a fixnum if it'll fit */
  /* presume b is a bignum */
  if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b)))
    {
      scm_t_inum val = mpz_get_si (SCM_I_BIG_MPZ (b));
      if (SCM_FIXABLE (val))
        b = SCM_I_MAKINUM (val);
    }
  return b;
}

#define MAKE_WRAPPER(name, func) \
  SCM name( SCM x ) \
      { \
      return scm_from_double( func ( scm_to_double( x ) ) ); \
      }

#define MAKE_WRAPPER2(name, func) \
  SCM name( SCM x1, SCM x2 )             \
      { \
      return scm_from_double( func ( scm_to_double( x1 ), \
                                     scm_to_double( x2 ) ) );   \
      }

static SCM gl_b_fast_fma = SCM_UNSPECIFIED;
static SCM gl_i_ilogb0 = SCM_UNSPECIFIED;
static SCM gl_i_ilogbnan = SCM_UNSPECIFIED;
static SCM gl_i_nan = SCM_UNSPECIFIED;
static SCM gl_i_infinite = SCM_UNSPECIFIED;
static SCM gl_i_zero = SCM_UNSPECIFIED;
static SCM gl_i_subnormal = SCM_UNSPECIFIED;
static SCM gl_i_normal = SCM_UNSPECIFIED;

SCM is_integer_wrapper( SCM x )
  {
  return ( SCM_I_INUMP( x ) || SCM_BIGP( x ) ) ? SCM_BOOL_T : SCM_BOOL_F;
  }

SCM is_real_wrapper( SCM x )
  {
  return SCM_REALP( x ) ? SCM_BOOL_T : SCM_BOOL_F;
  }

SCM real_to_integer_wrapper( SCM x )
  {
  if( SCM_LIKELY( SCM_REALP( x ) ) )
    {
    double r = SCM_REAL_VALUE( x );
    if( r == floor( r ) )
      {
      return scm_integer_from_double( r );
      }
    else
      {
      scm_misc_error(
        "real_to_integer_wrapper", "argument not integer valued", SCM_EOL
        );
      }
    }
  else
    {
    scm_wrong_type_arg_msg( "real_to_integer_wrapper", SCM_ARG1, x, "real");
    }
  }

MAKE_WRAPPER( r_round_wrapper, round )
MAKE_WRAPPER( r_trunc_wrapper, trunc )
MAKE_WRAPPER( r_floor_wrapper, floor )
MAKE_WRAPPER( r_ceil_wrapper, ceil )

MAKE_WRAPPER( r_sqrt_wrapper, sqrt )
MAKE_WRAPPER( r_sin_wrapper, sin )
MAKE_WRAPPER( r_cos_wrapper, cos )
MAKE_WRAPPER( r_tan_wrapper, tan )
MAKE_WRAPPER( r_asin_wrapper, asin )
MAKE_WRAPPER( r_acos_wrapper, acos )
MAKE_WRAPPER( r_atan_wrapper, atan )
MAKE_WRAPPER( r_exp_wrapper, exp )
MAKE_WRAPPER( r_log_wrapper, log )
MAKE_WRAPPER( r_log10_wrapper, log10 )
MAKE_WRAPPER( r_sinh_wrapper, sinh )
MAKE_WRAPPER( r_cosh_wrapper, cosh )
MAKE_WRAPPER( r_tanh_wrapper, tanh )
MAKE_WRAPPER( r_asinh_wrapper, asinh )
MAKE_WRAPPER( r_acosh_wrapper, acosh )
MAKE_WRAPPER( r_atanh_wrapper, atanh )

MAKE_WRAPPER2( r_expt_wrapper, pow )
MAKE_WRAPPER2( r_atan2_wrapper, atan2 )

#ifdef STANDARD_MATH
MAKE_WRAPPER2( fmod_wrapper, fmod )
MAKE_WRAPPER2( r_remainder_wrapper, remainder );
MAKE_WRAPPER2( fmin_wrapper, fmin )
MAKE_WRAPPER2( fmax_wrapper, fmax )
MAKE_WRAPPER2( fdim_wrapper, fdim )
MAKE_WRAPPER( r_exp2_wrapper, exp2 )
MAKE_WRAPPER( r_expm1_wrapper, expm1 )
MAKE_WRAPPER( r_log2_wrapper, log2 )
MAKE_WRAPPER( r_log1p_wrapper, log1p )
MAKE_WRAPPER( logb_wrapper, logb )
MAKE_WRAPPER( r_cbrt_wrapper, cbrt )
MAKE_WRAPPER2( r_hypot_wrapper, hypot )
MAKE_WRAPPER( r_erf_wrapper, erf )
MAKE_WRAPPER( r_erfc_wrapper, erfc )
MAKE_WRAPPER( r_lgamma_wrapper, lgamma )
MAKE_WRAPPER( r_tgamma_wrapper, tgamma )
MAKE_WRAPPER( r_nearbyint_wrapper, nearbyint )
MAKE_WRAPPER( rint_wrapper, rint )
MAKE_WRAPPER2( r_nextafter_wrapper, nextafter )
MAKE_WRAPPER2( r_copysign_wrapper, copysign )
#endif

#ifdef STANDARD_MATH
SCM frexp_wrapper( SCM x )
  {
  double r = scm_to_double( x );
  int iExp = 0;
  double rFraction = frexp( r, &iExp );
  SCM xExp = scm_from_signed_integer( iExp );
  return scm_cons( scm_from_double( rFraction ), xExp );
  }
#endif

#ifdef STANDARD_MATH
SCM ldexp_wrapper( SCM x1, SCM x2 )
  {
  double rBase = scm_to_double( x1 );
  int iExp = convert_to_integer( x2 );
  return scm_from_double( ldexp( rBase, iExp ) );
  }
#endif

#ifdef STANDARD_MATH
SCM modf_wrapper( SCM x )
  {
  double r = scm_to_double( x );
  double rInt = 0;
  double rFraction = modf( r, &rInt );
  return scm_cons( scm_from_double( rFraction ), scm_from_double( rInt ) );
  }
#endif

#ifdef STANDARD_MATH
SCM fma_wrapper( SCM x1, SCM x2, SCM x3 )
  {
  return scm_from_double( fma( scm_to_double( x1 ),
                               scm_to_double( x2 ),
                               scm_to_double( x3 ) ) );
  }
#endif

#ifdef STANDARD_MATH
SCM fast_fma_wrapper( void )
  {
  return gl_b_fast_fma;
  }
#endif

#ifdef STANDARD_MATH
SCM ilogb_wrapper( SCM x )
  {
  int iResult = ilogb( scm_to_double( x ) );
  return scm_from_signed_integer( iResult );
  }
#endif

#ifdef STANDARD_MATH
SCM ilogb0_wrapper( void )
  {
  return gl_i_ilogb0;
  }
#endif

#ifdef STANDARD_MATH
SCM ilogbnan_wrapper( void )
  {
  return gl_i_ilogbnan;
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_wrapper( SCM x )
  {
  int iResult = fpclassify( scm_to_double( x ) );
  return scm_from_signed_integer( iResult );
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_nan_wrapper( void )
  {
  return gl_i_nan;
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_infinite_wrapper( void )
  {
  return gl_i_infinite;
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_zero_wrapper( void )
  {
  return gl_i_zero;
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_subnormal_wrapper( void )
  {
  return gl_i_subnormal;
  }
#endif

#ifdef STANDARD_MATH
SCM fpclassify_normal_wrapper( void )
  {
  return gl_i_normal;
  }
#endif

#ifdef STANDARD_MATH
SCM r_isnormal_wrapper( SCM x )
  {
  return isnormal( scm_to_double( x ) ) ? SCM_BOOL_T : SCM_BOOL_F;
  }
#endif

#ifdef STANDARD_MATH
SCM r_signbit_wrapper( SCM x )
  {
  return ( signbit( scm_to_double( x ) ) != 0 ) ? SCM_INUM1 : SCM_INUM0;
  }
#endif

#ifdef POSIX_MATH
MAKE_WRAPPER( r_j0_wrapper, j0 )
MAKE_WRAPPER( r_j1_wrapper, j1 )
MAKE_WRAPPER( r_y0_wrapper, y0 )
MAKE_WRAPPER( r_y1_wrapper, y1 )
#endif

#ifdef POSIX_MATH
SCM r_jn_wrapper( SCM xN, SCM xArg )
  {
  int i = convert_to_integer( xN );
  double r = scm_to_double( xArg );
  return scm_from_double( jn( i, r ) );
  }
#endif

#ifdef POSIX_MATH
SCM r_yn_wrapper( SCM xN, SCM xArg )
  {
  int i = convert_to_integer( xN );
  double r = scm_to_double( xArg );
  return scm_from_double( yn( i, r ) );
  }
#endif

void init_theme_d_support( void )
  {
#ifdef STANDARD_MATH
#ifdef FP_FAST_FMA
  gl_b_fast_fma = SCM_BOOL_T;
#else
  gl_b_fast_fma = SCM_BOOL_F;
#endif
  gl_i_ilogb0 = scm_from_signed_integer( FP_ILOGB0 );
  gl_i_ilogbnan = scm_from_signed_integer( FP_ILOGBNAN );
  gl_i_nan = scm_from_signed_integer( FP_NAN ); 
  gl_i_infinite = scm_from_signed_integer( FP_INFINITE ); 
  gl_i_zero = scm_from_signed_integer( FP_ZERO ); 
  gl_i_subnormal = scm_from_signed_integer( FP_SUBNORMAL ); 
  gl_i_normal = scm_from_signed_integer( FP_NORMAL ); 
#endif
  scm_c_define_gsubr( "is-integer?", 1, 0, 0, is_integer_wrapper );
  scm_c_define_gsubr( "is-real?", 1, 0, 0, is_real_wrapper );
  scm_c_define_gsubr( "real->integer", 1, 0, 0, real_to_integer_wrapper );
  scm_c_define_gsubr( "r-round", 1, 0, 0, r_round_wrapper );
  scm_c_define_gsubr( "r-truncate", 1, 0, 0, r_trunc_wrapper );
  scm_c_define_gsubr( "r-floor", 1, 0, 0, r_floor_wrapper );
  scm_c_define_gsubr( "r-ceiling", 1, 0, 0, r_ceil_wrapper );
  scm_c_define_gsubr( "r-sqrt", 1, 0, 0, r_sqrt_wrapper );
  scm_c_define_gsubr( "r-sin", 1, 0, 0, r_sin_wrapper );
  scm_c_define_gsubr( "r-cos", 1, 0, 0, r_cos_wrapper );
  scm_c_define_gsubr( "r-tan", 1, 0, 0, r_tan_wrapper );
  scm_c_define_gsubr( "r-asin", 1, 0, 0, r_asin_wrapper );
  scm_c_define_gsubr( "r-acos", 1, 0, 0, r_acos_wrapper );
  scm_c_define_gsubr( "r-atan", 1, 0, 0, r_atan_wrapper );
  scm_c_define_gsubr( "r-exp", 1, 0, 0, r_exp_wrapper );
  scm_c_define_gsubr( "r-log", 1, 0, 0, r_log_wrapper );
  scm_c_define_gsubr( "r-log10", 1, 0, 0, r_log10_wrapper );
  scm_c_define_gsubr( "r-sinh", 1, 0, 0, r_sinh_wrapper );
  scm_c_define_gsubr( "r-cosh", 1, 0, 0, r_cosh_wrapper );
  scm_c_define_gsubr( "r-tanh", 1, 0, 0, r_tanh_wrapper );
  scm_c_define_gsubr( "r-asinh", 1, 0, 0, r_asinh_wrapper );
  scm_c_define_gsubr( "r-acosh", 1, 0, 0, r_acosh_wrapper );
  scm_c_define_gsubr( "r-atanh", 1, 0, 0, r_atanh_wrapper );
  scm_c_define_gsubr( "r-expt", 2, 0, 0, r_expt_wrapper );
  scm_c_define_gsubr( "r-atan2", 2, 0, 0, r_atan2_wrapper );
#ifdef STANDARD_MATH
  scm_c_define_gsubr( "fmod", 2, 0, 0, fmod_wrapper );
  scm_c_define_gsubr( "r-remainder", 2, 0, 0, r_remainder_wrapper );
  scm_c_define_gsubr( "fma", 3, 0, 0, fma_wrapper );
  scm_c_define_gsubr( "fast-fma?", 0, 0, 0, fast_fma_wrapper );
  scm_c_define_gsubr( "fmin", 2, 0, 0, fmin_wrapper );
  scm_c_define_gsubr( "fmax", 2, 0, 0, fmax_wrapper );
  scm_c_define_gsubr( "fdim", 2, 0, 0, fdim_wrapper );
  scm_c_define_gsubr( "r-exp2", 1, 0, 0, r_exp2_wrapper ); 
  scm_c_define_gsubr( "r-expm1", 1, 0, 0, r_expm1_wrapper ); 
  scm_c_define_gsubr( "r-log2", 1, 0, 0, r_log2_wrapper ); 
  scm_c_define_gsubr( "r-log1p", 1, 0, 0, r_log1p_wrapper ); 
  scm_c_define_gsubr( "logb", 1, 0, 0, logb_wrapper ); 
  scm_c_define_gsubr( "ilogb", 1, 0, 0, ilogb_wrapper );
  scm_c_define_gsubr( "ilogb0", 0, 0, 0, ilogb0_wrapper );
  scm_c_define_gsubr( "ilogbnan", 0, 0, 0, ilogbnan_wrapper );
  scm_c_define_gsubr( "r-cbrt", 1, 0, 0, r_cbrt_wrapper );
  scm_c_define_gsubr( "r-hypot", 2, 0, 0, r_hypot_wrapper );
  scm_c_define_gsubr( "r-erf", 1, 0, 0, r_erf_wrapper );
  scm_c_define_gsubr( "r-erfc", 1, 0, 0, r_erfc_wrapper );
  scm_c_define_gsubr( "r-lgamma", 1, 0, 0, r_lgamma_wrapper );
  scm_c_define_gsubr( "r-tgamma", 1, 0, 0, r_tgamma_wrapper );
  scm_c_define_gsubr( "r-nearbyint", 1, 0, 0, r_nearbyint_wrapper );
  scm_c_define_gsubr( "rint", 1, 0, 0, rint_wrapper );
  scm_c_define_gsubr( "frexp", 1, 0, 0, frexp_wrapper );
  scm_c_define_gsubr( "ldexp", 2, 0, 0, ldexp_wrapper );
  scm_c_define_gsubr( "modf", 1, 0, 0, modf_wrapper );
  scm_c_define_gsubr( "r-nextafter", 2, 0, 0, r_nextafter_wrapper );
  scm_c_define_gsubr( "r-copysign", 2, 0, 0, r_copysign_wrapper );
  scm_c_define_gsubr( "fpclassify", 1, 0, 0, fpclassify_wrapper );
  scm_c_define_gsubr( "fpclassify-nan", 0, 0, 0, fpclassify_nan_wrapper );
  scm_c_define_gsubr( "fpclassify-infinite", 0, 0, 0,
                      fpclassify_infinite_wrapper );
  scm_c_define_gsubr( "fpclassify-zero", 0, 0, 0, fpclassify_zero_wrapper );
  scm_c_define_gsubr( "fpclassify-subnormal", 0, 0, 0,
                      fpclassify_subnormal_wrapper );
  scm_c_define_gsubr( "fpclassify-normal", 0, 0, 0,
                      fpclassify_normal_wrapper );
  scm_c_define_gsubr( "r-isnormal?", 1, 0, 0, r_isnormal_wrapper );
  scm_c_define_gsubr( "r-signbit", 1, 0, 0, r_signbit_wrapper );
#endif
#ifdef POSIX_MATH
  scm_c_define_gsubr( "r-j0", 1, 0, 0, r_j0_wrapper );
  scm_c_define_gsubr( "r-j1", 1, 0, 0, r_j1_wrapper );
  scm_c_define_gsubr( "r-jn", 2, 0, 0, r_jn_wrapper );
  scm_c_define_gsubr( "r-y0", 1, 0, 0, r_y0_wrapper );
  scm_c_define_gsubr( "r-y1", 1, 0, 0, r_y1_wrapper );
  scm_c_define_gsubr( "r-yn", 2, 0, 0, r_yn_wrapper );
#endif
  }