/****************************************************************************
**
**  This file is part of GAP, a system for computational discrete algebra.
**
**  Copyright of GAP belongs to its developers, whose names are too numerous
**  to list here. Please refer to the COPYRIGHT file for details.
**
**  SPDX-License-Identifier: GPL-2.0-or-later
**
**  This file contains the functions for the macfloat package.
**
**  Machine floating point values, aka macfloats, are stored as bags
**  containing a 64 bit value.
*/

// glibc only declares exp10 in its headers if we define _GNU_SOURCE
#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif

#include <math.h>

#include "macfloat.h"

#include "ariths.h"
#include "bool.h"
#include "error.h"
#include "integer.h"
#include "io.h"
#include "modules.h"
#include "plist.h"
#include "saveload.h"
#include "stringobj.h"
#include "sysstr.h"

#include "config.h"

#include <stdio.h>      // for snprintf
#include <stdlib.h>     // for strtod


#define RequireMacFloat(funcname, op) \
    RequireArgumentCondition(funcname, op, IS_MACFLOAT(op), \
        "must be a macfloat")


Double VAL_MACFLOAT(Obj obj)
{
    Double val;
    memcpy(&val, CONST_ADDR_OBJ(obj), sizeof(Double));
    return val;
}

void SET_VAL_MACFLOAT(Obj obj, Double val)
{
    memcpy(ADDR_OBJ(obj), &val, sizeof(Double));
}

BOOL IS_MACFLOAT(Obj obj)
{
    return TNUM_OBJ(obj) == T_MACFLOAT;
}


/****************************************************************************
**
*F  TypeMacfloat( <macfloat> )  . . . . . . . . . .  type of a macfloat value
**
**  'TypeMacfloat' returns the type of macfloatean values.
**
**  'TypeMacfloat' is the function in 'TypeObjFuncs' for macfloatean values.
*/
static Obj TYPE_MACFLOAT;

static Obj TypeMacfloat(Obj val)
{
    return TYPE_MACFLOAT;
}


// helper function for printing a "decimal" representation of a macfloat
// into a buffer.
static void
PrintMacfloatToBuf(char * buf, size_t bufsize, Double val, int precision)
{
    // handle printing of NaN and infinities ourselves, to ensure
    // they are printed uniformly across all platforms
    if (isnan(val)) {
        strcpy(buf, "nan");
    }
    else if (isinf(val)) {
        if (val > 0)
            strcpy(buf, "inf");
        else
            strcpy(buf, "-inf");
    }
    else {
        snprintf(buf, bufsize, "%.*" PRINTFFORMAT, precision, val);
        // check if a period is in the output; this is not always the case,
        // e.g. if the value is an integer
        if (strchr(buf, '.'))
            return;    // everything is fine
        // we need to insert a '.'; either at the end, or before an exponent
        // (e.g. "7e10" -> "7.e10"). For this we need 1 extra byte of storage,
        // plus of course 1 byte for the string terminator; check if the
        // buffer is big enough
        if (strlen(buf) + 2 <= bufsize) {
            char * loc = strchr(buf, 'e');
            if (loc) {
                memmove(loc + 1, loc, strlen(loc) + 1);
                loc[0] = '.';
            }
            else {
                strxcat(buf, ".", bufsize);
            }
        }
    }
}


/****************************************************************************
**
*F  PrintMacfloat( <macfloat> ) . . . . . . . . . . .  print a macfloat value
**
**  'PrintMacfloat' prints the macfloating value <macfloat>.
*/
static void PrintMacfloat(Obj x)
{
    Char buf[1024];
    // TODO: should we use PRINTFDIGITS instead of 16?
    PrintMacfloatToBuf(buf, sizeof(buf), VAL_MACFLOAT(x), 16);
    Pr("%s", (Int)buf, 0);
}


/****************************************************************************
**
*F  EqMacfloat( <macfloatL>, <macfloatR> )  . . . . . . . . .  test if <macfloatL> =  <macfloatR>
**
**  'EqMacfloat' returns 'True' if the two macfloatean values <macfloatL> and <macfloatR> are
**  equal, and 'False' otherwise.
*/
static Int EqMacfloat(Obj macfloatL, Obj macfloatR)
{
  return VAL_MACFLOAT(macfloatL) == VAL_MACFLOAT(macfloatR);
}

static Obj FuncEQ_MACFLOAT(Obj self, Obj macfloatL, Obj macfloatR)
{
  return EqMacfloat(macfloatL,macfloatR) ? True : False;
}


/****************************************************************************
**
*F  LtMacfloat( <macfloatL>, <macfloatR> )  . . . . . . . . .  test if <macfloatL> <  <macfloatR>
**
*/
static Int LtMacfloat(Obj macfloatL, Obj macfloatR)
{
  return VAL_MACFLOAT(macfloatL) < VAL_MACFLOAT(macfloatR);
}


/****************************************************************************
**
*F  SaveMacfloat( <macfloat> ) . . . . . . . . . . . . . . save a Macfloatean
**
*/
#ifdef GAP_ENABLE_SAVELOAD
static void SaveMacfloat(Obj obj)
{
    const UInt1 *data = (const UInt1 *)CONST_ADDR_OBJ(obj);
    for (UInt i = 0; i < sizeof(Double); i++)
        SaveUInt1(data[i]);
}
#endif


/****************************************************************************
**
*F  LoadMacfloat( <macfloat> ) . . . . . . . . . . . . . . load a Macfloatean
**
*/
#ifdef GAP_ENABLE_SAVELOAD
static void LoadMacfloat(Obj obj)
{
    UInt1 *data = (UInt1 *)ADDR_OBJ(obj);
    for (UInt i = 0; i < sizeof(Double); i++)
        data[i] = LoadUInt1();
}
#endif


Obj NEW_MACFLOAT( Double val )
{
    Obj f = NewBag(T_MACFLOAT, sizeof(Double));
    SET_VAL_MACFLOAT(f, val);
    return f;
}

/****************************************************************************
**
*F  ZeroMacfloat(<macfloat> ) . . . . . . . . . . . . . . . . . . . return the zero
**
*/


static Obj ZeroMacfloat(Obj f)
{
  return NEW_MACFLOAT((Double)0.0);
}

/****************************************************************************
**
*F  AInvMacfloat( <macfloat> ) . . . . . . . . . . . . . . . . .  unary minus
**
*/


static Obj AInvMacfloat(Obj f)
{
  return NEW_MACFLOAT(-VAL_MACFLOAT(f));
}

/****************************************************************************
**
*F  OneMacfloat(<macfloat> ) . . . . . . . . . . . . . . . . . . . return the one
**
*/


static Obj OneMacfloat(Obj f)
{
  return NEW_MACFLOAT((Double)1.0);
}

/****************************************************************************
**
*F  InvMacfloat(<macfloat> ) . . . . . . . . . . . . . . . . . . . reciprocal
**
*/


static Obj InvMacfloat(Obj f)
{
  return NEW_MACFLOAT((Double)1.0/VAL_MACFLOAT(f));
}

/****************************************************************************
**
*F  ProdMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . . . product
**
*/


static Obj ProdMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(VAL_MACFLOAT(fl)*VAL_MACFLOAT(fr));
}

/****************************************************************************
**
*F  PowMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . . exponentiation
**
*/


static Obj PowMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(MATH(pow)(VAL_MACFLOAT(fl),VAL_MACFLOAT(fr)));
}

/****************************************************************************
**
*F  SumMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . .  sum
**
*/


static Obj SumMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(VAL_MACFLOAT(fl)+VAL_MACFLOAT(fr));
}

/****************************************************************************
**
*F  DiffMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . . difference
**
*/


static Obj DiffMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(VAL_MACFLOAT(fl)-VAL_MACFLOAT(fr));
}

/****************************************************************************
**
*F  QuoMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . . quotient
**
*/


static Obj QuoMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(VAL_MACFLOAT(fl)/VAL_MACFLOAT(fr));
}

/****************************************************************************
**
*F  LQuoMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . .left quotient
**
*/


static Obj LQuoMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(VAL_MACFLOAT(fr)/VAL_MACFLOAT(fl));
}

/****************************************************************************
**
*F  ModMacfloat(<macfloatl>, <macfloatr> ) . . . . . . . . . . . . . . .mod
**
*/


static Obj ModMacfloat(Obj fl, Obj fr)
{
  return NEW_MACFLOAT(MATH(fmod)(VAL_MACFLOAT(fl),VAL_MACFLOAT(fr)));
}


/****************************************************************************
**
*F  FuncMACFLOAT_INT(<int>) . . . . . . . . . . . . . . . conversion
**
*/

static Obj FuncMACFLOAT_INT(Obj self, Obj i)
{
  if (!IS_INTOBJ(i))
    return Fail;
  else
    return NEW_MACFLOAT((Double)INT_INTOBJ(i));
}

/****************************************************************************
**
*F  FuncMACFLOAT_STRING(<string>) . . . . . . . . . . . . . . . conversion
**
*/

static Obj FuncMACFLOAT_STRING(Obj self, Obj s)
{
    RequireStringRep(SELF_NAME, s);

  char * endptr;
  UChar *sp = CHARS_STRING(s);
  Obj res= NEW_MACFLOAT((Double) STRTOD((char *)sp,&endptr));
  if ((UChar *)endptr != sp + GET_LEN_STRING(s))
    return Fail;
  return res;
}

/****************************************************************************
**
*F SumIntMacfloat( <int>, <macfloat> )
**
*/

static Obj SumIntMacfloat(Obj i, Obj f)
{
  return NEW_MACFLOAT( (Double)(INT_INTOBJ(i)) + VAL_MACFLOAT(f));
}


/****************************************************************************
**
*F FuncSIN_MACFLOAT( <self>, <macfloat> ) . .The sin function from the math library
**
*/

#define MAKEMATHPRIMITIVE(NAME, name)                                        \
    static Obj Func##NAME##_MACFLOAT(Obj self, Obj f)                        \
    {                                                                        \
        return NEW_MACFLOAT(MATH(name)(VAL_MACFLOAT(f)));                    \
    }

#define MAKEMATHPRIMITIVE2(NAME, name)                                       \
    static Obj Func##NAME##_MACFLOAT(Obj self, Obj f, Obj g)                 \
    {                                                                        \
        return NEW_MACFLOAT(MATH(name)(VAL_MACFLOAT(f), VAL_MACFLOAT(g)));   \
    }

MAKEMATHPRIMITIVE(COS,cos)
MAKEMATHPRIMITIVE(SIN,sin)
MAKEMATHPRIMITIVE(TAN,tan)
MAKEMATHPRIMITIVE(ACOS,acos)
MAKEMATHPRIMITIVE(ASIN,asin)
MAKEMATHPRIMITIVE(ATAN,atan)

MAKEMATHPRIMITIVE(COSH,cosh)
MAKEMATHPRIMITIVE(SINH,sinh)
MAKEMATHPRIMITIVE(TANH,tanh)
MAKEMATHPRIMITIVE(ACOSH,acosh)
MAKEMATHPRIMITIVE(ASINH,asinh)
MAKEMATHPRIMITIVE(ATANH,atanh)

MAKEMATHPRIMITIVE(LOG,log)
MAKEMATHPRIMITIVE(LOG2,log2)
MAKEMATHPRIMITIVE(LOG10,log10)
MAKEMATHPRIMITIVE(LOG1P,log1p)

MAKEMATHPRIMITIVE(EXP,exp)
MAKEMATHPRIMITIVE(EXP2,exp2)
MAKEMATHPRIMITIVE(EXPM1,expm1)
#ifdef HAVE_EXP10
MAKEMATHPRIMITIVE(EXP10,exp10)
#endif

MAKEMATHPRIMITIVE(SQRT,sqrt)
MAKEMATHPRIMITIVE(CBRT,cbrt)
MAKEMATHPRIMITIVE(RINT,rint)
MAKEMATHPRIMITIVE(FLOOR,floor)
MAKEMATHPRIMITIVE(CEIL,ceil)
MAKEMATHPRIMITIVE(ABS,fabs)
MAKEMATHPRIMITIVE2(ATAN2,atan2)
MAKEMATHPRIMITIVE2(HYPOT,hypot)

MAKEMATHPRIMITIVE(ERF,erf)
MAKEMATHPRIMITIVE(GAMMA,tgamma)


static Obj FuncSIGN_MACFLOAT(Obj self, Obj f)
{
  Double vf = VAL_MACFLOAT(f);

  return vf == 0. ? INTOBJ_INT(0) : signbit(vf) ? INTOBJ_INT(-1) : INTOBJ_INT(1);
}

static Obj FuncSIGNBIT_MACFLOAT(Obj self, Obj f)
{
  return signbit(VAL_MACFLOAT(f)) ? True : False;
}


static Obj FuncINTFLOOR_MACFLOAT(Obj self, Obj macfloat)
{
    RequireMacFloat(SELF_NAME, macfloat);

    Double f = VAL_MACFLOAT(macfloat);
    if (isnan(f))
        ErrorQuit("cannot convert float nan to integer", 0, 0);
    if (isinf(f))
        ErrorQuit("cannot convert float %s to integer", (Int)(f > 0 ? "inf" : "-inf"), 0);

#ifdef HAVE_TRUNC
  f = trunc(f);
#else
  if (f >= 0.0)
    f = floor(f);
  else
    f = -floor(-f);
#endif


  if (fabs(f) < (Double)((Int)1 << NR_SMALL_INT_BITS))
    return INTOBJ_INT((Int)f);

  int str_len = (int) (log(fabs(f)) / log(16.0)) + 3;

  Obj str = NEW_STRING(str_len);
  char *s = CSTR_STRING(str), *p = s+str_len-1;
  if (f < 0.0) {
    f = -f;
    s[0] = '-';
  }
  while (p > s || (p == s && s[0] != '-')) {
    int d = (int) fmod(f,16.0);
    *p-- = d < 10 ? '0'+d : 'a'+d-10;
    f /= 16.0;
  }
  return IntHexString(str);
}

static Obj FuncSTRING_DIGITS_MACFLOAT(Obj self, Obj gapprec, Obj f)
{
  Char buf[1024];
  Obj str;
  int prec = INT_INTOBJ(gapprec);
  if (prec > 40) // too much anyways, and would risk buffer overrun
    prec = 40;
  PrintMacfloatToBuf(buf, sizeof(buf), VAL_MACFLOAT(f), prec);
  str = MakeString(buf);
  return str;
}

static Obj FuncLDEXP_MACFLOAT(Obj self, Obj f, Obj i)
{
  return NEW_MACFLOAT(ldexp(VAL_MACFLOAT(f),INT_INTOBJ(i)));
}

static Obj FuncFREXP_MACFLOAT(Obj self, Obj f)
{
  int i;
  Obj d = NEW_MACFLOAT(frexp (VAL_MACFLOAT(f), &i));
  Obj l = NEW_PLIST(T_PLIST,2);
  SET_ELM_PLIST(l,1,d);
  SET_ELM_PLIST(l,2,INTOBJ_INT(i));
  SET_LEN_PLIST(l,2);
  return l;
}

/****************************************************************************
**
*F * * * * * * * * * * * * * initialize module * * * * * * * * * * * * * * *
*/


/****************************************************************************
**
*V  BagNames  . . . . . . . . . . . . . . . . . . . . . . . list of bag names
*/
static StructBagNames BagNames[] = {
  { T_MACFLOAT, "macfloat" },
  { -1,    "" }
};


/****************************************************************************
**
*V  GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
*/
static StructGVarFunc GVarFuncs [] = {
    GVAR_FUNC_1ARGS(MACFLOAT_INT, int),
    GVAR_FUNC_1ARGS(MACFLOAT_STRING, string),

    GVAR_FUNC_1ARGS(SIN_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(COS_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(TAN_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ASIN_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ACOS_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ATAN_MACFLOAT, macfloat),

    GVAR_FUNC_1ARGS(SINH_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(COSH_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(TANH_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ASINH_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ACOSH_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ATANH_MACFLOAT, macfloat),

    GVAR_FUNC_2ARGS(ATAN2_MACFLOAT, real, imag),
    GVAR_FUNC_2ARGS(HYPOT_MACFLOAT, real, imag),

    GVAR_FUNC_1ARGS(LOG_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(LOG2_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(LOG10_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(LOG1P_MACFLOAT, macfloat),

    GVAR_FUNC_1ARGS(EXP_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(EXP2_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(EXPM1_MACFLOAT, macfloat),
#ifdef HAVE_EXP10
    GVAR_FUNC_1ARGS(EXP10_MACFLOAT, macfloat),
#endif

    GVAR_FUNC_2ARGS(LDEXP_MACFLOAT, macfloat, int),
    GVAR_FUNC_1ARGS(FREXP_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(SQRT_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(CBRT_MACFLOAT, macfloat),

    GVAR_FUNC_1ARGS(ERF_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(GAMMA_MACFLOAT, macfloat),

    GVAR_FUNC_1ARGS(RINT_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(INTFLOOR_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(FLOOR_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(CEIL_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(ABS_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(SIGN_MACFLOAT, macfloat),
    GVAR_FUNC_1ARGS(SIGNBIT_MACFLOAT, macfloat),

    GVAR_FUNC_2ARGS(STRING_DIGITS_MACFLOAT, digits, macfloat),
    GVAR_FUNC_2ARGS(EQ_MACFLOAT, x, y),
    { 0, 0, 0, 0, 0 }
};


/****************************************************************************
**
*F  InitKernel( <module> )  . . . . . . . . initialise kernel data structures
*/
static Int InitKernel (
    StructInitInfo *    module )
{
    // set the bag type names (for error messages and debugging)
    InitBagNamesFromTable( BagNames );

    // install the marking functions for macfloatean values
    InitMarkFuncBags( T_MACFLOAT, MarkNoSubBags );
#ifdef HPCGAP
    MakeBagTypePublic( T_MACFLOAT );
#endif

    // init functions
    InitHdlrFuncsFromTable( GVarFuncs );

    // install the type function
    ImportGVarFromLibrary( "TYPE_MACFLOAT", &TYPE_MACFLOAT );
    TypeObjFuncs[ T_MACFLOAT ] = TypeMacfloat;

#ifdef GAP_ENABLE_SAVELOAD
    // install the saving functions
    SaveObjFuncs[ T_MACFLOAT ] = SaveMacfloat;

    // install the loading functions
    LoadObjFuncs[ T_MACFLOAT ] = LoadMacfloat;
#endif

    // install the printer for macfloatean values
    PrintObjFuncs[ T_MACFLOAT ] = PrintMacfloat;

    // install the comparison functions
    EqFuncs[ T_MACFLOAT ][ T_MACFLOAT ] = EqMacfloat;
    LtFuncs[ T_MACFLOAT ][ T_MACFLOAT ] = LtMacfloat;

    // allow method selection to protest against comparisons of float and int
    for (int t = T_INT; t <= T_CYC; t++)
        EqFuncs[T_MACFLOAT][t] = EqFuncs[t][T_MACFLOAT] = EqObject;

    // install the unary arithmetic methods
    ZeroSameMutFuncs[T_MACFLOAT] = ZeroMacfloat;
    ZeroMutFuncs[ T_MACFLOAT ] = ZeroMacfloat;
    AInvMutFuncs[ T_MACFLOAT ] = AInvMacfloat;
    OneFuncs [ T_MACFLOAT ] = OneMacfloat;
    OneSameMut[T_MACFLOAT] = OneMacfloat;
    InvFuncs [ T_MACFLOAT ] = InvMacfloat;

    // install binary arithmetic methods
    ProdFuncs[ T_MACFLOAT ][ T_MACFLOAT ] = ProdMacfloat;
    PowFuncs [ T_MACFLOAT ][ T_MACFLOAT ] = PowMacfloat;
    SumFuncs[ T_MACFLOAT ][ T_MACFLOAT ] = SumMacfloat;
    DiffFuncs [ T_MACFLOAT ][ T_MACFLOAT ] = DiffMacfloat;
    QuoFuncs [ T_MACFLOAT ][ T_MACFLOAT ] = QuoMacfloat;
    LQuoFuncs [ T_MACFLOAT ][ T_MACFLOAT ] = LQuoMacfloat;
    ModFuncs [ T_MACFLOAT ][ T_MACFLOAT ] = ModMacfloat;
    SumFuncs [ T_INT ][ T_MACFLOAT ] = SumIntMacfloat;

    // Probably support mixed ops with small ints in the kernel as well
    // on any reasonable system, all small ints should have macfloat
    // equivalents

    // Anything else, like mixed ops with rationals, we can leave to the
    // library at least for a while

    return 0;
}


/****************************************************************************
**
*F  InitLibrary( <module> ) . . . . . . .  initialise library data structures
*/
static Int InitLibrary (
    StructInitInfo *    module )
{
    // init functions
    InitGVarFuncsFromTable( GVarFuncs );

    return 0;
}


/****************************************************************************
**
*F  InitInfoMacfloat()  . . . . . . . . . . . . . . . . . table of init functions
*/
static StructInitInfo module = {
    // init struct using C99 designated initializers; for a full list of
    // fields, please refer to the definition of StructInitInfo
    .type = MODULE_BUILTIN,
    .name = "macfloat",
    .initKernel = InitKernel,
    .initLibrary = InitLibrary,
};

StructInitInfo * InitInfoMacfloat ( void )
{
    return &module;
}