// SPDX-License-Identifier: MPL-2.0
// (c) Hare authors <https://harelang.org>

// Sections of the code below, in particular log() and exp(), are based on Go's
// implementation, which is, in turn, based on FreeBSD's. The original C code,
// as well as the respective comments and constants are from
// /usr/src/lib/msun/src/{e_log,e_exp}.c.
//
// The FreeBSD copyright notice:
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
//
// Developed at SunPro, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.
// ====================================================
//
// The Go copyright notice:
// ====================================================
// Copyright (c) 2009 The Go Authors. All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
//    * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//    * Redistributions in binary form must reproduce the above
// copyright notice, this list of conditions and the following disclaimer
// in the documentation and/or other materials provided with the
// distribution.
//    * Neither the name of Google Inc. nor the names of its
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
// ====================================================

use types;

// The standard tolerance used by [[isclosef32]] and [[isclosef64]], which is
// just an arbitrary way to measure whether two floating-point numbers are
// "sufficiently" close to each other.
export def STANDARD_TOL = 1e-14;

// Returns whether x and y are within tol of each other.
export fn isclosef64(x: f64, y: f64, tol: f64 = STANDARD_TOL) bool = {
	if (isnan(x) || isnan(y) || isnan(tol)) {
		return false;
	};
	return absf64(x - y) < tol;
};

// Returns whether x and y are within tol of each other.
export fn isclosef32(x: f32, y: f32, tol: f32 = STANDARD_TOL) bool = {
	if (isnan(x) || isnan(y) || isnan(tol)) {
		return false;
	};
	return absf32(x - y) < tol;
};

// e - https://oeis.org/A001113
export def E: f64 = 2.71828182845904523536028747135266249775724709369995957496696763;
// pi - https://oeis.org/A000796
export def PI: f64 = 3.14159265358979323846264338327950288419716939937510582097494459;
// tau - https://oeis.org/A019692
export def TAU: f64 = 6.2831853071795864769252867665590057683943387987502116419498892;
// phi - https://oeis.org/A001622
export def PHI: f64 = 1.61803398874989484820458683436563811772030917980576286213544862;
// sqrt(2) - https://oeis.org/A002193
export def SQRT_2: f64 = 1.41421356237309504880168872420969807856967187537694807317667974;
// sqrt(e) - https://oeis.org/A019774
export def SQRT_E: f64 = 1.64872127070012814684865078781416357165377610071014801157507931;
// sqrt(pi) - https://oeis.org/A002161
export def SQRT_PI: f64 = 1.77245385090551602729816748334114518279754945612238712821380779;
// sqrt(phi) - https://oeis.org/A139339
export def SQRT_PHI: f64 = 1.27201964951406896425242246173749149171560804184009624861664038;
// ln(2) - https://oeis.org/A002162
export def LN_2: f64 = 0.693147180559945309417232121458176568075500134360255254120680009;
// ln(2) - https://oeis.org/A002162
export def LN2_HI: f64 = 6.93147180369123816490e-01;
// ln(2) - https://oeis.org/A002162
export def LN2_LO: f64 = 1.90821492927058770002e-10;
// log_{2}(e)
export def LOG2_E: f64 = 1f64 / LN_2;
// ln(10) - https://oeis.org/A002392
export def LN_10: f64 = 2.30258509299404568401799145468436420760110148862877297603332790;
// log_{10}(e)
export def LOG10_E: f64 = 1f64 / LN_10;

// __ieee754_log(x)
// Return the logarithm of x
//
// Method :
//   1. Argument Reduction: find k and f such that
//			x = 2**k * (1+f),
//	   where  sqrt(2)/2 < 1+f < sqrt(2) .
//
//   2. Approximation of log(1+f).
//	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
//		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
//	     	 = 2s + s*R
//      We use a special Reme algorithm on [0,0.1716] to generate
//	a polynomial of degree 14 to approximate R.  The maximum error
//	of this polynomial approximation is bounded by 2**-58.45. In
//	other words,
//		        2      4      6      8      10      12      14
//	    R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s  +L6*s  +L7*s
//	(the values of L1 to L7 are listed in the program) and
//	    |      2          14          |     -58.45
//	    | L1*s +...+L7*s    -  R(z) | <= 2
//	    |                             |
//	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
//	In order to guarantee error in log below 1ulp, we compute log by
//		log(1+f) = f - s*(f - R)		(if f is not too large)
//		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
//
//	3. Finally,  log(x) = k*Ln2 + log(1+f).
//			    = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo)))
//	   Here Ln2 is split into two floating point number:
//			Ln2_hi + Ln2_lo,
//	   where n*Ln2_hi is always exact for |n| < 2000.
//
// Special cases:
//	log(x) is NaN with signal if x < 0 (including -INF) ;
//	log(+INF) is +INF; log(0) is -INF with signal;
//	log(NaN) is that NaN with no signal.
//
// Accuracy:
//	according to an error analysis, the error is always less than
//	1 ulp (unit in the last place).
//
// Constants:
// The hexadecimal values are the intended ones for the following
// constants. The decimal values may be used, provided that the
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.

// Returns the natural logarithm of x.
export fn logf64(x: f64) f64 = {
	const L1 = 6.666666666666735130e-01; // 3fe55555 55555593
	const L2 = 3.999999999940941908e-01; // 3fd99999 9997fa04
	const L3 = 2.857142874366239149e-01; // 3fd24924 94229359
	const L4 = 2.222219843214978396e-01; // 3fcc71c5 1d8e78af
	const L5 = 1.818357216161805012e-01; // 3fc74664 96cb03de
	const L6 = 1.531383769920937332e-01; // 3fc39a09 d078c69f
	const L7 = 1.479819860511658591e-01; // 3fc2f112 df3e5244

	// special cases
	if (isnan(x) || x == INF) {
		return x;
	} else if (x < 0f64) {
		return NAN;
	} else if (x == 0f64) {
		return -INF;
	};

	// Reduce
	const (f1, ki) = frexpf64(x);
	if (f1 < (SQRT_2 / 2f64)) {
		f1 *= 2f64;
		ki -= 1i64;
	};
	let f = f1 - 1f64;
	let k = (ki: f64);

	// Compute
	const s = f / (2f64 + f);
	const s2 = s * s;
	const s4 = s2 * s2;
	const t1 = s2 * (L1 + s4 * (L3 + s4 * (L5 + s4 * L7)));
	const t2 = s4 * (L2 + s4 * (L4 + s4 * L6));
	const R = t1 + t2;
	const hfsq = 0.5f64 * f * f;
	return k * LN2_HI - ((hfsq - (s * (hfsq + R) + k * LN2_LO)) - f);
};

// Returns the decimal logarithm of x.
export fn log10f64(x: f64) f64 = {
	return logf64(x) * (1f64 / LN_10);
};

// Returns the binary logarithm of x.
export fn log2f64(x: f64) f64 = {
	const (frac, exp) = frexpf64(x);
	// Make sure exact powers of two give an exact answer.
	// Don't depend on log(0.5) * (1 / LN_2) + exp being exactly exp - 1.
	if (frac == 0.5f64) {
		return ((exp - 1): f64);
	};
	return logf64(frac) * (1f64 / LN_2) + (exp: f64);
};

// double log1p(double x)
//
// Method :
//   1. Argument Reduction: find k and f such that
//                      1+x = 2**k * (1+f),
//         where  sqrt(2)/2 < 1+f < sqrt(2) .
//
//      Note. If k=0, then f=x is exact. However, if k!=0, then f
//      may not be representable exactly. In that case, a correction
//      term is need. Let u=1+x rounded. Let c = (1+x)-u, then
//      log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
//      and add back the correction term c/u.
//      (Note: when x > 2**53, one can simply return log(x))
//
//   2. Approximation of log1p(f).
//      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
//               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
//               = 2s + s*R
//      We use a special Reme algorithm on [0,0.1716] to generate
//      a polynomial of degree 14 to approximate R The maximum error
//      of this polynomial approximation is bounded by 2**-58.45. In
//      other words,
//                      2      4      6      8      10      12      14
//          R(z) ~ LP1*s +LP2*s +LP3*s +LP4*s +LP5*s  +LP6*s  +LP7*s
//      (the values of LP1 to LP7 are listed in the program)
//      and
//          |      2          14          |     -58.45
//          | LP1*s +...+LP7*s    -  R(z) | <= 2
//          |                             |
//      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
//      In order to guarantee error in log below 1ulp, we compute log
//      by
//              log1p(f) = f - (hfsq - s*(hfsq+R)).
//
//   3. Finally, log1p(x) = k*ln2 + log1p(f).
//                        = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
//      Here ln2 is split into two floating point number:
//                   ln2_hi + ln2_lo,
//      where n*ln2_hi is always exact for |n| < 2000.
//
// Special cases:
//      log1p(x) is NaN with signal if x < -1 (including -INF) ;
//      log1p(+INF) is +INF; log1p(-1) is -INF with signal;
//      log1p(NaN) is that NaN with no signal.
//
// Accuracy:
//      according to an error analysis, the error is always less than
//      1 ulp (unit in the last place).
//
// Constants:
// The hexadecimal values are the intended ones for the following
// constants. The decimal values may be used, provided that the
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.
//
// Note: Assuming log() return accurate answer, the following
//       algorithm can be used to compute log1p(x) to within a few ULP:
//
//              u = 1+x;
//              if(u==1.0) return x ; else
//                         return log(u)*(x/(u-1.0));
//
//       See HP-15C Advanced Functions Handbook, p.193.

// Returns the natural logarithm of 1 plus its argument x.
// It is more accurate than log(1 + x) when x is near zero.
export fn log1pf64(x: f64) f64 = {
	// sqrt(2) - 1
	const SQRT2M1 = 4.142135623730950488017e-01; // 0x3fda827999fcef34
	// sqrt(2) / 2 - 1
	const SQRT2HALFM1 = -2.928932188134524755992e-01; // 0xbfd2bec333018866
	const SMALL = 1f64 / ((1i64 << 29): f64); // 2**-29
	const TINY = 1f64 / ((1i64 << 54): f64); // 2**-54
	const TWO53 = ((1i64 << 53): f64); // 2**53
	const LN2HI = 6.93147180369123816490e-01; // 3fe62e42fee00000
	const LN2LO = 1.90821492927058770002e-10; // 3dea39ef35793c76
	const LP1 = 6.666666666666735130e-01; // 3fe5555555555593
	const LP2 = 3.999999999940941908e-01; // 3fd999999997fa04
	const LP3 = 2.857142874366239149e-01; // 3fd2492494229359
	const LP4 = 2.222219843214978396e-01; // 3fcc71c51d8e78af
	const LP5 = 1.818357216161805012e-01; // 3fc7466496cb03de
	const LP6 = 1.531383769920937332e-01; // 3fc39a09d078c69f
	const LP7 = 1.479819860511658591e-01; // 3fc2f112df3e5244

	if (x < -1f64 || isnan(x)) {
		return NAN;
	} else if (x == -1f64) {
		return -INF;
	} else if (x == INF) {
		return INF;
	};

	const absx = absf64(x);

	let f = 0f64;
	let iu = 0u64;
	let k = 1i64;
	if (absx < SQRT2M1) { //  |x| < Sqrt(2)-1
		if (absx < SMALL) { // |x| < 2**-29
			if (absx < TINY) { // |x| < 2**-54
				return x;
			};
			return x - (x * x * 0.5f64);
		};
		if (x > SQRT2HALFM1) { // Sqrt(2)/2-1 < x
			// (Sqrt(2)/2-1) < x < (Sqrt(2)-1)
			k = 0;
			f = x;
			iu = 1;
		};
	};
	let c = 0f64;
	if (k != 0) {
		let u = 0f64;
		if (absx < TWO53) { // 1<<53
			u = 1.0 + x;
			iu = f64bits(u);
			k = (((iu >> 52) - 1023): i64);
			// Correction term
			if (k > 0) {
				c = 1f64 - (u - x);
			} else {
				c = x - (u - 1f64);
			};
			c /= u;
		} else {
			u = x;
			iu = f64bits(u);
			k = (((iu >> 52) - 1023): i64);
			c = 0f64;
		};
		iu &= 0x000fffffffffffff;
		if (iu < 0x0006a09e667f3bcd) { // Mantissa of Sqrt(2)
			// Normalize u
			u = f64frombits(iu | 0x3ff0000000000000);
		} else {
			k += 1;
			// Normalize u/2
			u = f64frombits(iu | 0x3fe0000000000000);
			iu = (0x0010000000000000 - iu) >> 2;
		};
		f = u - 1f64; // Sqrt(2)/2 < u < Sqrt(2)
	};
	const hfsq = 0.5 * f * f;
	let s = 0f64;
	let R = 0f64;
	let z = 0f64;
	if (iu == 0) { // |f| < 2**-20
		if (f == 0f64) {
			if (k == 0) {
				return 0f64;
			};
			c += (k: f64) * LN2LO;
			return (k: f64) * LN2HI + c;
		};
		R = hfsq * (1.0 - 0.66666666666666666 * f); // Avoid division
		if (k == 0) {
			return f - R;
		};
		return (k: f64) * LN2HI - ((R - ((k: f64) * LN2LO + c)) - f);
	};
	s = f / (2f64 + f);
	z = s * s;
	R = z * (LP1 +
		z * (LP2 +
		z * (LP3 +
		z * (LP4 +
		z * (LP5 +
		z * (LP6 +
		z * LP7))))));
	if (k == 0) {
		return f - (hfsq - s * (hfsq + R));
	};
	return (k: f64) * LN2HI -
		((hfsq - (s * (hfsq + R) + ((k: f64) * LN2LO + c))) - f);
};

// exp(x)
// Returns the exponential of x.
//
// Method
//   1. Argument reduction:
//      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
//      Given x, find r and integer k such that
//
//               x = k*ln2 + r,  |r| <= 0.5*ln2.
//
//      Here r will be represented as r = hi-lo for better
//      accuracy.
//
//   2. Approximation of exp(r) by a special rational function on
//      the interval [0,0.34658]:
//      Write
//          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
//      We use a special Remez algorithm on [0,0.34658] to generate
//      a polynomial of degree 5 to approximate R. The maximum error
//      of this polynomial approximation is bounded by 2**-59. In
//      other words,
//          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
//      (where z=r*r, and the values of P1 to P5 are listed below)
//      and
//          |                  5          |     -59
//          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
//          |                             |
//      The computation of exp(r) thus becomes
//                             2*r
//              exp(r) = 1 + -------
//                            R - r
//                                 r*R1(r)
//                     = 1 + r + ----------- (for better accuracy)
//                                2 - R1(r)
//      where
//                               2       4             10
//              R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
//
//   3. Scale back to obtain exp(x):
//      From step 1, we have
//         exp(x) = 2**k * exp(r)
//
// Special cases:
//      exp(INF) is INF, exp(NaN) is NaN;
//      exp(-INF) is 0, and
//      for finite argument, only exp(0)=1 is exact.
//
// Accuracy:
//      according to an error analysis, the error is always less than
//      1 ulp (unit in the last place).
//
// Misc. info.
//      For IEEE double
//          if x >  7.09782712893383973096e+02 then exp(x) overflow
//          if x < -7.45133219101941108420e+02 then exp(x) underflow
//
// Constants:
// The hexadecimal values are the intended ones for the following
// constants. The decimal values may be used, provided that the
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.

// Returns e^r * 2^k where r = hi - lo and |r| <= (ln(2) / 2).
export fn expmultif64(hi: f64, lo: f64, k: i64) f64 = {
	const P1 = 1.66666666666666657415e-01; // 0x3fc55555; 0x55555555
	const P2 = -2.77777777770155933842e-03; // 0xbf66c16c; 0X16bebd9n
	const P3 = 6.61375632143793436117e-05; // 0x3f11566a; 0Xaf25de2c
	const P4 = -1.65339022054652515390e-06; // 0xbebbbd41; 0Xc5d26bf1
	const P5 = 4.13813679705723846039e-08; // 0x3e663769; 0X72bea4d0

	let r = hi - lo;
	let t = r * r;
	let c = r - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
	const y = 1f64 - ((lo - (r * c) / (2f64 - c)) - hi);
	return ldexpf64(y, k);
};

// Returns e^x.
export fn expf64(x: f64) f64 = {
	const overflow = 7.09782712893383973096e+02;
	const underflow = -7.45133219101941108420e+02;
	const near_zero = 1f64 / ((1i64 << 28i64): f64);

	if (isnan(x) || x == INF) {
		return x;
	} else if (x == -INF) {
		return 0f64;
	} else if (x > overflow) {
		return INF;
	} else if (x < underflow) {
		return 0f64;
	} else if (-near_zero < x && x < near_zero) {
		return 1f64 + x;
	};

	// Reduce; computed as r = hi - lo for extra precision.
	let k = 0i64;
	if (x < 0f64) {
		k = (((LOG2_E * x) - 0.5): i64);
	} else if (x > 0f64) {
		k = (((LOG2_E * x) + 0.5): i64);
	};
	const hi = x - ((k: f64) * LN2_HI);
	const lo = (k: f64) * LN2_LO;

	// Compute
	return expmultif64(hi, lo, k);
};

// Returns 2^x.
export fn exp2f64(x: f64) f64 = {
	const overflow = 1.0239999999999999e+03;
	const underflow = -1.0740e+03;

	if (isnan(x) || x == INF) {
		return x;
	} else if (x == -INF) {
		return 0f64;
	} else if (x > overflow) {
		return INF;
	} else if (x < underflow) {
		return 0f64;
	};

	// Argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
	// Computed as r = hi - lo for extra precision.
	let k = 0i64;
	if (x > 0f64) {
		k = ((x + 0.5): i64);
	} else if (x < 0f64) {
		k = ((x - 0.5): i64);
	};
	const t = x - (k: f64);
	const hi = t * LN2_HI;
	const lo = -t * LN2_LO;

	// Compute
	return expmultif64(hi, lo, k);
};

// __ieee754_sqrt(x)
// Return correctly rounded sqrt.
//           -----------------------------------------
//           | Use the hardware sqrt if you have one |
//           -----------------------------------------
// Method:
//   Bit by bit method using integer arithmetic. (Slow, but portable)
//   1. Normalization
//      Scale x to y in [1,4) with even powers of 2:
//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
//              sqrt(x) = 2**k * sqrt(y)
//   2. Bit by bit computation
//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
//           i                                                   0
//                                     i+1         2
//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
//           i      i            i                 i
//
//      To compute q    from q , one checks whether
//                  i+1       i
//
//                            -(i+1) 2
//                      (q + 2      )  <= y.                     (2)
//                        i
//                                                            -(i+1)
//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
//                             i+1   i             i+1   i
//
//      With some algebraic manipulation, it is not difficult to see
//      that (2) is equivalent to
//                             -(i+1)
//                      s  +  2       <= y                       (3)
//                       i                i
//
//      The advantage of (3) is that s  and y  can be computed by
//                                    i      i
//      the following recurrence formula:
//          if (3) is false
//
//          s     =  s  ,       y    = y   ;                     (4)
//           i+1      i          i+1    i
//
//      otherwise,
//                         -i                      -(i+1)
//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
//           i+1      i          i+1    i     i
//
//      One may easily use induction to prove (4) and (5).
//      Note. Since the left hand side of (3) contain only i+2 bits,
//            it is not necessary to do a full (53-bit) comparison
//            in (3).
//   3. Final rounding
//      After generating the 53 bits result, we compute one more bit.
//      Together with the remainder, we can decide whether the
//      result is exact, bigger than 1/2ulp, or less than 1/2ulp
//      (it will never equal to 1/2ulp).
//      The rounding mode can be detected by checking whether
//      huge + tiny is equal to huge, and whether huge - tiny is
//      equal to huge for some floating point number "huge" and "tiny".

// Returns the square root of x.
export fn sqrtf64(x: f64) f64 = {
	if (x == 0f64) {
		return x;
	} else if (isnan(x) || x == INF) {
		return x;
	} else if (x < 0f64) {
		return NAN;
	};

	let bits = f64bits(x);

	// Normalize x
	let exp = (((bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK): i64);
	if (exp == 0i64) {
		// Subnormal x
		for (bits & (1 << F64_MANTISSA_BITS) == 0) {
			bits <<= 1;
			exp -= 1;
		};
		exp += 1;
	};
	// Unbias exponent
	exp -= (F64_EXPONENT_BIAS: i64);
	bits = bits & ~(F64_EXPONENT_MASK << F64_MANTISSA_BITS);
	bits = bits | (1u64 << (F64_MANTISSA_BITS: u64));
	// Odd exp, double x to make it even
	if (exp & 1i64 == 1i64) {
		bits <<= 1;
	};
	// exp = exp/2, exponent of square root
	exp >>= 1;
	// Generate sqrt(x) bit by bit
	bits <<= 1;
	// q = sqrt(x)
	let q = 0u64;
	let s = 0u64;
	// r = moving bit from MSB to LSB
	let r = ((1u64 << (F64_MANTISSA_BITS + 1u64)): u64);
	for (r != 0) {
		const t = s + r;
		if (t <= bits) {
			s = t + r;
			bits -= t;
			q += r;
		};
		bits <<= 1u64;
		r >>= 1u64;
	};
	// Final rounding
	if (bits != 0) {
		// Remainder, result not exact
		// Round according to extra bit
		q += q & 1;
	};
	// significand + biased exponent
	bits = (q >> 1) + (
		((exp - 1i64 + (F64_EXPONENT_BIAS: i64)): u64) <<
		F64_MANTISSA_BITS);
	return f64frombits(bits);
};

fn is_f64_odd_int(x: f64) bool = {
	const (x_int, x_frac) = modfracf64(x);
	const has_no_frac = (x_frac == 0f64);
	const is_odd = ((x_int: i64 & 1i64) == 1i64);
	return has_no_frac && is_odd;
};

// Returns x^p.
export fn powf64(x: f64, p: f64) f64 = {
	if (x == 1f64 || p == 0f64) {
		return 1f64;
	} else if (p == 1f64) {
		return x;
	} else if (isnan(x)) {
		return NAN;
	} else if (isnan(p)) {
		return NAN;
	} else if (x == 0f64) {
		if (p < 0f64) {
			if (is_f64_odd_int(p)) {
				return copysignf64(INF, x);
			} else {
				return INF;
			};
		} else if (p > 0f64) {
			if (is_f64_odd_int(p)) {
				return x;
			} else {
				return 0f64;
			};
		};
	} else if (isinf(p)) {
		if (x == -1f64) {
			return 1f64;
		} else if ((absf64(x) < 1f64) == (p == INF)) {
			return 0f64;
		};
		return INF;
	} else if (isinf(x)) {
		if (x == -INF) {
			return powf64(-0f64, -p);
		} else if (p < 0f64) {
			return 0f64;
		} else if (p > 0f64) {
			return INF;
		};
	} else if (p == 0.5f64) {
		return sqrtf64(x);
	} else if (p == -0.5f64) {
		return 1f64 / sqrtf64(x);
	};

	let (p_int, p_frac) = modfracf64(absf64(p));
	if (p_frac != 0f64 && x < 0f64) {
		return NAN;
	};
	if (p_int > types::I64_MAX: f64) {
		if (x == -1f64) {
			return 1f64;
		} else if ((absf64(x) < 1f64) == (p > 0f64)) {
			return 0f64;
		} else {
			return INF;
		};
	};

	let res_mantissa = 1f64;
	let res_exp = 0i64;

	// The method used later in this function doesn't apply to fractional
	// powers, so we have to handle these separately with
	// x^p = e^{p * ln(x)}
	if (p_frac != 0f64) {
		if (p_frac > 0.5f64) {
			p_frac -= 1f64;
			p_int += 1f64;
		};
		res_mantissa = expf64(p_frac * logf64(x));
	};

	// Repeatedly square our number x, for each bit in our power p.
	// If the current bit is 1 in p, add the respective power of x to our
	// result.
	let (x_mantissa, x_exp) = frexpf64(x);
	for (let i = p_int: i64; i != 0; i >>= 1) {
		// Check for over/underflow.
		if (x_exp <= -1i64 << (F64_EXPONENT_BITS: i64)) {
			return 0f64;
		};
		if (x_exp >= 1i64 << (F64_EXPONENT_BITS: i64)) {
			return INF;
		};
		// Perform squaring.
		if (i & 1i64 == 1i64) {
			res_mantissa *= x_mantissa;
			res_exp += x_exp;
		};
		x_mantissa *= x_mantissa;
		x_exp <<= 1;
		// Correct mantisa to be in [0.5, 1).
		if (x_mantissa < 0.5f64) {
			x_mantissa += x_mantissa;
			x_exp -= 1;
		};
	};

	if (p < 0f64) {
		res_mantissa = 1f64 / res_mantissa;
		res_exp = -res_exp;
	};

	return ldexpf64(res_mantissa, res_exp);
};

// Returns the greatest integer value less than or equal to x.
export fn floorf64(x: f64) f64 = {
	if (x == 0f64 || isnan(x) || isinf(x)) {
		return x;
	};
	if (x < 0f64) {
		let (int_part, frac_part) = modfracf64(-x);
		if (frac_part != 0f64) {
			int_part += 1f64;
		};
		return -int_part;
	};
	return modfracf64(x).0;
};

// Returns the least integer value greater than or equal to x.
export fn ceilf64(x: f64) f64 = -floorf64(-x);

// Returns the integer value of x.
export fn truncf64(x: f64) f64 = {
	if (x == 0f64 || isnan(x) || isinf(x)) {
		return x;
	};
	return modfracf64(x).0;
};

// Returns the nearest integer, rounding half away from zero.
export fn roundf64(x: f64) f64 = {
	let bits = f64bits(x);
	let e = (bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK;
	if (e < F64_EXPONENT_BIAS) {
		// Round abs(x) < 1 including denormals.
		bits &= F64_SIGN_MASK; // +-0
		if (e == F64_EXPONENT_BIAS - 1) {
			bits |= F64_ONE; // +-1
		};
	} else if (e < F64_EXPONENT_BIAS + F64_MANTISSA_BITS) {
		// Round any abs(x) >= 1 containing a fractional component
		// [0,1).
		// Numbers with larger exponents are returned unchanged since
		// they must be either an integer, infinity, or NaN.
		const half = 1u64 << (F64_MANTISSA_BITS - 1);
		e -= F64_EXPONENT_BIAS;
		bits += half >> e;
		bits = bits & ~(F64_MANTISSA_MASK >> e);
	};
	return f64frombits(bits);
};

// Returns the floating-point remainder of x / y. The magnitude of the result
// is less than y and its sign agrees with that of x.
export fn modf64(x: f64, y: f64) f64 = {
	if (y == 0f64) {
		return NAN;
	};
	if (isinf(x) || isnan(x) || isnan(y)) {
		return NAN;
	};

	y = absf64(y);

	const (y_frac, y_exp) = frexpf64(y);
	let r = x;
	if (x < 0f64) {
		r = -x;
	};

	for (r >= y) {
		const (r_frac, r_exp) = frexpf64(r);
		if (r_frac < y_frac) {
			r_exp -= 1i64;
		};
		r = r - ldexpf64(y, r_exp - y_exp);
	};
	if (x < 0f64) {
		r = -r;
	};
	return r;
};