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

// Sections of the code below are based on Go's implementation, which is, in
// turn, based on:
// * the Cephes Mathematical Library (cephes/cmath/{sin,..}), available from
//   http://www.netlib.org/cephes/cmath.tgz.
// * FreeBSD's /usr/src/lib/msun/src/{s_asinh.c,...}
// The original C code, as well as the respective comments and constants are
// from these libraries.
//
// The Cephes copyright notice:
// ====================================================
// Cephes Math Library Release 2.8:  June, 2000
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
//
// The readme file at http://netlib.sandia.gov/cephes/ says:
//    Some software in this archive may be from the book _Methods and
// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
// International, 1989) or from the Cephes Mathematical Library, a
// commercial product. In either event, it is copyrighted by the author.
// What you see here may be used freely but it comes with no support or
// guarantee.
//
//   The two known misprints in the book are repaired here in the
// source listings for the gamma function and the incomplete beta
// integral.
//
//   Stephen L. Moshier
//   moshier@na-net.ornl.gov
// ====================================================
//
// 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.
// ====================================================

// sin coefficients
const SIN_CF: [_]f64 = [
	1.58962301576546568060e-10, // 0x3de5d8fd1fd19ccd
	-2.50507477628578072866e-8, // 0xbe5ae5e5a9291f5d
	2.75573136213857245213e-6, // 0x3ec71de3567d48a1
	-1.98412698295895385996e-4, // 0xbf2a01a019bfdf03
	8.33333333332211858878e-3, // 0x3f8111111110f7d0
	-1.66666666666666307295e-1, // 0xbfc5555555555548
];

// cos coefficients
const COS_CF: [_]f64 = [
	-1.13585365213876817300e-11, // 0xbda8fa49a0861a9b
	2.08757008419747316778e-9, // 0x3e21ee9d7b4e3f05
	-2.75573141792967388112e-7, // 0xbe927e4f7eac4bc6
	2.48015872888517045348e-5, // 0x3efa01a019c844f5
	-1.38888888888730564116e-3, // 0xbf56c16c16c14f91
	4.16666666666665929218e-2, // 0x3fa555555555554b
];

// PI / 4 split into three parts
def PI4A: f64 = 7.85398125648498535156e-1; // 0x3fe921fb40000000
def PI4B: f64 = 3.77489470793079817668e-8; // 0x3e64442d00000000
def PI4C: f64 = 2.69515142907905952645e-15; // 0x3ce8469898cc5170

// reduce_threshold is the maximum value of x where the reduction using PI/4
// in 3 float64 parts still gives accurate results. This threshold
// is set by y*C being representable as a float64 without error
// where y is given by y = floor(x * (4 / PI)) and C is the leading partial
// terms of 4/PI. Since the leading terms (PI4A and PI4B in sin.go) have 30
// and 32 trailing zero bits, y should have less than 30 significant bits.
//	y < 1<<30  -> floor(x*4/PI) < 1<<30 -> x < (1<<30 - 1) * PI/4
// So, conservatively we can take x < 1<<29.
// Above this threshold Payne-Hanek range reduction must be used.
def REDUCE_THRESHOLD: f64 = ((1u64 << 29): f64);

// MPI4 is the binary digits of 4/pi as a uint64 array,
// that is, 4/pi = Sum MPI4[i]*2^(-64*i)
// 19 64-bit digits and the leading one bit give 1217 bits
// of precision to handle the largest possible float64 exponent.
const MPI4: [_]u64 = [
	0x0000000000000001,
	0x45f306dc9c882a53,
	0xf84eafa3ea69bb81,
	0xb6c52b3278872083,
	0xfca2c757bd778ac3,
	0x6e48dc74849ba5c0,
	0x0c925dd413a32439,
	0xfc3bd63962534e7d,
	0xd1046bea5d768909,
	0xd338e04d68befc82,
	0x7323ac7306a673e9,
	0x3908bf177bf25076,
	0x3ff12fffbc0b301f,
	0xde5e2316b414da3e,
	0xda6cfd9e4f96136e,
	0x9e8c7ecd3cbfd45a,
	0xea4f758fd7cbe2f6,
	0x7a0e73ef14a525d4,
	0xd7f6bf623f1aba10,
	0xac06608df8f6d757,
];

// trig_reduce implements Payne-Hanek range reduction by PI/4 for x > 0. It
// returns the integer part mod 8 (j) and the fractional part (z) of x / (PI/4).
fn trig_reduce(x: f64) (u64, f64) = {
	// The implementation is based on: "ARGUMENT REDUCTION FOR HUGE
	// ARGUMENTS: Good to the Last Bit" K. C. Ng et al, March 24, 1992
	// The simulated multi-precision calculation of x * B uses 64-bit
	// integer arithmetic.
	const PI4 = PI / 4f64;
	if (x < PI4) {
		return (0u64, x);
	};
	// Extract out the integer and exponent such that x = ix * 2 ** exp
	let ix = f64bits(x);
	const exp =
		((ix >> F64_MANTISSA_BITS &
			F64_EXPONENT_MASK): i64) -
		(F64_EXPONENT_BIAS: i64) -
		(F64_MANTISSA_BITS: i64);
	ix = ix & ~(F64_EXPONENT_MASK << F64_MANTISSA_BITS);
	ix |= 1 << F64_MANTISSA_BITS;
	// Use the exponent to extract the 3 appropriate uint64 digits from
	// MPI4, B ~ (z0, z1, z2), such that the product leading digit has the
	// exponent -61. Note, exp >= -53 since x >= PI4 and exp < 971 for
	// maximum float64.
	const digit = ((exp + 61): u64) / 64;
	const bitshift = ((exp + 61): u64) % 64;
	const z0 = (MPI4[digit] << bitshift) |
		(MPI4[digit + 1] >> (64 - bitshift));
	const z1 = (MPI4[digit + 1] << bitshift) |
		(MPI4[digit + 2] >> (64 - bitshift));
	const z2 = (MPI4[digit + 2] << bitshift) |
		(MPI4[digit + 3] >> (64 - bitshift));
	// Multiply mantissa by the digits and extract the upper two digits
	// (hi, lo).
	const (z2hi, _) = mulu64(z2, ix);
	const (z1hi, z1lo) = mulu64(z1, ix);
	const z0lo = z0 * ix;
	const lo = z1lo + z2hi;
	let hi = z0lo + z1hi;
	hi += (z1lo >> 63) & (z2hi >> 63); // carry from lo
	// The top 3 bits are j.
	let j = hi >> 61;
	// Extract the fraction and find its magnitude.
	hi = hi << 3 | lo >> 61;
	const lz = ((leading_zeros_u64(hi)): uint);
	const e = ((F64_EXPONENT_BIAS - (lz + 1)): u64);
	// Clear implicit mantissa bit and shift into place.
	hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1)));
	hi >>= 64 - F64_MANTISSA_BITS;
	// Include the exponent and convert to a float.
	hi |= e << F64_MANTISSA_BITS;
	let z = f64frombits(hi);
	// Map zeros to origin.
	if (j & 1 == 1) {
		j += 1;
		j &= 7;
		z -= 1f64;
	};
	// Multiply the fractional part by pi/4.
	return (j, z * PI4);
};

//      cos.c
//
//      Circular cosine
//
// SYNOPSIS:
//
// double x, y, cos();
// y = cos( x );
//
// DESCRIPTION:
//
// Range reduction is into intervals of pi/4.  The reduction error is nearly
// eliminated by contriving an extended precision modular arithmetic.
//
// Two polynomial approximating functions are employed.
// Between 0 and pi/4 the cosine is approximated by
//      1  -  x**2 Q(x**2).
// Between pi/4 and pi/2 the sine is represented as
//      x  +  x**3 P(x**2).
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain      # trials      peak         rms
//    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
//    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18

// Returns the cosine of x, in radians.
export fn cosf64(x: f64) f64 = {
	if (isnan(x) || isinf(x)) {
		return NAN;
	};

	// Make argument positive
	let is_negative = false;
	x = absf64(x);

	let j = 0u64;
	let y = 0f64;
	let z = 0f64;
	if (x >= REDUCE_THRESHOLD) {
		const reduce_res = trig_reduce(x);
		j = reduce_res.0;
		z = reduce_res.1;
	} else {
		// Integer part of x/(PI/4), as integer for tests on the phase
		// angle
		j = ((x * (4f64 / PI)): i64: u64);
		// Integer part of x/(PI/4), as float
		y = (j: i64: f64);

		// Map zeros to origin
		if (j & 1 == 1) {
			j += 1;
			y += 1f64;
		};
		// Octant modulo 2PI radians (360 degrees)
		j &= 7;
		// Extended precision modular arithmetic
		z = ((x - (y * PI4A)) - (y * PI4B)) - (y * PI4C);
	};

	if (j > 3) {
		j -= 4;
		is_negative = !is_negative;
	};
	if (j > 1) {
		is_negative = !is_negative;
	};

	const zz = z * z;
	if (j == 1 || j == 2) {
		y = z + z * zz * ((((((SIN_CF[0] * zz) +
			SIN_CF[1]) * zz +
			SIN_CF[2]) * zz +
			SIN_CF[3]) * zz +
			SIN_CF[4]) * zz +
			SIN_CF[5]);
	} else {
		y = 1.0 - 0.5 * zz + zz * zz * ((((((COS_CF[0] * zz) +
			COS_CF[1]) * zz +
			COS_CF[2]) * zz +
			COS_CF[3]) * zz +
			COS_CF[4]) * zz +
			COS_CF[5]);
	};
	if (is_negative) {
		y = -y;
	};
	return y;
};

//      sin.c
//
//      Circular sine
//
// SYNOPSIS:
//
// double x, y, sin();
// y = sin( x );
//
// DESCRIPTION:
//
// Range reduction is into intervals of pi/4.  The reduction error is nearly
// eliminated by contriving an extended precision modular arithmetic.
//
// Two polynomial approximating functions are employed.
// Between 0 and pi/4 the sine is approximated by
//      x  +  x**3 P(x**2).
// Between pi/4 and pi/2 the cosine is represented as
//      1  -  x**2 Q(x**2).
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain      # trials      peak         rms
//    DEC       0, 10       150000       3.0e-17     7.8e-18
//    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
//
// Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9.  The loss
// is not gradual, but jumps suddenly to about 1 part in 10e7.  Results may
// be meaningless for x > 2**49 = 5.6e14.

// Returns the sine of x, in radians.
export fn sinf64(x: f64) f64 = {
	if (x == 0f64 || isnan(x)) {
		return x;
	} else if (isinf(x)) {
		return NAN;
	};

	// Make argument positive but save the sign
	let is_negative = false;
	if (x < 0f64) {
		x = -x;
		is_negative = true;
	};

	let j = 0u64;
	let y = 0f64;
	let z = 0f64;
	if (x >= REDUCE_THRESHOLD) {
		const reduce_res = trig_reduce(x);
		j = reduce_res.0;
		z = reduce_res.1;
	} else {
		// Integer part of x/(PI/4), as integer for tests on the phase
		// angle
		j = ((x * (4f64 / PI)): i64: u64);
		// Integer part of x/(PI/4), as float
		y = (j: i64: f64);

		// Map zeros to origin
		if (j & 1 == 1) {
			j += 1;
			y += 1f64;
		};
		// Octant modulo 2PI radians (360 degrees)
		j &= 7;
		// Extended precision modular arithmetic
		z = ((x - (y * PI4A)) - (y * PI4B)) - (y * PI4C);
	};

	// Reflect in x axis
	if (j > 3) {
		j -= 4;
		is_negative = !is_negative;
	};

	const zz = z * z;
	if (j == 1 || j == 2) {
		y = 1.0 - 0.5 * zz + zz * zz *
			((((((COS_CF[0] * zz) +
				COS_CF[1]) * zz +
				COS_CF[2]) * zz +
				COS_CF[3]) * zz +
				COS_CF[4]) * zz +
				COS_CF[5]);
	} else {
		y = z + z * zz *
			((((((SIN_CF[0] * zz) +
				SIN_CF[1]) * zz +
				SIN_CF[2]) * zz +
				SIN_CF[3]) * zz +
				SIN_CF[4]) * zz +
				SIN_CF[5]);
	};
	if (is_negative) {
		y = -y;
	};
	return y;
};

//      tan.c
//
//      Circular tangent
//
// SYNOPSIS:
//
// double x, y, tan();
// y = tan( x );
//
// DESCRIPTION:
//
// Returns the circular tangent of the radian argument x.
//
// Range reduction is modulo pi/4.  A rational function
//       x + x**3 P(x**2)/Q(x**2)
// is employed in the basic interval [0, pi/4].
//
// ACCURACY:
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    DEC      +-1.07e9      44000      4.1e-17     1.0e-17
//    IEEE     +-1.07e9      30000      2.9e-16     8.1e-17
//
// Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9.  The loss
// is not gradual, but jumps suddenly to about 1 part in 10e7.  Results may
// be meaningless for x > 2**49 = 5.6e14.
// [Accuracy loss statement from sin.go comments.]

// tan coefficients
const TAN_P: [_]f64 = [
	-1.30936939181383777646e4, // 0xc0c992d8d24f3f38
	1.15351664838587416140e6, // 0x413199eca5fc9ddd
	-1.79565251976484877988e7, // 0xc1711fead3299176
];
const TAN_Q: [_]f64 = [
	1.00000000000000000000e0,
	1.36812963470692954678e4, // 0x40cab8a5eeb36572
	-1.32089234440210967447e6, // 0xc13427bc582abc96
	2.50083801823357915839e7, // 0x4177d98fc2ead8ef
	-5.38695755929454629881e7, // 0xc189afe03cbe5a31
];

// Returns the tangent of x, in radians.
export fn tanf64(x: f64) f64 = {
	if (x == 0f64 || isnan(x)) {
		return x;
	} else if (isinf(x)) {
		return NAN;
	};

	// Make argument positive but save the sign
	let is_negative = false;
	if (x < 0f64) {
		x = -x;
		is_negative = true;
	};
	let j = 0u64;
	let y = 0f64;
	let z = 0f64;
	if (x >= REDUCE_THRESHOLD) {
		const reduce_res = trig_reduce(x);
		j = reduce_res.0;
		z = reduce_res.1;
	} else {
		// Integer part of x/(PI/4), as integer for tests on the phase
		// angle
		j = ((x * (4f64 / PI)): i64: u64);
		// Integer part of x/(PI/4), as float
		y = (j: i64: f64);

		// Map zeros and singularities to origin
		if (j & 1 == 1) {
			j += 1;
			y += 1f64;
		};

		z = ((x - (y * PI4A)) - (y * PI4B)) - (y * PI4C);
	};
	const zz = z * z;

	if (zz > 1e-14) {
		y = z + z * (zz *
			(((TAN_P[0] * zz) + TAN_P[1]) * zz + TAN_P[2]) /
			((((zz + TAN_Q[1]) * zz +
				TAN_Q[2]) * zz +
				TAN_Q[3]) * zz +
				TAN_Q[4]));
	} else {
		y = z;
	};
	if (j & 2 == 2) {
		y = -1f64 / y;
	};
	if (is_negative) {
		y = -y;
	};
	return y;
};

// Evaluates a series valid in the range [0, 0.66].
fn xatan(x: f64) f64 = {
	const P0 = -8.750608600031904122785e-01;
	const P1 = -1.615753718733365076637e+01;
	const P2 = -7.500855792314704667340e+01;
	const P3 = -1.228866684490136173410e+02;
	const P4 = -6.485021904942025371773e+01;
	const Q0 =  2.485846490142306297962e+01;
	const Q1 =  1.650270098316988542046e+02;
	const Q2 =  4.328810604912902668951e+02;
	const Q3 =  4.853903996359136964868e+02;
	const Q4 =  1.945506571482613964425e+02;
	let z = x * x;
	z = z * ((((P0 * z + P1) * z + P2) * z + P3) * z + P4) /
		(((((z + Q0) * z + Q1) * z + Q2) * z + Q3) * z + Q4);
	z = (x * z) + x;
	return z;
};

// Reduces argument (known to be positive) to the range [0, 0.66] and calls
// xatan.
fn satan(x: f64) f64 = {
	// pi / 2 = PIO2 + morebits
	const morebits = 6.123233995736765886130e-17;
	// tan(3 * pi / 8)
	const tan3pio8 = 2.41421356237309504880;
	if (x <= 0.66) {
		return xatan(x);
	};
	if (x > tan3pio8) {
		return (PI / 2f64) - xatan(1f64 / x) + morebits;
	};
	return (PI / 4f64) +
		xatan((x - 1f64) / (x + 1f64)) +
		(0.5f64 * morebits);
};

// Returns the arcsine, in radians, of x.
export fn asinf64(x: f64) f64 = {
	if (x == 0f64) {
		return x;
	};
	let is_negative = false;
	if (x < 0.064) {
		x = -x;
		is_negative = true;
	};
	if (x > 1f64) {
		return NAN;
	};
	let temp = sqrtf64(1f64 - x * x);
	if (x > 0.7f64) {
		temp = PI / 2f64 - satan(temp / x);
	} else {
		temp = satan(x / temp);
	};

	if (is_negative) {
		temp = -temp;
	};
	return temp;
};

// Returns the arccosine, in radians, of x.
export fn acosf64(x: f64) f64 = {
	return PI / 2f64 - asinf64(x);
};

// atan.c
// Inverse circular tangent (arctangent)
//
// SYNOPSIS:
// double x, y, atan();
// y = atan( x );
//
// DESCRIPTION:
// Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
//
// Range reduction is from three intervals into the interval from zero to 0.66.
// The approximant uses a rational function of degree 4/5 of the form
// x + x**3 P(x)/Q(x).
//
// ACCURACY:
//                      Relative error:
// arithmetic   domain    # trials  peak     rms
//    DEC       -10, 10   50000     2.4e-17  8.3e-18
//    IEEE      -10, 10   10^6      1.8e-16  5.0e-17

// Returns the arctangent, in radians, of x.
export fn atanf64(x: f64) f64 = {
	if (x == 0f64) {
		return x;
	};
	if (x > 0f64) {
		return satan(x);
	};
	return -satan(-x);
};

// Floating-point hyperbolic sine and cosine.
// The exponential func is called for arguments greater in magnitude than 0.5.
// A series is used for arguments smaller in magnitude than 0.5.
// Cosh(x) is computed from the exponential func for all arguments.

// Returns the hyperbolic sine of x.
export fn sinhf64(x: f64) f64 = {
	// The coefficients are #2029 from Hart & Cheney. (20.36D)
	const P0 = -0.6307673640497716991184787251e+6;
	const P1 = -0.8991272022039509355398013511e+5;
	const P2 = -0.2894211355989563807284660366e+4;
	const P3 = -0.2630563213397497062819489e+2;
	const Q0 = -0.6307673640497716991212077277e+6;
	const Q1 = 0.1521517378790019070696485176e+5;
	const Q2 = -0.173678953558233699533450911e+3;

	let is_negative = false;
	if (x < 0f64) {
		x = -x;
		is_negative = true;
	};

	let temp = 0f64;
	if (x > 21f64) {
		temp = expf64(x) * 0.5f64;
	} else if (x > 0.5f64) {
		const ex = expf64(x);
		temp = (ex - (1f64 / ex)) * 0.5f64;
	} else {
		const sq = x * x;
		temp = (((P3 * sq + P2) * sq + P1) * sq + P0) * x;
		temp = temp / (((sq + Q2) * sq + Q1) * sq + Q0);
	};

	if (is_negative) {
		temp = -temp;
	};

	return temp;
};

// Returns the hyperbolic cosine of x.
export fn coshf64(x: f64) f64 = {
	x = absf64(x);
	if (x > 21f64) {
		return expf64(x) * 0.5f64;
	};
	const ex = expf64(x);
	return (ex + 1f64 / ex) * 0.5f64;
};

//      tanh.c
//
//      Hyperbolic tangent
//
// SYNOPSIS:
//
// double x, y, tanh();
//
// y = tanh( x );
//
// DESCRIPTION:
//
// Returns hyperbolic tangent of argument in the range MINLOG to MAXLOG.
//      MAXLOG = 8.8029691931113054295988e+01 = log(2**127)
//      MINLOG = -8.872283911167299960540e+01 = log(2**-128)
//
// A rational function is used for |x| < 0.625.  The form
// x + x**3 P(x)/Q(x) of Cody & Waite is employed.
// Otherwise,
//      tanh(x) = sinh(x)/cosh(x) = 1  -  2/(exp(2x) + 1).
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    IEEE      -2,2        30000       2.5e-16     5.8e-17

// tanh coefficients
const TANH_P: [_]f64 = [
	-9.64399179425052238628e-1,
	-9.92877231001918586564e1,
	-1.61468768441708447952e3,
];
const TANH_Q: [_]f64 = [
	1.12811678491632931402e2,
	2.23548839060100448583e3,
	4.84406305325125486048e3,
];

// Returns the hyperbolic tangent of x.
export fn tanhf64(x: f64) f64 = {
	const MAXLOG = 8.8029691931113054295988e+01; // log(2**127)
	let z = absf64(x);
	if (z > 0.5f64 * MAXLOG) {
		if (x < 0f64) {
			return -1f64;
		};
		return 1f64;
	} else if (z >= 0.625f64) {
		const s = expf64(2f64 * z);
		z = 1f64 - 2f64 / (s + 1f64);
		if (x < 0f64) {
			z = -z;
		};
	} else {
		if (x == 0f64) {
			return x;
		};
		const s = x * x;
		z = x + x * s * ((TANH_P[0] * s + TANH_P[1]) * s + TANH_P[2]) /
			(((s + TANH_Q[0]) * s + TANH_Q[1]) * s + TANH_Q[2]);
	};
	return z;
};

// asinh(x)
// Method :
//	Based on
//	        asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
//	we have
//	asinh(x) := x  if  1+x*x=1,
//	         := sign(x)*(log(x)+ln2)) for large |x|, else
//	         := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
//	         := sign(x)*log1p(|x| + x**2/(1 + sqrt(1+x**2)))
//

// Returns the inverse hyperbolic sine of x.
export fn asinhf64(x: f64) f64 = {
	const NEAR_ZERO = 1f64 / ((1i64 << 28): f64);
	const LARGE = ((1i64 << 28): f64);

	if (isnan(x) || isinf(x)) {
		return x;
	};

	let is_negative = false;
	if (x < 0f64) {
		x = -x;
		is_negative = true;
	};

	let temp = 0f64;

	if (x > LARGE) {
		// |x| > 2**28
		temp = logf64(x) + LN_2;
	} else if (x > 2f64) {
		// 2**28 > |x| > 2.0
		temp = logf64(2f64 * x +
			1f64 / (sqrtf64(x * x + 1f64) + x));
	} else if (x < NEAR_ZERO) {
		// |x| < 2**-28
		temp = x;
	} else {
		// 2.0 > |x| > 2**-28
		temp = log1pf64(x + x * x /
			(1f64 + sqrtf64(1f64 + x * x)));
	};
	if (is_negative) {
		temp = -temp;
	};
	return temp;
};

// __ieee754_acosh(x)
// Method :
//	Based on
//	        acosh(x) = log [ x + sqrt(x*x-1) ]
//	we have
//	        acosh(x) := log(x)+ln2,	if x is large; else
//	        acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
//	        acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
//
// Special cases:
//	acosh(x) is NaN with signal if x<1.
//	acosh(NaN) is NaN without signal.
//

// Returns the inverse hyperbolic cosine of x.
export fn acoshf64(x: f64) f64 = {
	const LARGE = ((1i64 << 28): f64);

	if (x < 1f64 || isnan(x)) {
		return NAN;
	} else if (x == 1f64) {
		return 0f64;
	} else if (x >= LARGE) {
		// x > 2**28
		return logf64(x) + LN_2;
	} else if (x > 2f64) {
		// 2**28 > x > 2
		return logf64(2f64 * x - 1f64 /
			(x + sqrtf64(x * x - 1f64)));
	};
	const t = x - 1f64;
	// 2 >= x > 1
	return log1pf64(t + sqrtf64(2f64 * t + t * t));
};

// __ieee754_atanh(x)
// Method :
//	1. Reduce x to positive by atanh(-x) = -atanh(x)
//	2. For x>=0.5
//	            1              2x                          x
//	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
//	            2             1 - x                      1 - x
//
//	For x<0.5
//	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
//
// Special cases:
//	atanh(x) is NaN if |x| > 1 with signal;
//	atanh(NaN) is that NaN with no signal;
//	atanh(+-1) is +-INF with signal.
//

// Returns the inverse hyperbolic tangent of x.
export fn atanhf64(x: f64) f64 = {
	const NEAR_ZERO = 1f64 / ((1i64 << 28): f64);

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

	let is_negative = false;

	if (x < 0f64) {
		x = -x;
		is_negative = true;
	};

	let temp = 0f64;

	if (x < NEAR_ZERO) {
		temp = x;
	} else if (x < 0.5f64) {
		temp = x + x;
		temp = 0.5f64 * log1pf64(temp + temp * x / (1f64 - x));
	} else {
		temp = 0.5f64 * log1pf64((x + x) / (1f64 - x));
	};
	if (is_negative) {
		temp = -temp;
	};
	return temp;
};

// Returns the arctangent, in radians, of y / x.
export fn atan2f64(y: f64, x: f64) f64 = {
	if (isnan(y) || isnan(x)) {
		return NAN;
	} else if (y == 0f64) {
		x = if (x >= 0f64 && signf64(x) > 0) 0f64 else PI;
		return copysignf64(x, y);
	} else if (x == 0f64) {
		return copysignf64(PI / 2f64, y);
	} else if (isinf(x)) {
		if (signf64(x) > 0) {
			x = if (isinf(y)) PI / 4f64 else 0f64;
			return copysignf64(x, y);
		} else {
			x = if (isinf(y)) 3f64 * PI / 4f64 else PI;
			return copysignf64(x, y);
		};
	} else if (isinf(y)) {
		return copysignf64(PI / 2f64, y);
	};

	const q = atanf64(y / x);
	if (x < 0f64) {
		return if (q <= 0f64) q + PI else q - PI;
	};
	return q;
};

// Returns the square root of a*a + b*b, taking care to avoid unnecessary
// overflow and underflow.
export fn hypotf64(a: f64, b: f64) f64 = {
	if (isinf(a) || isinf(b)) {
		return INF;
	} else if (isnan(a) || isnan(b)) {
		return NAN;
	};
	a = absf64(a);
	b = absf64(b);
	if (a < b) {
		const temp = a;
		a = b;
		b = temp;
	};
	if (a == 0f64) {
		return 0f64;
	};
	b = b / a;
	return a * sqrtf64(1f64 + b * b);
};