diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java index bb9bd203e..4dc471172 100644 --- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java +++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java @@ -516,7 +516,7 @@ public Complex negate() { *

\( z \) projects to \( z \), except that all complex infinities (even those * with one infinite part and one NaN part) project to positive infinity on the real axis. * - * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to: + *

If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to: * *

return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());
* diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java index a624392de..e4d783d1b 100644 --- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java +++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java @@ -64,7 +64,7 @@ public final class ComplexFunctions { private static final double PI_OVER_4 = 0.25 * Math.PI; /** Natural logarithm of 2 (ln(2)). */ private static final double LN_2 = Math.log(2); - /** Base 10 logarithm of e divided by 2. */ + /** Base 10 logarithm of e divided by 2 (log10(e)/2). */ private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2; /** Base 10 logarithm of 2 (log10(2)). */ private static final double LOG10_2 = Math.log10(2); @@ -184,6 +184,42 @@ public final class ComplexFunctions { private ComplexFunctions() { } + /** + * Returns the absolute value of the complex number. This is also called complex norm, modulus, + * or magnitude. + * + *

\[ \text{abs}(x + i y) = \sqrt{(x^2 + y^2)} \] + * + *

Special cases: + * + *

+ * + *

The cases ensure that if either component is infinite then the result is positive + * infinity. If either component is NaN and this is not {@link #isInfinite(double, double) infinite} then + * the result is NaN. + * + *

This method follows the + * ISO C Standard, Annex G, + * in calculating the returned value without intermediate overflow or underflow. + * + *

The computed result will be within 1 ulp of the exact result. + * + * @param real Real part \( a \) of the complex number \( (a +ib) \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @return The absolute value. + * @see #isInfinite(double, double) + * @see #isNaN(double, double) + * @see Complex modulus + */ + public static double abs(double real, double imaginary) { + return computeAbs(real, imaginary); + } + /** * Returns the absolute value of the complex number. *

abs(x + i y) = sqrt(x^2 + y^2)
@@ -205,7 +241,7 @@ private ComplexFunctions() { * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). * @return The absolute value. */ - public static double abs(double real, double imaginary) { + private static double computeAbs(double real, double imaginary) { // Specialised implementation of hypot. // See NUMBERS-143 return hypot(real, imaginary); @@ -365,7 +401,7 @@ public static R negate(double real, double imaginary, ComplexSink action) *

\( z \) projects to \( z \), except that all complex infinities (even those * with one infinite part and one NaN part) project to positive infinity on the real axis. * - * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to: + *

If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to: * *

return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());
* @@ -388,7 +424,7 @@ public static R proj(double real, double imaginary, ComplexSink action) { /** * Returns the * - * exponential function of complex number using it's real and + * exponential function of complex number using its real and * imaginary part. * *

\[ \exp(z) = e^z \] @@ -452,7 +488,7 @@ public static R exp(double real, double imaginary, ComplexSink action) { zeroOrInf = real; } return action.apply(zeroOrInf * Math.cos(imaginary), - zeroOrInf * Math.sin(imaginary)); + zeroOrInf * Math.sin(imaginary)); } else if (Double.isNaN(real)) { // (NaN + i0) returns (NaN + i0) // (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception @@ -477,7 +513,7 @@ public static R exp(double real, double imaginary, ComplexSink action) { return action.apply(exp, imaginary); } return action.apply(exp * Math.cos(imaginary), - exp * Math.sin(imaginary)); + exp * Math.sin(imaginary)); } /** @@ -665,6 +701,61 @@ private static R log(DoubleUnaryOperator log, return action.apply(re, arg(real, imaginary)); } + /** + * Returns the + * + * square root of the complex number. + * + *

\[ \sqrt{x + iy} = \frac{1}{2} \sqrt{2} \left( \sqrt{ \sqrt{x^2 + y^2} + x } + i\ \text{sgn}(y) \sqrt{ \sqrt{x^2 + y^2} - x } \right) \] + * + *

The square root of \( z \) is in the range \( [0, +\infty) \) along the real axis and + * is unbounded along the imaginary axis. The imaginary part of the square root has a + * branch cut along the negative real axis \( (-infty,0) \). Special cases: + * + *

    + *
  • {@code z.conj().sqrt() == z.sqrt().conj()}. + *
  • If {@code z} is ±0 + i0, returns +0 + i0. + *
  • If {@code z} is x + i∞ for all x (including NaN), returns +∞ + i∞. + *
  • If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is −∞ + iy for finite positive-signed y, returns +0 + i∞. + *
  • If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0. + *
  • If {@code z} is −∞ + iNaN, returns NaN ± i∞ (where the sign of the imaginary part of the result is unspecified). + *
  • If {@code z} is +∞ + iNaN, returns +∞ + iNaN. + *
  • If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

Implements the following algorithm to compute \( \sqrt{x + iy} \): + *

    + *
  1. Let \( t = \sqrt{2 (|x| + |x + iy|)} \) + *
  2. if \( x \geq 0 \) return \( \frac{t}{2} + i \frac{y}{t} \) + *
  3. else return \( \frac{|y|}{t} + i\ \text{sgn}(y) \frac{t}{2} \) + *
+ * where: + *
    + *
  • \( |x| =\ \){@link Math#abs(double) abs}(x) + *
  • \( |x + y i| =\ \){@link Complex#abs} + *
  • \( \text{sgn}(y) =\ \){@link Math#copySign(double,double) copySign}(1.0, y) + *
+ * + *

The implementation is overflow and underflow safe based on the method described in:

+ *
+ * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994) + * Implementing complex elementary functions using exception handling. + * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244. + *
+ * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \). + * @param action Consumer for the square root of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see Sqrt + */ + public static R sqrt(double real, double imaginary, ComplexSink action) { + return computeSqrt(real, imaginary, action); + } + /** * Returns the square root of the complex number using its real * and imaginary parts {@code sqrt(x + i y)}. @@ -675,7 +766,7 @@ private static R log(DoubleUnaryOperator log, * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R sqrt(double real, double imaginary, ComplexSink action) { + private static R computeSqrt(double real, double imaginary, ComplexSink action) { // Handle NaN if (Double.isNaN(real) || Double.isNaN(imaginary)) { // Check for infinite @@ -862,6 +953,53 @@ public static R tan(double real, double imaginary, ComplexSink action) { return tanh(-imaginary, real, (r, i) -> multiplyNegativeI(r, i, action)); } + /** + * Returns the + * + * inverse sine of the complex number. + * + *

\[ \sin^{-1}(z) = - i \left(\ln{iz + \sqrt{1 - z^2}}\right) \] + * + *

The inverse sine of \( z \) is unbounded along the imaginary axis and + * in the range \( [-\pi, \pi] \) along the real axis. Special cases are handled + * as if the operation is implemented using \( \sin^{-1}(z) = -i \sinh^{-1}(iz) \). + * + *

The inverse sine is a multivalued function and requires a branch cut in + * the complex plane; the cut is conventionally placed at the line segments + * \( (\infty,-1) \) and \( (1,\infty) \) of the real axis. + * + *

This is implemented using real \( x \) and imaginary \( y \) parts: + * + *

\[ \begin{aligned} + * \sin^{-1}(z) &= \sin^{-1}(B) + i\ \text{sgn}(y)\ln \left(A + \sqrt{A^2-1} \right) \\ + * A &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\ + * B &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \] + * + *

where \( \text{sgn}(y) \) is the sign function implemented using + * {@link Math#copySign(double,double) copySign(1.0, y)}. + * + *

The implementation is based on the method described in:

+ *
+ * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997) + * Implementing the complex Arcsine and Arccosine Functions using Exception Handling. + * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335. + *
+ * + *

The code has been adapted from the Boost + * {@code c++} implementation {@code }. + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the inverse sine of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see ArcSin + */ + public static R asin(final double real, final double imaginary, + ComplexSink action) { + return computeAsin(real, imaginary, action); + } + /** * Returns the inverse sine of the complex number. * @@ -883,8 +1021,8 @@ public static R tan(double real, double imaginary, ComplexSink action) { * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R asin(final double real, final double imaginary, - ComplexSink action) { + private static R computeAsin(final double real, final double imaginary, + ComplexSink action) { // Compute with positive values and determine sign at the end final double x = Math.abs(real); final double y = Math.abs(imaginary); @@ -998,6 +1136,68 @@ public static R asin(final double real, final double imaginary, changeSign(im, imaginary)); } + /** + * Returns the + * + * inverse cosine of the complex number. + * + *

\[ \cos^{-1}(z) = \frac{\pi}{2} + i \left(\ln{iz + \sqrt{1 - z^2}}\right) \] + * + *

The inverse cosine of \( z \) is in the range \( [0, \pi) \) along the real axis and + * unbounded along the imaginary axis. Special cases: + * + *

    + *
  • {@code z.conj().acos() == z.acos().conj()}. + *
  • If {@code z} is ±0 + i0, returns π/2 − i0. + *
  • If {@code z} is ±0 + iNaN, returns π/2 + iNaN. + *
  • If {@code z} is x + i∞ for finite x, returns π/2 − i∞. + *
  • If {@code z} is x + iNaN, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is −∞ + iy for positive-signed finite y, returns π − i∞. + *
  • If {@code z} is +∞ + iy for positive-signed finite y, returns +0 − i∞. + *
  • If {@code z} is −∞ + i∞, returns 3π/4 − i∞. + *
  • If {@code z} is +∞ + i∞, returns π/4 − i∞. + *
  • If {@code z} is ±∞ + iNaN, returns NaN ± i∞ where the sign of the imaginary part of the result is unspecified. + *
  • If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + i∞, returns NaN − i∞. + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

The inverse cosine is a multivalued function and requires a branch cut in + * the complex plane; the cut is conventionally placed at the line segments + * \( (-\infty,-1) \) and \( (1,\infty) \) of the real axis. + * + *

This function is implemented using real \( x \) and imaginary \( y \) parts: + * + *

\[ \begin{aligned} + * \cos^{-1}(z) &= \cos^{-1}(B) - i\ \text{sgn}(y) \ln\left(A + \sqrt{A^2-1}\right) \\ + * A &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\ + * B &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \] + * + *

where \( \text{sgn}(y) \) is the sign function implemented using + * {@link Math#copySign(double,double) copySign(1.0, y)}. + * + *

The implementation is based on the method described in:

+ *
+ * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997) + * Implementing the complex Arcsine and Arccosine Functions using Exception Handling. + * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335. + *
+ * + *

The code has been adapted from the Boost + * {@code c++} implementation {@code }. + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the inverse cosine of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see ArcCos + */ + public static R acos(final double real, final double imaginary, + final ComplexSink action) { + return computeAcos(real, imaginary, action); + } + /** * Returns the inverse cosine of the complex number. * @@ -1019,7 +1219,7 @@ public static R asin(final double real, final double imaginary, * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R acos(final double real, final double imaginary, + private static R computeAcos(final double real, final double imaginary, final ComplexSink action) { // Compute with positive values and determine sign at the end final double x = Math.abs(real); @@ -1164,6 +1364,49 @@ public static R atan(double real, double imaginary, ComplexSink action) { return atanh(-imaginary, real, (r, i) -> multiplyNegativeI(r, i, action)); } + /** + * Returns the + * + * hyperbolic sine of the complexnumber. + * + *

\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \] + * + *

The hyperbolic sine of \( z \) is an entire function in the complex plane + * and is periodic with respect to the imaginary component with period \( 2\pi i \). + * Special cases: + * + *

    + *
  • {@code z.conj().sinh() == z.sinh().conj()}. + *
  • This is an odd function: \( \sinh(z) = -\sinh(-z) \). + *
  • If {@code z} is +0 + i0, returns +0 + i0. + *
  • If {@code z} is +0 + i∞, returns ±0 + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation). + *
  • If {@code z} is +0 + iNaN, returns ±0 + iNaN (where the sign of the real part of the result is unspecified). + *
  • If {@code z} is x + i∞ for positive finite x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is +∞ + i0, returns +∞ + i0. + *
  • If {@code z} is +∞ + iy for positive finite y, returns +∞ cis(y) (see ofCis(double) in Complex class). + *
  • If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation). + *
  • If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified). + *
  • If {@code z} is NaN + i0, returns NaN + i0. + *
  • If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

This is implemented using real \( x \) and imaginary \( y \) parts: + * + *

\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \] + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the hyperbolic sine of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action.. + * @see Sinh + */ + public static R sinh(double real, double imaginary, ComplexSink action) { + return computeSinh(real, imaginary, action); + } + /** * Returns the hyperbolic sine of the complex number using its real * and imaginary parts. @@ -1177,7 +1420,7 @@ public static R atan(double real, double imaginary, ComplexSink action) { * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R sinh(double real, double imaginary, ComplexSink action) { + public static R computeSinh(double real, double imaginary, ComplexSink action) { if (Double.isInfinite(real) && !Double.isFinite(imaginary)) { return action.apply(real, Double.NaN); } @@ -1205,7 +1448,50 @@ public static R sinh(double real, double imaginary, ComplexSink action) { } // No overflow of sinh/cosh return action.apply(Math.sinh(real) * Math.cos(imaginary), - Math.cosh(real) * Math.sin(imaginary)); + Math.cosh(real) * Math.sin(imaginary)); + } + + /** + * Returns the + * + * hyperbolic cosine of the complexnumber. + * + *

\[ \cosh(z) = \frac{1}{2} \left( e^{z} + e^{-z} \right) \] + * + *

The hyperbolic cosine of \( z \) is an entire function in the complex plane + * and is periodic with respect to the imaginary component with period \( 2\pi i \). + * Special cases: + * + *

    + *
  • {@code z.conj().cosh() == z.cosh().conj()}. + *
  • This is an even function: \( \cosh(z) = \cosh(-z) \). + *
  • If {@code z} is +0 + i0, returns 1 + i0. + *
  • If {@code z} is +0 + i∞, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified; "invalid" floating-point operation). + *
  • If {@code z} is +0 + iNaN, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified). + *
  • If {@code z} is x + i∞ for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is +∞ + i0, returns +∞ + i0. + *
  • If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y) (see ofCis(double) in Complex class). + *
  • If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified). + *
  • If {@code z} is +∞ + iNaN, returns +∞ + iNaN. + *
  • If {@code z} is NaN + i0, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified). + *
  • If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

This is implemented using real \( x \) and imaginary \( y \) parts: + * + *

\[ \cosh(x + iy) = \cosh(x)\cos(y) + i \sinh(x)\sin(y) \] + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the hyperbolic tangent of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see Cosh + */ + public static R cosh(double real, double imaginary, ComplexSink action) { + return computeCosh(real, imaginary, action); } /** @@ -1221,7 +1507,7 @@ public static R sinh(double real, double imaginary, ComplexSink action) { * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R cosh(double real, double imaginary, ComplexSink action) { + public static R computeCosh(double real, double imaginary, ComplexSink action) { // ISO C99: Preserve the even function by mapping to positive // f(z) = f(-z) if (Double.isInfinite(real) && !Double.isFinite(imaginary)) { @@ -1255,7 +1541,7 @@ public static R cosh(double real, double imaginary, ComplexSink action) { } // No overflow of sinh/cosh return action.apply(Math.cosh(real) * Math.cos(imaginary), - Math.sinh(real) * Math.sin(imaginary)); + Math.sinh(real) * Math.sin(imaginary)); } /** @@ -1319,6 +1605,58 @@ private static R coshsinh(double x, double real, double imaginary, boolean s return action.apply(re, im); } + /** + * Returns the + * + * hyperbolic tangent of the complexnumber. + * + *

\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \] + * + *

The hyperbolic tangent of \( z \) is an entire function in the complex plane + * and is periodic with respect to the imaginary component with period \( \pi i \) + * and has poles of the first order along the imaginary line, at coordinates + * \( (0, \pi(\frac{1}{2} + n)) \). + * Note that the {@code double} floating-point representation is unable to exactly represent + * \( \pi/2 \) and there is no value for which a pole error occurs. Special cases: + * + *

    + *
  • {@code z.conj().tanh() == z.tanh().conj()}. + *
  • This is an odd function: \( \tanh(z) = -\tanh(-z) \). + *
  • If {@code z} is +0 + i0, returns +0 + i0. + *
  • If {@code z} is 0 + i∞, returns 0 + iNaN. + *
  • If {@code z} is x + i∞ for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is 0 + iNaN, returns 0 + iNAN. + *
  • If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is +∞ + iy for positive-signed finite y, returns 1 + i0 sin(2y). + *
  • If {@code z} is +∞ + i∞, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified). + *
  • If {@code z} is +∞ + iNaN, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified). + *
  • If {@code z} is NaN + i0, returns NaN + i0. + *
  • If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

Special cases include the technical corrigendum + * + * DR 471: Complex math functions cacosh and ctanh. + * + *

This is defined using real \( x \) and imaginary \( y \) parts: + * + *

\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \] + * + *

The implementation uses double-angle identities to avoid overflow of {@code 2x} + * and {@code 2y}. + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the hyperbolic tangent of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see Tanh + */ + public static R tanh(double real, double imaginary, ComplexSink action) { + return computeTanh(real, imaginary, action); + } + /** * Returns the hyperbolic tangent of the complex number using its real * and imaginary parts. @@ -1332,7 +1670,7 @@ private static R coshsinh(double x, double real, double imaginary, boolean s * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R tanh(double real, double imaginary, ComplexSink action) { + public static R computeTanh(double real, double imaginary, ComplexSink action) { // Cache the absolute real value final double x = Math.abs(real); @@ -1424,7 +1762,7 @@ public static R tanh(double real, double imaginary, ComplexSink action) { final double cosy = Math.cos(imaginary); final double divisor = sinhx * sinhx + cosy * cosy; return action.apply(sinhx * coshx / divisor, - siny * cosy / divisor); + siny * cosy / divisor); } /** @@ -1545,6 +1883,59 @@ public static R acosh(double real, double imaginary, ComplexSink action) ); } + /** + * Returns the + * + * inverse hyperbolic tangent of the complex number. + * + *

\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \] + * + *

The inverse hyperbolic tangent of \( z \) is unbounded along the real axis and + * in the range \( [-\pi/2, \pi/2] \) along the imaginary axis. Special cases: + * + *

    + *
  • {@code z.conj().atanh() == z.atanh().conj()}. + *
  • This is an odd function: \( \tanh^{-1}(z) = -\tanh^{-1}(-z) \). + *
  • If {@code z} is +0 + i0, returns +0 + i0. + *
  • If {@code z} is +0 + iNaN, returns +0 + iNaN. + *
  • If {@code z} is +1 + i0, returns +∞ + i0 ("divide-by-zero" floating-point operation). + *
  • If {@code z} is x + i∞ for finite positive-signed x, returns +0 + iπ/2. + *
  • If {@code z} is x+iNaN for nonzero finite x, returns NaN+iNaN ("invalid" floating-point operation). + *
  • If {@code z} is +∞ + iy for finite positive-signed y, returns +0 + iπ/2. + *
  • If {@code z} is +∞ + i∞, returns +0 + iπ/2. + *
  • If {@code z} is +∞ + iNaN, returns +0 + iNaN. + *
  • If {@code z} is NaN+iy for finite y, returns NaN+iNaN ("invalid" floating-point operation). + *
  • If {@code z} is NaN + i∞, returns ±0 + iπ/2 (where the sign of the real part of the result is unspecified). + *
  • If {@code z} is NaN + iNaN, returns NaN + iNaN. + *
+ * + *

The inverse hyperbolic tangent is a multivalued function and requires a branch cut in + * the complex plane; the cut is conventionally placed at the line segments + * \( (\infty,-1] \) and \( [1,\infty) \) of the real axis. + * + *

This is implemented using real \( x \) and imaginary \( y \) parts: + * + *

\[ \tanh^{-1}(z) = \frac{1}{4} \ln \left(1 + \frac{4x}{(1-x)^2+y^2} \right) + \\ + * i \frac{1}{2} \left( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) + \frac{\pi}{2} \left(\text{sgn}(x^2+y^2-1)+1 \right) \text{sgn}(y) \right) \] + * + *

The imaginary part is computed using {@link Math#atan2(double, double)} to ensure the + * correct quadrant is returned from \( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) \). + * + *

The code has been adapted from the Boost + * {@code c++} implementation {@code }. + * + * @param real Real part \( a \) of the complex number \( (a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \). + * @param action Consumer for the inverse hyperbolic tangent of the complex number. + * @param the return type of the supplied action. + * @return the object returned by the supplied action. + * @see ArcTanh + */ + public static R atanh(final double real, final double imaginary, + final ComplexSink action) { + return computeAtanh(real, imaginary, action); + } + /** * Returns the inverse hyperbolic tangent of the complex number using its * real and imaginary parts. @@ -1567,8 +1958,8 @@ public static R acosh(double real, double imaginary, ComplexSink action) * @param the return type of the supplied action. * @return the object returned by the supplied action. */ - public static R atanh(final double real, final double imaginary, - final ComplexSink action) { + private static R computeAtanh(final double real, final double imaginary, + final ComplexSink action) { // Compute with positive values and determine sign at the end double x = Math.abs(real); double y = Math.abs(imaginary); diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java index 6d7196915..62dc190ee 100644 --- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java +++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java @@ -229,9 +229,9 @@ static void assertComplex(Complex c1, Complex c2, * @param maxUlps Maximum units of least precision between the two values. */ private static void assertOperation(String name, - UnaryOperator operation1, - ComplexUnaryOperator operation2, - long maxUlps) { + UnaryOperator operation1, + ComplexUnaryOperator operation2, + long maxUlps) { final List data = loadTestData(name); final long ulps = getTestUlps(maxUlps); for (final Complex[] pair : data) { diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java index c02d07bcd..2ac487b05 100644 --- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java +++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java @@ -218,31 +218,6 @@ private static void assertOperation(Complex c1, Complex c2, () -> String.format("%s expected: %s %s %s = %s", expectedName, c1, operationName, c2, z)); } - /** - * Assert the operation on the two complex numbers. - * - *

Assert the condition operation on the complex number is exactly equal to the condition - * operation on complex real and imaginary parts. - * - * @param c1 First input complex number. - * @param c2 Second input complex number. - * @param operation operation on the Complex objects. - * @param condition1 Condition operation on the Complex object. - * @param condition2 Condition operation on the complex real and imaginary parts. - * @param operationName the operation name - * @param expectedName the expected name - */ - private static void assertOperation(Complex c1, Complex c2, - BiFunction operation, - String operationName, - Predicate condition1, - TestUtils.DoubleBinaryPredicate condition2, - String expectedName) { - final Complex z = operation.apply(c1, c2); - Assertions.assertTrue(TestUtils.assertSame(z, expectedName, condition1, condition2), - () -> String.format("%s expected: %s %s %s = %s", expectedName, c1, operationName, c2, z)); - } - /** * Assert the operation on the complex number satisfies the conjugate equality. * @@ -278,8 +253,8 @@ private static void assertOperation(Complex c1, Complex c2, * @param operation2 Operation on the complex real and imaginary parts. */ private static void assertConjugateEquality(String name, - UnaryOperator operation1, - ComplexUnaryOperator operation2) { + UnaryOperator operation1, + ComplexUnaryOperator operation2) { // Edge cases. Inf/NaN are specifically handled in the C99 test cases // but are repeated here to enforce the conjugate equality even when the C99 // standard does not specify a sign. This may be revised in the future. @@ -322,9 +297,9 @@ private static void assertConjugateEquality(String name, * @param name Operation name. */ private static void assertConjugateEquality(Complex z, - String name, UnaryOperator operation1, - ComplexUnaryOperator operation2, - UnspecifiedSign sign) { + String name, UnaryOperator operation1, + ComplexUnaryOperator operation2, + UnspecifiedSign sign) { final Complex zConj = TestUtils.assertSame(z, "conj", Complex::conj, ComplexFunctions::conj); final Complex c1 = TestUtils.assertSame(zConj, name, operation1, operation2); @@ -362,9 +337,9 @@ private static void assertConjugateEquality(Complex z, * @param type the type */ private static void assertFunctionType(String name, - UnaryOperator operation1, - ComplexUnaryOperator operation2, - FunctionType type) { + UnaryOperator operation1, + ComplexUnaryOperator operation2, + FunctionType type) { // Note: It may not be possible to satisfy the conjugate equality // and be an odd/even function with regard to zero. // The C99 standard allows for these cases to have unspecified sign. @@ -412,10 +387,10 @@ private static void assertFunctionType(String name, * @param sign the sign specification */ private static void assertFunctionType(Complex z, - String name, - UnaryOperator operation1, - ComplexUnaryOperator operation2, - FunctionType type, UnspecifiedSign sign) { + String name, + UnaryOperator operation1, + ComplexUnaryOperator operation2, + FunctionType type, UnspecifiedSign sign) { final Complex c1 = TestUtils.assertSame(z, name, operation1, operation2); Complex c2 = TestUtils.assertSame(z.negate(), name, operation1, operation2); if (type == FunctionType.ODD) { @@ -454,9 +429,9 @@ private static void assertFunctionType(Complex z, * @param expected Expected complex number. */ private static void assertComplex(Complex z, - String name, UnaryOperator operation1, - ComplexUnaryOperator operation2, - Complex expected) { + String name, UnaryOperator operation1, + ComplexUnaryOperator operation2, + Complex expected) { assertComplex(z, name, operation1, operation2, expected, FunctionType.NONE, UnspecifiedSign.NONE); } @@ -481,9 +456,9 @@ private static void assertComplex(Complex z, * @param sign the sign specification */ private static void assertComplex(Complex z, - String name, UnaryOperator operation1, - ComplexUnaryOperator operation2, - Complex expected, UnspecifiedSign sign) { + String name, UnaryOperator operation1, + ComplexUnaryOperator operation2, + Complex expected, UnspecifiedSign sign) { assertComplex(z, name, operation1, operation2, expected, FunctionType.NONE, sign); } @@ -520,9 +495,9 @@ private static void assertComplex(Complex z, * @param sign Function sign specification. */ private static void assertComplex(Complex z, - String name, UnaryOperator operation1, - ComplexUnaryOperator operation2, - Complex expected, FunctionType type, UnspecifiedSign sign) { + String name, UnaryOperator operation1, + ComplexUnaryOperator operation2, + Complex expected, FunctionType type, UnspecifiedSign sign) { // Developer note: Set the sign specification to UnspecifiedSign.NONE // to see which equalities fail. They should be for input defined // in ISO C99 with an unspecified output sign, e.g. @@ -588,11 +563,11 @@ private static void assertComplex(Complex z, * @param type Function type. */ private static void assertComplex(Complex z, - String name, - UnaryOperator operation1, - ComplexUnaryOperator operation2, - Complex expected, - FunctionType type) { + String name, + UnaryOperator operation1, + ComplexUnaryOperator operation2, + Complex expected, + FunctionType type) { assertComplex(z, name, operation1, operation2, expected, type, UnspecifiedSign.NONE); } @@ -755,7 +730,7 @@ private static Complex complex(double real, double imaginary) { */ private static ArrayList createInfinites() { final double[] values = {0, 1, inf, negInf, nan}; - return createCombinations(values, "Inf", Complex::isInfinite, ComplexFunctions::isInfinite); + return createCombinations(values, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite); } /** @@ -785,7 +760,7 @@ private static ArrayList createZeroFinites() { */ private static ArrayList createNaNs() { final double[] values = {0, 1, nan}; - return createCombinations(values, "NaN", Complex::isNaN, ComplexFunctions::isNaN); + return createCombinations(values, "isNaN", Complex::isNaN, ComplexFunctions::isNaN); } /** @@ -869,18 +844,18 @@ void testMultiply() { // infinity, then the result of the * operator is an infinity;" for (final Complex z : infinites) { for (final Complex w : infinites) { - assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf"); } for (final Complex w : nonZeroFinites) { - assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf"); } // C.99 refers to non-zero finites. // Infer that Complex multiplication of zero with infinites is not defined. for (final Complex w : zeroFinites) { - assertOperation(z, w, Complex::multiply, "*", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); - assertOperation(w, z, Complex::multiply, "*", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); + assertOperation(z, w, Complex::multiply, "*", Complex::isNaN, "NaN"); + assertOperation(w, z, Complex::multiply, "*", Complex::isNaN, "NaN"); } } @@ -896,16 +871,16 @@ void testMultiply() { // Check multiply with (NaN,NaN) is not corrected final double[] values = {0, 1, inf, negInf, nan}; for (final Complex z : createCombinations(values, c -> true)) { - assertOperation(z, NAN, Complex::multiply, "*", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); - assertOperation(NAN, z, Complex::multiply, "*", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); + assertOperation(z, NAN, Complex::multiply, "*", Complex::isNaN, "NaN"); + assertOperation(NAN, z, Complex::multiply, "*", Complex::isNaN, "NaN"); } // Test multiply cases which result in overflow are corrected to infinity - assertOperation(maxMax, maxMax, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(maxNan, maxNan, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(nanMax, maxNan, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(maxNan, nanMax, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); - assertOperation(nanMax, nanMax, Complex::multiply, "*", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(maxMax, maxMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(maxNan, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(nanMax, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(maxNan, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(nanMax, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); } /** @@ -924,14 +899,14 @@ void testDivide() { // result of the / operator is an infinity;" for (final Complex z : infinites) { for (final Complex w : nonZeroFinites) { - assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); } for (final Complex w : zeroFinites) { - assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); } // Check inf/inf cannot be done for (final Complex w : infinites) { - assertOperation(z, w, Complex::divide, "/", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); + assertOperation(z, w, Complex::divide, "/", Complex::isNaN, "NaN"); } } @@ -947,10 +922,10 @@ void testDivide() { // a zero, then the result of the / operator is an infinity." for (final Complex w : zeroFinites) { for (final Complex z : nonZeroFinites) { - assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); } for (final Complex z : infinites) { - assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, ComplexFunctions::isInfinite, "Inf"); + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); } } @@ -968,10 +943,10 @@ void testDivide() { // Check (NaN,NaN) divide is not corrected for the edge case of divide by zero or infinite for (final Complex w : zeroFinites) { - assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); + assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN"); } for (final Complex w : infinites) { - assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, ComplexFunctions::isNaN, "NaN"); + assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN"); } } diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java index 1ea989da3..35ebd50b9 100644 --- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java +++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java @@ -442,9 +442,9 @@ void testNumberType() { */ private static void assertNumberType(double real, double imaginary, NumberType type) { final Complex z = Complex.ofCartesian(real, imaginary); - final boolean isNaN = TestUtils.assertSame(z, "NaN", Complex::isNaN, ComplexFunctions::isNaN); - final boolean isInfinite = TestUtils.assertSame(z, "Inf", Complex::isInfinite, ComplexFunctions::isInfinite); - final boolean isFinite = TestUtils.assertSame(z, "Fin", Complex::isFinite, ComplexFunctions::isFinite); + final boolean isNaN = TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN); + final boolean isInfinite = TestUtils.assertSame(z, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite); + final boolean isFinite = TestUtils.assertSame(z, "isFinite", Complex::isFinite, ComplexFunctions::isFinite); // A number can be only one int count = isNaN ? 1 : 0; count += isInfinite ? 1 : 0; @@ -478,7 +478,7 @@ void testConjugate() { @Test void testConjugateNaN() { final Complex z = TestUtils.assertSame(NAN, "conj", Complex::conj, ComplexFunctions::conj); - Assertions.assertTrue(TestUtils.assertSame(z, "NaN", Complex::isNaN, ComplexFunctions::isNaN)); + Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN)); } @Test @@ -503,7 +503,7 @@ void testNegate() { @Test void testNegateNaN() { final Complex z = TestUtils.assertSame(NAN, "negate", Complex::negate, ComplexFunctions::negate); - Assertions.assertTrue(TestUtils.assertSame(z, "NaN", Complex::isNaN, ComplexFunctions::isNaN)); + Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN)); } @Test @@ -853,7 +853,7 @@ void testMultiply() { void testMultiplyInfInf() { final Complex z = infInf.multiply(infInf); // Assert.assertTrue(z.isNaN()); // MATH-620 - Assertions.assertTrue(TestUtils.assertSame(z, "Inf", Complex::isInfinite, ComplexFunctions::isInfinite)); + Assertions.assertTrue(TestUtils.assertSame(z, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite)); // Expected results from g++: Assertions.assertEquals(Complex.ofCartesian(nan, inf), infInf.multiply(infInf)); @@ -1117,7 +1117,7 @@ void testDivideZeroZero() { void testDivideNaN() { final Complex x = Complex.ofCartesian(3.0, 4.0); final Complex z = x.divide(NAN); - Assertions.assertTrue(TestUtils.assertSame(z, "NaN", Complex::isNaN, ComplexFunctions::isNaN)); + Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN)); } @Test