#ifndef __PLINK2_FLOAT_H__
#define __PLINK2_FLOAT_H__

// This library is part of PLINK 2.0, copyright (C) 2005-2026 Shaun Purcell,
// Christopher Chang.
//
// This program is free software: you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published by the
// Free Software Foundation, either version 3 of the License, or (at your
// option) any later version.
//
// This library 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 Lesser General Public License
// for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.


// Basic floating-point constants and portable definitions.

#include <math.h>

#include "plink2_base.h"

// apparently these aren't always defined in limits.h
#ifndef DBL_MAX
#  define DBL_MAX 1.7976931348623157e308
#endif
#ifndef DBL_MIN
#  define DBL_MIN 2.2250738585072014e-308
#endif
#ifndef FLT_MAX
#  define FLT_MAX S_CAST(float, 3.40282347e38)
#endif

#ifdef __cplusplus
namespace plink2 {
#endif

static const double FLT_MAX_D = 3.4028235677973362e38;
// INFINITY is usually a float, but apparently not in w64-mingw32-g++.
static const float INFINITY_F = S_CAST(float, INFINITY);
static const double INFINITY_D = S_CAST(double, INFINITY);

// Notes on float64 ('double') precision:
// * 1 + 2^{-52} is the smallest float64 greater than 1, and 1 - 2^{-53} is the
//   largest float64 less than 1.  Since plink2 flushes denormals to zero on
//   the main target platforms, a 'ULP' ("unit in the last place") can be
//   assumed to be between 2^{-52} and 2^{-53} times the total value.
//   When manually updating an epsilon value that's ever added to 1, multiples
//   of 2^{-52} should be used.  If epsilon is only ever subtracted from 1, it
//   is ok to use multiples of 2^{-53}.
// * With default rounding behavior, the maximum rounding error from a basic
//   arithmetic operation is 0.5 ULP, i.e. ~2^{-53} times the result.
//   plink2 functions are written to assume default rounding, and to avoid
//   assumptions about denormal handling.
// * Unfortunately, accuracy of functions like exp and log is currently
//   platform-dependent, and it isn't realistically worthwhile to purge that
//   source of platform variation when we're continuing to use different linear
//   algebra libraries on different platforms.
//   The best survey I've found of the current state of the world is
//     Gladman B, Innocente V, Mather J, Zimmermann P, (2024) Accuracy of
//     Mathematical Functions in Single, Double, Double Extended, and Quadruple
//     Precision.  hal-03141101v7.
//     https://inria.hal.science/hal-03141101v7/document
//   It reports the following maximum observed errors (not an exhaustive check
//   due to the size of the search space):
//     acos: 1.53 ULP (CUDA 12.2.1)
//     cbrt: 1.53e22 ULP (!) (AMD 4.2); next-worst is GNU libc 2.40 at 3.67 ULP
//           AMD failure mode only applies to denormals.
//     cos/sin: Inf (!) (LLVM 18.1.8); next-worst is CUDA 12.2.1 at 1.52 ULP
//              LLVM failure mode only applies to very large arguments.
//     exp: 1.5 ULP (MSVC 2022)
//     expm1: 3.06 ULP (MSVC 2022)
//     log: 0.946 ULP (Newlib 4.4.0, OpenLibm 0.8.3, FreeBSD 14.1)
//     log1p: 1.74 ULP (ArmPL 24.04)
//     sqrt: 0.5 ULP
static const double k2m24 = 1.0 / (1LL << 24);
static const double k2m25 = 1.0 / (1LL << 25);
static const double k2m32 = 1.0 / (1LL << 32);
static const double k2m50 = 1.0 / (1LL << 50);
static const double k2m52 = 1.0 / (1LL << 52);
static const double k2m53 = 1.0 / (1LL << 53);
static const double k2m54 = 1.0 / (1LL << 54);
static const double k2m64 = k2m54 / (1 << 10);
static const double k2p50 = 1.0 * (1LL << 50);
static const double k2p64 = 4.0 * (1LL << 62);
static const double k2p100 = k2p50 * k2p50;
static const double k2p200 = k2p100 * k2p100;
static const double k2p400 = k2p200 * k2p200;
static const double k2p800 = k2p400 * k2p400;
static const double kE = 2.7182818284590452;
static const double kLn2 = 0.6931471805599453;
static const double kLn10 = 2.3025850929940457;
static const double kLnNormalMin = -708.3964185322641;
static const double kLnSqrtPi = 0.5723649429247001;
static const double kLnSqrt2Pi = 0.91893853320467278056;
static const double kPi = 3.1415926535897932;
static const double k2Pi = 6.283185307179586;
static const double kRecipE = 0.36787944117144233;
static const double kRecipLn10 = 0.43429448190325176;
static const double kSqrt2 = 1.4142135623730951;
static const double kSqrtPi = 1.7724538509055159;

// Some more negative powers of 2, mostly used as tolerances for floating-point
// approximate-equality checks.
static const double k2m21 = 1.0 / (1 << 21);
static const double k2m30 = 1.0 / (1 << 30);
static const double k2m35 = 1.0 / (1LL << 35);
static const double k2m44 = 1.0 / (1LL << 44);
static const double k2m60 = 1.0 / (1LL << 60);
// Square of this value underflows float64 (even without denormal flushing).
static const double k2m537p5 = kSqrt2 / (k2p400 * k2p100 * (1LL << 38));

static const double kBigEpsilon = k2m21;  // must be >= sqrt(kSmallEpsilon)
static const double kEpsilon = k2m30;
static const double kSmallEpsilon = k2m44;

#if defined(__cplusplus)
#  if __cplusplus >= 201103L
HEADER_INLINE bool isfinite_f(float fxx) {
  using namespace std;
  return isfinite(fxx);
}

HEADER_INLINE bool isfinite_d(double dxx) {
  using namespace std;
  return isfinite(dxx);
}
#  else
#    ifdef isfinite
#      define isfinite_f isfinite
#      define isfinite_d isfinite
#    else
HEADER_INLINE bool isfinite_f(float fxx) {
  return (fxx == fxx) && (fxx != INFINITY_F) && (fxx != -INFINITY_F);
}

HEADER_INLINE bool isfinite_d(double dxx) {
  return (dxx == dxx) && (dxx != INFINITY_D)) && (dxx != -INFINITY_D));
}
#    endif
#  endif
#else
#  define isfinite_f isfinite
#  define isfinite_d isfinite
#endif

// Floating-point environment setup.
//
// - plink2 is designed to work properly when denormals are flushed to zero, so
//   setting that flag is an easy performance win.
// - On 32-bit systems, the implementation of plink2_highprec's double-double
//   type assumes intermediate results are rounded by the processor to 64 bits,
//   not 80.  ...honestly it's fine to just drop 32-bit support at this point,
//   but passing "-msse -mfpmath=sse" (which doesn't work with chips older than
//   ~2000) to gcc is a simple fix.

void flush_denormals();

// We treat minor floating-point differences across platforms as inevitable,
// since we don't want to lock ourselves to a primitive linear algebra library
// implementation.  But we try to avoid introducing such differences without
// good reason.
//
// There is a fused-multiply-add (FMA) instruction on ARM and x86_64 v3 (AVX2),
// but not older x86.  This can speed up some critical loops by ~30% and/or
// improve accuracy when used properly.  That can be enough benefit to justify
// the cost of more minor differences between x86_64 v1 and v3 floating-point
// results.  The C language standard allows compilers to automatically
// "contract" (x*y)+z in a single expression to fma(x,y,z) when the target has
// efficient hardware support for the latter, and both clang 14+ and gcc were
// doing this under some of our build configurations.
//
// Unfortunately, this can break correctness of code like that in
// plink2_highprec, and even the latest gcc does not honor the #pragma that is
// supposed to be able to shield such code.  That leaves two options:
// 1. Compile files like include/plink2_highprec.cc with -ffp-contract=off, and
//    allow contraction elsewhere.
// 2. Compile everything with -ffp-contract=off, and then insert allowed-FMAs
//    by hand.
// I've decided to go with option 2 since I'm generally trying to innovate at
// least a little bit on floating-point accuracy, the amount of time-critical
// floating-point code to review is relatively small, and my preexisting code
// was not written with much awareness of how automatic contractions could
// endanger accuracy.
// (Simple counterintuitive example: (x1*y1)-(x2*y2) can be contracted to
// fma(x1, y1, -x2*y2), potentially yielding a nonzero result when x1=x2 and
// y1=y2.  See also
//   https://stackoverflow.com/questions/73985098/clang-14-0-0-floating-point-optimizations
// .)
#if defined (__FMA__) || defined(__ARM_FEATURE_FMA)
#  define USE_FMA
// Note that this typically isn't worthwhile when multiplying by a power-of-2
// constant (doesn't require a regular multiply op), and that careful attention
// must be paid to dependency chains (if third argument isn't naturally
// available by the time the first two are, FMA is likely to slow things down;
// and in general, using multiple accumulators is a lot more helpful than
// inserting FP contractions).
HEADER_INLINE double prefer_fma(double a, double b, double c) {
  return fma(a, b, c);
}

HEADER_INLINE float prefer_fmaf(float a, float b, float c) {
  return fmaf(a, b, c);
}
#else
HEADER_INLINE double prefer_fma(double a, double b, double c) {
  return a * b + c;
}

HEADER_INLINE float prefer_fmaf(float a, float b, float c) {
  return a * b + c;
}
#endif

// General-purpose polynomial and rational-function evaluators, degree doesn't
// need to be known at compile time.
double poly_eval(const double* coefs, uint32_t degree, double xx);

double ratfun_eval_smallx(const double* numer_coefs, const double* denom_coefs, uint32_t degree, double xx);

double ratfun_eval_largex(const double* numer_coefs, const double* denom_coefs, uint32_t degree, double xx);

HEADER_INLINE double ratfun_eval(const double* numer_coefs, const double* denom_coefs, uint32_t degree, double xx) {
  if (xx <= 1) {
    return ratfun_eval_smallx(numer_coefs, denom_coefs, degree, xx);
  }
  return ratfun_eval_largex(numer_coefs, denom_coefs, degree, xx);
}

// Hardcoded sequences of prefer_fma() instructions for
// polynomial/rational-function  evaluation, for use when degree is known at
// compile-time.  (todo: check whether the largest hardcoded degree here still
// tends to be inlined by modern compilers; when that isn't the case,
// poly_eval/ratfun_eval should be used instead.)

// Common case is evaluation of a rational function, with both numerator and
// denominator polynomials.  Two FMA chains per polynomial should work well
// across a variety of platforms.
//
// Occasionally it is useful to call this form directly (see plink2_stats
// QuantileToZscore()).
HEADER_INLINE double _poly4(double xx, double x2, double c0, double c1, double c2, double c3, double c4) {
  const double even_terms =
    prefer_fma(
      prefer_fma(
        c4,
        x2,
        c2),
      x2,
      c0);
  const double odd_terms =
    prefer_fma(
      c3,
      x2,
      c1)
    * xx;
  return even_terms + odd_terms;
}

HEADER_INLINE double _poly5(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5) {
  return _poly4(xx, x2, c0, c1, c2, prefer_fma(c5, x2, c3), c4);
}

HEADER_INLINE double _poly6(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6) {
  return _poly5(xx, x2, c0, c1, c2, c3, prefer_fma(c6, x2, c4), c5);
}

HEADER_INLINE double _poly7(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7) {
  return _poly6(xx, x2, c0, c1, c2, c3, c4, prefer_fma(c7, x2, c5), c6);
}

HEADER_INLINE double _poly8(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {
  return _poly7(xx, x2, c0, c1, c2, c3, c4, c5, prefer_fma(c8, x2, c6), c7);
}

HEADER_INLINE double _poly9(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9) {
  return _poly8(xx, x2, c0, c1, c2, c3, c4, c5, c6, prefer_fma(c9, x2, c7), c8);
}

HEADER_INLINE double _poly10(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10) {
  return _poly9(xx, x2, c0, c1, c2, c3, c4, c5, c6, c7, prefer_fma(c10, x2, c8), c9);
}

HEADER_INLINE double _poly11(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10, double c11) {
  return _poly10(xx, x2, c0, c1, c2, c3, c4, c5, c6, c7, c8, prefer_fma(c11, x2, c9), c10);
}

HEADER_INLINE double _poly12(double xx, double x2, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10, double c11, double c12) {
  return _poly11(xx, x2, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, prefer_fma(c12, x2, c10), c11);
}

// In the relatively rare case where we only need to evaluate one polynomial
// and it has degree 8+, go up to 4 hardcoded FMA chains.
HEADER_INLINE double _poly8_solo(double xx, double x2, double x3, double x4, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {
  const double rem0 =
    prefer_fma(
      prefer_fma(
        c8,
        x4,
        c4),
      x4,
      c0);
  const double rem1 =
    prefer_fma(
      c5,
      x4,
      c1)
    * xx;
  const double rem2 =
    prefer_fma(
      c6,
      x4,
      c2)
    * x2;
  const double rem3 =
    prefer_fma(
      c7,
      x4,
      c3)
    * x3;
  return rem0 + rem1 + (rem2 + rem3);
}

HEADER_INLINE double _poly9_solo(double xx, double x2, double x3, double x4, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9) {
  return _poly8_solo(xx, x2, x3, x4, c0, c1, c2, c3, c4, prefer_fma(c9, x4, c5), c6, c7, c8);
}

HEADER_INLINE double POLY2(double xx, double c0, double c1, double c2) {
  return
    prefer_fma(
      prefer_fma(
        c2,
        xx,
        c1),
      xx,
      c0);
}

HEADER_INLINE double POLY3(double xx, double c0, double c1, double c2, double c3) {
  return
    prefer_fma(
      prefer_fma(
        prefer_fma(
          c3,
          xx,
          c2),
        xx,
        c1),
      xx,
      c0);
}

HEADER_INLINE double POLY4(double xx, double c0, double c1, double c2, double c3, double c4) {
  return _poly4(xx, xx * xx, c0, c1, c2, c3, c4);
}

HEADER_INLINE double POLY5(double xx, double c0, double c1, double c2, double c3, double c4, double c5) {
  return _poly5(xx, xx * xx, c0, c1, c2, c3, c4, c5);
}

HEADER_INLINE double POLY6(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6) {
  return _poly6(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6);
}

HEADER_INLINE double POLY7(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7) {
  return _poly7(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7);
}

HEADER_INLINE double POLY8(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {
  return _poly8(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7, c8);
}

HEADER_INLINE double POLY9(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9) {
  return _poly9(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9);
}

HEADER_INLINE double POLY10(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10) {
  return _poly10(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10);
}

HEADER_INLINE double POLY11(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10, double c11) {
  return _poly11(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11);
}

HEADER_INLINE double POLY12(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10, double c11, double c12) {
  return _poly12(xx, xx * xx, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12);
}


HEADER_INLINE double POLY8_SOLO(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {
  const double x2 = xx * xx;
  return _poly8_solo(xx, x2, xx * x2, x2 * x2, c0, c1, c2, c3, c4, c5, c6, c7, c8);
}

HEADER_INLINE double POLY9_SOLO(double xx, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9) {
  const double x2 = xx * xx;
  return _poly9_solo(xx, x2, xx * x2, x2 * x2, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9);
}



HEADER_INLINE uint64_t float64bits(double xx) {
  double* xx_ptr = &xx;
  return *R_CAST(uint64_t*, xx_ptr);
}

HEADER_INLINE double float64frombits(uint64_t ullii) {
  uint64_t* ullii_ptr = &ullii;
  return *R_CAST(double*, ullii_ptr);
}

HEADER_INLINE double prev_float64(double xx) {
  return float64frombits(float64bits(xx) - 1);
}

HEADER_INLINE double next_float64(double xx) {
  return float64frombits(float64bits(xx) + 1);
}

// Returns log(exp(xx) + exp(yy)).
HEADER_INLINE double lnsum(double xx, double yy) {
  // log(exp(xx) + exp(yy))
  // = log(exp(max(xx,yy)) * (1 + exp(min(xx,yy) - max(xx,yy))))
  // = max(xx,yy) + log(1 + exp(min(xx,yy) - max(xx,yy)))
  // If max(xx,yy) >> min(xx,yy), quickly return max(xx,yy).
  // If xx and yy are close enough to each other that exp(min(xx,yy) -
  // max(xx,yy)) rounds to 1, nothing bad happens; we can focus on maximizing
  // accuracy in the exp(min(xx,yy) - max(xx,yy)) near zero case.
  const double max_arg = MAXV(xx, yy);
  const double ln_ratio = MINV(xx, yy) - max_arg;
  if (ln_ratio < -54 * kLn2) {
    return max_arg;
  }
  return max_arg + log1p(exp(ln_ratio));
}

// Returns log(exp(xx) - exp(yy)), assuming xx > yy.
HEADER_INLINE double lndiff(double xx, double yy) {
  // log(exp(xx) - exp(yy))
  // = log(exp(xx) * (1 - exp(yy-xx)))
  // = xx + log(1 - exp(yy-xx))
  // Subcases of interest:
  // * xx >> yy: exp(yy-xx) indistinguishable from zero, quickly return xx.
  // * (yy-xx) very small: if exp(yy-xx) rounds to 1, we have a problem.
  //   xx + log(-expm1(yy-xx)) avoids the problem.
  const double ln_ratio = yy - xx;
  if (ln_ratio < -54 * kLn2) {
    // yy is too small to matter.
    return xx;
  }
  return xx + log(-expm1(ln_ratio));
}

HEADER_INLINE double ceil_limit(double xx, double limit) {
  if (xx >= limit) {
    return limit;
  }
  return ceil(xx);
}

HEADER_INLINE double flush_if_denormal(double xx) {
  if (fabs(xx) < DBL_MIN) {
    return 0.0;
  }
  return xx;
}

HEADER_INLINE double exp_flush(double xx) {
  return flush_if_denormal(exp(xx));
}

#ifdef __cplusplus
}  // namespace plink2
#endif

#endif  // __PLINK2_FLOAT_H__
