hare

complex.ha (21042B)

---

 // 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 Cephes Math Library. The original C code can be found at
 // http://netlib.sandia.gov/cephes/c9x-complex/.
 //
 // 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 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 math;
 use math::checked;
 
 // A complex number containing a real component and an imaginary component,
 // represented as two single-precision floating point numbers.
 export type c64 = (f32, f32);
 
 // A complex number containing a real component and an imaginary component,
 // represented as two double-precision floating point numbers.
 export type c128 = (f64, f64);
 
 // A tagged union of all complex types.
 export type complex = (c64 | c128);
 
 // Converts a [[c64]] to a [[c128]].
 export fn c64to128(z: c64) c128 = (z.0: f64, z.1: f64);
 
 // Truncates a [[c128]] to a [[c64]]. Precision may be lost.
 export fn c128to64(z: c128) c64 = (z.0: f32, z.1: f32);
 
 // Adds two complex numbers
 export fn addc64(a: c64, b: c64) c64 = (a.0 + b.0, a.1 + b.1);
 
 // Adds two complex numbers.
 export fn addc128(a: c128, b: c128) c128 = (a.0 + b.0, a.1 + b.1);
 
 // Subtracts two complex numbers.
 export fn subc64(a: c64, b: c64) c64 = (a.0 - b.0, a.1 - b.1);
 
 // Subtracts two complex numbers.
 export fn subc128(a: c128, b: c128) c128 = (a.0 - b.0, a.1 - b.1);
 
 // Multiplies two complex numbers.
 export fn mulc64(a: c64, b: c64) c64 =
 	(a.0 * b.0 - a.1 * b.1, a.1 * b.0 + a.0 * b.1);
 
 // Multiplies two complex numbers.
 export fn mulc128(a: c128, b: c128) c128 =
 	(a.0 * b.0 - a.1 * b.1, a.1 * b.0 + a.0 * b.1);
 
 // Divides two complex numbers.
 export fn divc64(a: c64, b: c64) c64 = {
 	const denom = b.0 * b.0 + b.1 * b.1;
 	return (
 		(a.0 * b.0 + a.1 * b.1) / denom,
 		(a.1 * b.0 - a.0 * b.1) / denom,
 	);
 };
 
 // Divides two complex numbers.
 export fn divc128(a: c128, b: c128) c128 = {
 	const denom = b.0 * b.0 + b.1 * b.1;
 	return (
 		(a.0 * b.0 + a.1 * b.1) / denom,
 		(a.1 * b.0 - a.0 * b.1) / denom,
 	);
 };
 
 // Takes the conjugate of a complex number by negating the imaginary component.
 export fn conjc64(z: c64) c64 = (z.0, -z.1);
 
 // Takes the conjugate of a complex number by negating the imaginary component.
 export fn conjc128(z: c128) c128 = (z.0, -z.1);
 
 // Takes the absolute value of a complex number.
 export fn absc128(z: c128) f64 = math::hypotf64(z.0, z.1);
 
 // Gets the argument, or phase, of a complex number.
 export fn argc128(z: c128) f64 = math::atan2f64(z.1, z.0);
 
 // Checks if two complex numbers are equal. Be sure to take floating point
 // round-off errors into account.
 export fn equalc64(a: c64, b: c64) bool = a.0 == b.0 && a.1 == b.1;
 
 // Checks if two complex numbers are equal. Be sure to take floating point
 // round-off errors into account.
 export fn equalc128(a: c128, b: c128) bool = a.0 == b.0 && a.1 == b.1;
 
 // Checks if two complex numbers are equal. Be sure to take floating point
 // round-off errors into account.
 export fn equal(a: complex, b: complex) bool = {
 	match (a) {
 	case let a: c64 =>
 		return equalc64(a, b as c64);
 	case let a: c128 =>
 		return equalc128(a, b as c128);
 	};
 };
 
 // Returns [[math::E]] raised to the power of z.
 export fn expc128(z: c128) c128 = {
 	if (math::isinf(z.0)) {
 		if (z.0 > 0f64 && z.1 == 0f64) {
 			return z;
 		};
 		if (math::isinf(z.1) || math::isnan(z.1)) {
 			if (z.0 < 0f64) {
 				return (0f64, math::copysignf64(0f64, z.1));
 			} else {
 				return (math::INF, math::NAN);
 			};
 		};
 	} else if (math::isnan(z.0) && z.1 == 0f64) {
 		return (math::NAN, z.1);
 	};
 	return rectc128(math::expf64(z.0), z.1);
 };
 
 // Returns true if the given complex number is infinite.
 export fn isinf(z: c128) bool = math::isinf(z.0) || math::isinf(z.1);
 
 // Returns true if the given complex number is NaN -- that is -- either
 // component is NaN and neither component is an infinity.
 export fn isnan(z: c128) bool =
 	!isinf(z) && (math::isnan(z.0) || math::isnan(z.1));
 
 // Returns the natural logarithm of z.
 export fn logc128(z: c128) c128 = (math::logf64(absc128(z)), argc128(z));
 
 // Negates z.
 export fn negc64(z: c64) c64 = (-z.0, -z.1);
 
 // Negates z.
 export fn negc128(z: c128) c128 = (-z.0, -z.1);
 
 // Creates a new [[c128]] from the polar coordinates (r, theta).
 export fn rectc128(r: f64, theta: f64) c128 =
 	(r * math::cosf64(theta), r * math::sinf64(theta));
 
 // Returns the polar coordinates of z.
 export fn polarc128(z: c128) (f64, f64) = (absc128(z), argc128(z));
 
 // Returns a raised to the power of b.
 export fn powc128(a: c128, b: c128) c128 = {
 	if (a.0 == 0f64 && a.1 == 0f64) {
 		if (isnan(b)) {
 			return (math::NAN, math::NAN);
 		} else if (b.0 == 0f64) {
 			return (1f64, 0f64);
 		} else if (b.0 < 0f64) {
 			return (math::INF, if (b.1 == 0f64) 0f64 else math::INF);
 		} else {
 			assert(b.0 > 0f64);
 			return (0f64, 0f64);
 		};
 	};
 	const mod = absc128(a);
 	if (mod == 0f64) {
 		return (0f64, 0f64);
 	};
 	let r = math::powf64(mod, b.0);
 	const phase = argc128(a);
 	let theta = b.0 * phase;
 	if (b.1 != 0f64) {
 		r *= math::expf64(-b.1 * phase);
 		theta += b.1 * math::logf64(mod);
 	};
 	return rectc128(r, theta);
 };
 
 // Projects z onto the surface of a Riemann Sphere. If z is finite, it projects
 // to itself. If z is infinite, it projects to positive infinity on the real
 // axis.
 export fn projc64(z: c64) c64 =
 	if (!isinf(c64to128(z))) z else (math::INF, math::copysignf32(0f32, z.1));
 
 // Projects z onto the surface of a Riemann Sphere. If z is finite, it projects
 // to itself. If z is infinite, it projects to positive infinity on the real
 // axis.
 export fn projc128(z: c128) c128 =
 	if (!isinf(z)) z else (math::INF, math::copysignf64(0f64, z.1));
 
 // Returns the square root of z.
 export fn sqrtc128(z: c128) c128 = {
 	if (z.1 == 0f64) {
 		if (z.0 == 0f64) {
 			return (0f64, z.1);
 		};
 		if (z.0 < 0f64) {
 			return (0f64, math::copysignf64(math::sqrtf64(-z.0), z.1));
 		};
 		return (math::sqrtf64(z.0), z.1);
 	};
 	if (math::isinf(z.1)) {
 		return (math::INF, z.1);
 	};
 	if (z.0 == 0f64) {
 		if (z.1 < 0f64) {
 			const r = math::sqrtf64(-0.5 * z.1);
 			return (r, -r);
 		} else {
 			const r = math::sqrtf64(0.5 * z.1);
 			return (r, r);
 		};
 	};
 	let a = z.0, b = z.1;
 	const scale = if (math::absf64(a) > 4f64 || math::absf64(b) > 4f64) {
 		a *= 0.25;
 		b *= 0.25;
 		yield 2f64;
 	} else {
 		a *= 1.8014398509481984e16; // 2**54
 		b *= 1.8014398509481984e16;
 		yield 7.450580596923828125e-9; // 2**-27
 	};
 	let r = math::hypotf64(a, b);
 	const t = if (a > 0f64) {
 		const t = math::sqrtf64(0.5 * r + 0.5 * a);
 		r = scale * math::absf64(0.5 * b / t);
 		yield t * scale;
 	} else {
 		r = math::sqrtf64(0.5 * r - 0.5 * a);
 		const t = scale * math::absf64(0.5 * b / r);
 		r *= scale;
 		yield t;
 	};
 	return (t, if (b < 0f64) -r else r);
 };
 
 // Returns the sine of z, in radians.
 export fn sinc128(z: c128) c128 = {
 	if (z.1 == 0f64 && (math::isinf(z.0) || math::isnan(z.0))) {
 		return (math::NAN, z.1);
 	} else if (math::isinf(z.1)) {
 		if (z.0 == 0f64) {
 			return z;
 		} else if (math::isinf(z.0) || math::isnan(z.0)) {
 			return (math::NAN, z.1);
 		};
 	} else if (z.0 == 0f64 && math::isnan(z.1)) {
 		return z;
 	};
 	const shch = sinhcosh(z.1);
 	return (math::sinf64(z.0) * shch.1, math::cosf64(z.0) * shch.0);
 };
 
 // Returns the hyperbolic sine of z.
 export fn sinhc128(z: c128) c128 = {
 	if (z.0 == 0f64 && (math::isinf(z.1) || math::isnan(z.1))) {
 		return (z.0, math::NAN);
 	} else if (math::isinf(z.0)) {
 		if (z.1 == 0f64) {
 			return z;
 		} else if (math::isinf(z.1) || math::isnan(z.1)) {
 			return (z.0, math::NAN);
 		};
 	} else if (z.1 == 0f64 && math::isnan(z.0)) {
 		return (math::NAN, z.1);
 	};
 	const shch = sinhcosh(z.0);
 	return (math::cosf64(z.1) * shch.0, math::sinf64(z.1) * shch.1);
 };
 
 // Returns the arcsine, in radians, of z.
 export fn asinc128(z: c128) c128 = {
 	if (z.1 == 0f64 && math::absf64(z.0) <= 1f64) {
 		return (math::asinf64(z.0), z.1);
 	} else if (z.0 == 0f64 && math::absf64(z.1) <= 1f64) {
 		return (z.0, math::asinhf64(z.1));
 	} else if (math::isnan(z.1)) {
 		if (z.0 == 0f64) {
 			return (z.0, math::NAN);
 		} else if (math::isinf(z.0)) {
 			return (math::NAN, z.0);
 		} else {
 			return (math::NAN, math::NAN);
 		};
 	} else if (math::isinf(z.1)) {
 		if (math::isnan(z.0)) {
 			return z;
 		} else if (math::isinf(z.0)) {
 			return (math::copysignf64(math::PI / 4f64, z.0), z.1);
 		} else {
 			return (math::copysignf64(0f64, z.0), z.1);
 		};
 	} else if (math::isinf(z.0)) {
 		return (math::copysignf64(math::PI / 2f64, z.0),
 			math::copysignf64(z.0, z.1));
 	};
 	const ct = (-z.1, z.0); // i * z
 	const zz = mulc128(z, z);
 	const z1 = (1f64 - zz.0, -zz.1); // 1 - z * z
 	const z2 = sqrtc128(z1); // z2 = sqrt(1 - z * z)
 	const w = logc128(addc128(ct, z2));
 	return (w.1, -w.0); // -i * w
 };
 
 // Returns the inverse hyperbolic sine of z.
 export fn asinhc128(z: c128) c128 = {
 	if (z.1 == 0f64 && math::absf64(z.0) <= 1f64) {
 		return (math::asinhf64(z.0), z.1);
 	} else if (z.0 == 0f64 && math::absf64(z.1) <= 1f64) {
 		return (z.0, math::asinf64(z.1));
 	} else if (math::isinf(z.0)) {
 		if (math::isinf(z.1)) {
 			return (z.0, math::copysignf64(math::PI / 4f64, z.1));
 		} else if (math::isnan(z.1)) {
 			return z;
 		} else {
 			return (z.0, math::copysignf64(0f64, z.1));
 		};
 	} else if (math::isnan(z.0)) {
 		if (z.1 == 0f64) {
 			return z;
 		} else if (math::isinf(z.1)) {
 			return (z.1, z.0);
 		} else {
 			return (math::NAN, math::NAN);
 		};
 	} else if (math::isinf(z.1)) {
 		return (math::copysignf64(z.1, z.0),
 			math::copysignf64(math::PI / 2f64, z.1));
 	};
 	const zz = mulc128(z, z);
 	const z1 = (1f64 + zz.0, zz.1); // 1 + z * z
 	return logc128(addc128(z, sqrtc128(z1))); // log(x + sqrt(1 + x * x))
 };
 
 // Returns the cosine of z, in radians.
 export fn cosc128(z: c128) c128 = {
 	if (z.1 == 0f64 && (math::isinf(z.0) || math::isnan(z.0))) {
 		return (math::NAN, -z.1 * math::copysignf64(0f64, z.0));
 	} else if (math::isinf(z.1)) {
 		if (z.0 == 0f64) {
 			return (math::INF, -z.0 * math::copysignf64(0f64, z.1));
 		} else if (math::isinf(z.0) || math::isnan(z.0)) {
 			return (math::INF, math::NAN);
 		};
 	} else if (z.0 == 0f64 && math::isnan(z.1)) {
 		return (math::NAN, 0f64);
 	};
 	const shch = sinhcosh(z.1);
 	return (math::cosf64(z.0) * shch.1, -math::sinf64(z.0) * shch.0);
 };
 
 // Returns the hyperbolic cosine of z.
 export fn coshc128(z: c128) c128 = {
 	if (z.0 == 0f64 && (math::isinf(z.1) || math::isnan(z.1))) {
 		return (math::NAN, z.0 * math::copysignf64(0f64, z.1));
 	} else if (math::isinf(z.0)) {
 		if (z.1 == 0f64) {
 			return (math::INF, z.1 * math::copysignf64(0f64, z.0));
 		} else if (math::isinf(z.1) || math::isnan(z.1)) {
 			return (math::INF, math::NAN);
 		};
 	} else if (z.1 == 0f64 && math::isnan(z.0)) {
 		return (math::NAN, z.1);
 	};
 	const shch = sinhcosh(z.0);
 	return (math::cosf64(z.1) * shch.1, math::sinf64(z.1) * shch.0);
 };
 
 // Returns the arccosine, in radians, of z.
 export fn acosc128(z: c128) c128 = {
 	const w = asinc128(z);
 	return (math::PI / 2f64 - w.0, -w.1);
 };
 
 // Returns the inverse hyperbolic cosine of z.
 export fn acoshc128(z: c128) c128 = {
 	if (z.0 == 0f64 && z.1 == 0f64) {
 		return (0f64, math::copysignf64(math::PI / 2f64, z.1));
 	};
 	const w = acosc128(z);
 	if (w.1 <= 0f64) {
 		return (-w.1, w.0); // i * w
 	};
 	return (w.1, -w.0); // -i * w
 };
 
 fn sinhcosh(x: f64) (f64, f64) = {
 	if (math::absf64(x) <= 0.5) {
 		return (math::sinhf64(x), math::coshf64(x));
 	};
 	let e = math::expf64(x);
 	const ei = 0.5 / e;
 	e *= 0.5;
 	return (e - ei, e + ei);
 };
 
 // reducePi reduces the input argument x to the range (-Pi/2, Pi/2].
 // x must be greater than or equal to 0. For small arguments it
 // uses Cody-Waite reduction in 3 f64 parts based on:
 // "Elementary Function Evaluation: Algorithms and Implementation"
 // Jean-Michel Muller, 1997.
 // For very large arguments it uses Payne-Hanek range reduction based on:
 // "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit"
 // K. C. Ng et al, March 24, 1992.
 fn reducePi(x: f64) f64 = {
 	// reduceThreshold is the maximum value of x where the reduction using
 	// Cody-Waite reduction still gives accurate results. This threshold
 	// is set by t*PIn being representable as a f64 without error
 	// where t is given by t = floor(x * (1 / Pi)) and PIn are the leading partial
 	// terms of Pi. Since the leading terms, PI1 and PI2 below, have 30 and 32
 	// trailing zero bits respectively, t should have less than 30 significant bits.
 	//	t < 1<<30  -> floor(x*(1/Pi)+0.5) < 1<<30 -> x < (1<<30-1) * Pi - 0.5
 	// So, conservatively we can take x < 1<<30.
 	const reduceThreshold = (1u64 << 30): f64;
 	if (math::absf64(x) < reduceThreshold) {
 		// Use Cody-Waite reduction in three parts.
 		// PI1, PI2 and PI3 comprise an extended precision value of PI
 		// such that PI ~= PI1 + PI2 + PI3. The parts are chosen so
 		// that PI1 and PI2 have an approximately equal number of trailing
 		// zero bits. This ensures that t*PI1 and t*PI2 are exact for
 		// large integer values of t. The full precision PI3 ensures the
 		// approximation of PI is accurate to 102 bits to handle cancellation
 		// during subtraction.
 		const PI1 = 3.141592502593994;      // 0x400921fb40000000
 		const PI2 = 1.5099578831723193e-07; // 0x3e84442d00000000
 		const PI3 = 1.0780605716316238e-14; // 0x3d08469898cc5170
 		let t = x / math::PI;
 		t += 0.5;
 		t = (t: i64): f64;
 		return ((x - t*PI1) - t*PI2) - t*PI3;
 	};
 	// Must apply Payne-Hanek range reduction
 	const mask: u64 = 0x7FF;
 	const shift: u64 = 64 - 11 - 1;
 	const bias: u64 = 1023;
 	const fracMask: u64 = (1u64 << shift) - 1;
 
 	// Extract out the integer and exponent such that,
 	// x = ix * 2 ** exp.
 	let ix: u64 = math::f64bits(x);
 	let exp: u64 = (ix >> shift & mask) - bias - shift;
 	ix &= fracMask;
 	ix |= 1u64 << shift;
 
 	// mPi is the binary digits of 1/Pi as a u64 array,
 	// that is, 1/Pi = Sum mPi[i]*2^(-64*i).
 	// 19 64-bit digits give 1216 bits of precision
 	// to handle the largest possible float64 exponent.
 	let mPi: [_]u64 = [
 		0x0000000000000000,
 		0x517cc1b727220a94,
 		0xfe13abe8fa9a6ee0,
 		0x6db14acc9e21c820,
 		0xff28b1d5ef5de2b0,
 		0xdb92371d2126e970,
 		0x0324977504e8c90e,
 		0x7f0ef58e5894d39f,
 		0x74411afa975da242,
 		0x74ce38135a2fbf20,
 		0x9cc8eb1cc1a99cfa,
 		0x4e422fc5defc941d,
 		0x8ffc4bffef02cc07,
 		0xf79788c5ad05368f,
 		0xb69b3f6793e584db,
 		0xa7a31fb34f2ff516,
 		0xba93dd63f5f2f8bd,
 		0x9e839cfbc5294975,
 		0x35fdafd88fc6ae84,
 		0x2b0198237e3db5d5,
 	];
 	// Use the exponent to extract the 3 appropriate u64 digits from mPi,
 	// B ~ (z0, z1, z2), such that the product leading digit has the exponent -64.
 	// Note, exp >= 50 since x >= reduceThreshold and exp < 971 for maximum f64.
 	let digit: u64 = (exp + 64): u64 / 64;
 	let bitshift: u64 = (exp + 64): u64 % 64;
 	let z0: u64 = (mPi[digit] << bitshift) | (mPi[digit+1] >> (64 - bitshift));
 	let z1: u64 = (mPi[digit+1] << bitshift) | (mPi[digit+2] >> (64 - bitshift));
 	let z2: u64 = (mPi[digit+2] << bitshift) | (mPi[digit+3] >> (64 - bitshift));
 	// Multiply mantissa by the digits and extract the upper two digits (hi, lo).
 	let (z2hi, z2lo) = math::mulu64(z2, ix);
 	let (z1hi, z1lo) = math::mulu64(z1, ix);
 	let z0lo: u64 = z0 * ix;
 	let (lo, c) = checked::addu64(z1lo, z2hi);
 	let hi = z0lo + z1hi + (if (c) 1u64 else 0u64);
 	// Find the magnitude of the fraction.
 	let lz: u8 = math::leading_zeros_u64(hi);
 	let e: u64 = (bias - (lz + 1)): u64;
 	// Clear implicit mantissa bit and shift into place.
 	hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1)));
 	hi >>= 64 - shift;
 	// Include the exponent and convert to a float.
 	hi |= e << shift;
 	x = math::f64frombits(hi);
 	// map to (-Pi/2, Pi/2]
 	if (x > 0.5) {
 		x -= 1f64;
 	};
 	return math::PI * x;
 };
 
 // Taylor series expansion for cosh(2y) - cos(2x)
 fn tanSeries(z: c128) f64 = {
 	const MACHEP = 1f64/(1u64 << 53): f64;
 	let x = math::absf64(2f64 * z.0);
 	let y = math::absf64(2f64 * z.1);
 	x = reducePi(x);
 	x = x * x;
 	y = y * y;
 
 	let x2 = 1f64;
 	let y2 = 1f64;
 	let f = 1f64;
 	let rn = 0f64;
 	let d = 0f64;
 
 	for (true) {
 		rn += 1f64;
 		f *= rn;
 		rn += 1f64;
 		f *= rn;
 		x2 *= x;
 		y2 *= y;
 		let t = y2 + x2;
 		t /= f;
 		d += t;
 
 		rn += 1f64;
 		f *= rn;
 		rn += 1f64;
 		f *= rn;
 		x2 *= x;
 		y2 *= y;
 		t = y2 - x2;
 		t /= f;
 		d += t;
 		if (!(math::absf64(t/d) > MACHEP)) {
 			// Caution: Use ! and > instead of <= for correct behavior if t/d is NaN.
 			// See issue https://github.com/golang/go/issues/17577.
 			break;
 		};
 	};
 	return d;
 };
 
 // Returns the tangent of x.
 export fn tanc128(x: c128) c128 = {
 	if (math::isinf(x.1)) {
 		if (math::isinf(x.0) || math::isnan(x.0)) {
 			return (math::copysignf64(0f64, x.0), math::copysignf64(1f64, x.1));
 		};
 		return (math::copysignf64(0f64, math::sinf64(2f64*x.0)), math::copysignf64(1f64, x.1));
 	};
 	if (x.0 == 0f64 && math::isnan(x.1)) {
 		return x;
 	};
 	let d = math::cosf64(2f64*x.0) + math::coshf64(2f64*x.1);
 	if (math::absf64(d) < 0.25f64) {
 		d = tanSeries(x);
 	};
 	if (d == 0f64) {
 		return (math::INF, math::INF);
 	};
 	return (math::sinf64(2f64*x.0)/d, math::sinhf64(2f64*x.1)/d);
 };
 
 // Returns the hyperbolic tangent of x.
 export fn tanhc128(x: c128) c128 = {
 	if (math::isinf(x.0)) {
 		if (math::isinf(x.1) || math::isnan(x.1)) {
 			return (math::copysignf64(1f64, x.0), math::copysignf64(0f64, x.1));
 		};
 		return (math::copysignf64(1f64, x.0), math::copysignf64(0f64, math::sinf64(2f64*x.1)));
 	};
 	if (x.1 == 0f64 && math::isnan(x.0)) {
 		return x;
 	};
 	let d = math::coshf64(2f64*x.0) + math::cosf64(2f64*x.1);
 	if (d == 0f64) {
 		return (math::INF, math::INF);
 	};
 	return (math::sinhf64(2f64*x.0)/d, math::sinf64(2f64*x.1)/d);
 };
 
 // Returns the inverse tangent of x.
 export fn atanc128(x: c128) c128 = {
 	if (x.1 == 0f64) {
 		return (math::atanf64(x.0), x.1);
 	} else if (x.0 == 0f64 && math::absf64(x.1) <= 1f64) {
 		return (x.0, math::atanhf64(x.1));
 	} else if (math::isinf(x.1) || math::isinf(x.0)) {
 		if (math::isnan(x.0)) {
 			return (math::NAN, math::copysignf64(0f64, x.1));
 		};
 		return (math::copysignf64(math::PI/2f64, x.0), math::copysignf64(0f64, x.1));
 	} else if (math::isnan(x.0) || math::isnan(x.1)) {
 		return (math::NAN, math::NAN);
 	};
 	let x2 = x.0 * x.0;
 	let a = 1f64 - x2 - x.1 * x.1;
 	if (a == 0f64) {
 		return (math::NAN, math::NAN);
 	};
 	let t = 0.5f64 * math::atan2f64(2f64*x.0, a);
 	let w = reducePi(t);
 	t = x.1 - 1f64;
 	let b = x2 + t*t;
 	if (b == 0f64) {
 		return (math::NAN, math::NAN);
 	};
 	t = x.1 + 1f64;
 	let c = (x2 + t*t) / b;
 	return (w, 0.25f64*math::logf64(c));
 };
 
 // Returns the inverse hyperbolic tangent of x.
 export fn atanhc128(x: c128) c128 = {
 	let z = (-x.1, x.0); // z = i * x
 	z = atanc128(z);
 	return (z.1, -z.0); // z = -i * z
 };