1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.numbers.complex;
19
20 import java.io.Serializable;
21 import java.util.ArrayList;
22 import java.util.List;
23 import java.util.function.DoubleUnaryOperator;
24
25 /**
26 * Cartesian representation of a complex number. The complex number is expressed
27 * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \)
28 * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the
29 * complex number \( a + ib \), \( a \) is called the <em>real part</em> and
30 * \( b \) is called the <em>imaginary part</em>.
31 *
32 * <p>This class is immutable. All arithmetic will create a new instance for the
33 * result.</p>
34 *
35 * <p>Arithmetic in this class conforms to the C99 standard for complex numbers
36 * defined in ISO/IEC 9899, Annex G. Methods have been named using the equivalent
37 * method in ISO C99. The behavior for special cases is listed as defined in C99.</p>
38 *
39 * <p>For functions \( f \) which obey the conjugate equality \( conj(f(z)) = f(conj(z)) \),
40 * the specifications for the upper half-plane imply the specifications for the lower
41 * half-plane.</p>
42 *
43 * <p>For functions that are either odd, \( f(z) = -f(-z) \), or even, \( f(z) = f(-z) \),
44 * the specifications for the first quadrant imply the specifications for the other three
45 * quadrants.</p>
46 *
47 * <p>Special cases of <a href="http://mathworld.wolfram.com/BranchCut.html">branch cuts</a>
48 * for multivalued functions adopt the principle value convention from C99. Specials cases
49 * from C99 that raise the "invalid" or "divide-by-zero"
50 * <a href="https://en.cppreference.com/w/c/numeric/fenv/FE_exceptions">floating-point
51 * exceptions</a> return the documented value without an explicit mechanism to notify
52 * of the exception case, that is no exceptions are thrown during computations in-line with
53 * the convention of the corresponding single-valued functions in
54 * {@link Math Math}.
55 * These cases are documented in the method special cases as "invalid" or "divide-by-zero"
56 * floating-point operation.
57 * Note: Invalid floating-point exception cases will result in a complex number where the
58 * cardinality of NaN component parts has increased as a real or imaginary part could
59 * not be computed and is set to NaN.
60 *
61 * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards">
62 * ISO/IEC 9899 - Programming languages - C</a>
63 */
64 public final class Complex implements Serializable {
65 /**
66 * A complex number representing \( i \), the square root of \( -1 \).
67 *
68 * <p>\( (0 + i 1) \).
69 */
70 public static final Complex I = new Complex(0, 1);
71 /**
72 * A complex number representing one.
73 *
74 * <p>\( (1 + i 0) \).
75 */
76 public static final Complex ONE = new Complex(1, 0);
77 /**
78 * A complex number representing zero.
79 *
80 * <p>\( (0 + i 0) \).
81 */
82 public static final Complex ZERO = new Complex(0, 0);
83
84 /** A complex number representing {@code NaN + i NaN}. */
85 private static final Complex NAN = new Complex(Double.NaN, Double.NaN);
86 /** π/2. */
87 private static final double PI_OVER_2 = 0.5 * Math.PI;
88 /** π/4. */
89 private static final double PI_OVER_4 = 0.25 * Math.PI;
90 /** Natural logarithm of 2 (ln(2)). */
91 private static final double LN_2 = Math.log(2);
92 /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */
93 private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2;
94 /** Base 10 logarithm of 2 (log10(2)). */
95 private static final double LOG10_2 = Math.log10(2);
96 /** {@code 1/2}. */
97 private static final double HALF = 0.5;
98 /** {@code sqrt(2)}. */
99 private static final double ROOT2 = 1.4142135623730951;
100 /** {@code 1.0 / sqrt(2)}.
101 * This is pre-computed to the closest double from the exact result.
102 * It is 1 ULP different from 1.0 / Math.sqrt(2) but equal to Math.sqrt(2) / 2.
103 */
104 private static final double ONE_OVER_ROOT2 = 0.7071067811865476;
105 /** The bit representation of {@code -0.0}. */
106 private static final long NEGATIVE_ZERO_LONG_BITS = Double.doubleToLongBits(-0.0);
107 /** Exponent offset in IEEE754 representation. */
108 private static final int EXPONENT_OFFSET = 1023;
109 /**
110 * Largest double-precision floating-point number such that
111 * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
112 * bound on the relative error due to rounding real numbers to double
113 * precision floating-point numbers.
114 *
115 * <p>In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
116 * Copied from o.a.c.numbers.Precision.
117 *
118 * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
119 */
120 private static final double EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
121 /** Mask to remove the sign bit from a long. */
122 private static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL;
123 /** Mask to extract the 52-bit mantissa from a long representation of a double. */
124 private static final long MANTISSA_MASK = 0x000f_ffff_ffff_ffffL;
125 /** The multiplier used to split the double value into hi and low parts. This must be odd
126 * and a value of 2^s + 1 in the range {@code p/2 <= s <= p-1} where p is the number of
127 * bits of precision of the floating point number. Here {@code s = 27}.*/
128 private static final double MULTIPLIER = 1.34217729E8;
129
130 /**
131 * Crossover point to switch computation for asin/acos factor A.
132 * This has been updated from the 1.5 value used by Hull et al to 10
133 * as used in boost::math::complex.
134 * @see <a href="https://svn.boost.org/trac/boost/ticket/7290">Boost ticket 7290</a>
135 */
136 private static final double A_CROSSOVER = 10.0;
137 /** Crossover point to switch computation for asin/acos factor B. */
138 private static final double B_CROSSOVER = 0.6471;
139 /**
140 * The safe maximum double value {@code x} to avoid loss of precision in asin/acos.
141 * Equal to sqrt(M) / 8 in Hull, et al (1997) with M the largest normalised floating-point value.
142 */
143 private static final double SAFE_MAX = Math.sqrt(Double.MAX_VALUE) / 8;
144 /**
145 * The safe minimum double value {@code x} to avoid loss of precision/underflow in asin/acos.
146 * Equal to sqrt(u) * 4 in Hull, et al (1997) with u the smallest normalised floating-point value.
147 */
148 private static final double SAFE_MIN = Math.sqrt(Double.MIN_NORMAL) * 4;
149 /**
150 * The safe maximum double value {@code x} to avoid loss of precision in atanh.
151 * Equal to sqrt(M) / 2 with M the largest normalised floating-point value.
152 */
153 private static final double SAFE_UPPER = Math.sqrt(Double.MAX_VALUE) / 2;
154 /**
155 * The safe minimum double value {@code x} to avoid loss of precision/underflow in atanh.
156 * Equal to sqrt(u) * 2 with u the smallest normalised floating-point value.
157 */
158 private static final double SAFE_LOWER = Math.sqrt(Double.MIN_NORMAL) * 2;
159 /** The safe maximum double value {@code x} to avoid overflow in sqrt. */
160 private static final double SQRT_SAFE_UPPER = Double.MAX_VALUE / 8;
161 /**
162 * A safe maximum double value {@code m} where {@code e^m} is not infinite.
163 * This can be used when functions require approximations of sinh(x) or cosh(x)
164 * when x is large using exp(x):
165 * <pre>
166 * sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
167 * cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2 </pre>
168 *
169 * <p>This value can be used to approximate e^x using a product:
170 *
171 * <pre>
172 * e^x = product_n (e^m) * e^(x-nm)
173 * n = (int) x/m
174 * e.g. e^2000 = e^m * e^m * e^(2000 - 2m) </pre>
175 *
176 * <p>The value should be below ln(max_value) ~ 709.783.
177 * The value m is set to an integer for less error when subtracting m and chosen as
178 * even (m=708) as it is used as a threshold in tanh with m/2.
179 *
180 * <p>The value is used to compute e^x multiplied by a small number avoiding
181 * overflow (sinh/cosh) or a small number divided by e^x without underflow due to
182 * infinite e^x (tanh). The following conditions are used:
183 * <pre>
184 * 0.5 * e^m * Double.MIN_VALUE * e^m * e^m = Infinity
185 * 2.0 / e^m / e^m = 0.0 </pre>
186 */
187 private static final double SAFE_EXP = 708;
188 /**
189 * The value of Math.exp(SAFE_EXP): e^708.
190 * To be used in overflow/underflow safe products of e^m to approximate e^x where x > m.
191 */
192 private static final double EXP_M = Math.exp(SAFE_EXP);
193
194 /** 54 shifted 20-bits to align with the exponent of the upper 32-bits of a double. */
195 private static final int EXP_54 = 0x36_00000;
196 /** Represents an exponent of 500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */
197 private static final int EXP_500 = 0x5f3_00000;
198 /** Represents an exponent of 1024 in unbiased form (infinite or nan)
199 * shifted 20-bits to align with the upper 32-bits of a double. */
200 private static final int EXP_1024 = 0x7ff_00000;
201 /** Represents an exponent of -500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */
202 private static final int EXP_NEG_500 = 0x20b_00000;
203 /** 2^600. */
204 private static final double TWO_POW_600 = 0x1.0p+600;
205 /** 2^-600. */
206 private static final double TWO_POW_NEG_600 = 0x1.0p-600;
207
208 /** Serializable version identifier. */
209 private static final long serialVersionUID = 20180201L;
210
211 /**
212 * The size of the buffer for {@link #toString()}.
213 *
214 * <p>The longest double will require a sign, a maximum of 17 digits, the decimal place
215 * and the exponent, e.g. for max value this is 24 chars: -1.7976931348623157e+308.
216 * Set the buffer size to twice this and round up to a power of 2 thus
217 * allowing for formatting characters. The size is 64.
218 */
219 private static final int TO_STRING_SIZE = 64;
220 /** The minimum number of characters in the format. This is 5, e.g. {@code "(0,0)"}. */
221 private static final int FORMAT_MIN_LEN = 5;
222 /** {@link #toString() String representation}. */
223 private static final char FORMAT_START = '(';
224 /** {@link #toString() String representation}. */
225 private static final char FORMAT_END = ')';
226 /** {@link #toString() String representation}. */
227 private static final char FORMAT_SEP = ',';
228 /** The minimum number of characters before the separator. This is 2, e.g. {@code "(0"}. */
229 private static final int BEFORE_SEP = 2;
230
231 /** The imaginary part. */
232 private final double imaginary;
233 /** The real part. */
234 private final double real;
235
236 /**
237 * Define a constructor for a Complex.
238 * This is used in functions that implement trigonomic identities.
239 */
240 @FunctionalInterface
241 private interface ComplexConstructor {
242 /**
243 * Create a complex number given the real and imaginary parts.
244 *
245 * @param real Real part.
246 * @param imaginary Imaginary part.
247 * @return {@code Complex} object.
248 */
249 Complex create(double real, double imaginary);
250 }
251
252 /**
253 * Private default constructor.
254 *
255 * @param real Real part.
256 * @param imaginary Imaginary part.
257 */
258 private Complex(double real, double imaginary) {
259 this.real = real;
260 this.imaginary = imaginary;
261 }
262
263 /**
264 * Create a complex number given the real and imaginary parts.
265 *
266 * @param real Real part.
267 * @param imaginary Imaginary part.
268 * @return {@code Complex} number.
269 */
270 public static Complex ofCartesian(double real, double imaginary) {
271 return new Complex(real, imaginary);
272 }
273
274 /**
275 * Creates a complex number from its polar representation using modulus {@code rho} (\( \rho \))
276 * and phase angle {@code theta} (\( \theta \)).
277 *
278 * \[ \begin{aligned}
279 * x &= \rho \cos(\theta) \\
280 * y &= \rho \sin(\theta) \end{aligned} \]
281 *
282 * <p>Requires that {@code rho} is non-negative and non-NaN and {@code theta} is finite;
283 * otherwise returns a complex with NaN real and imaginary parts. A {@code rho} value of
284 * {@code -0.0} is considered negative and an invalid modulus.
285 *
286 * <p>A non-NaN complex number constructed using this method will satisfy the following
287 * to within floating-point error when {@code theta} is in the range
288 * \( -\pi\ \lt \theta \leq \pi \):
289 *
290 * <pre>
291 * Complex.ofPolar(rho, theta).abs() == rho
292 * Complex.ofPolar(rho, theta).arg() == theta</pre>
293 *
294 * <p>If {@code rho} is infinite then the resulting parts may be infinite or NaN
295 * following the rules for double arithmetic, for example:</p>
296 *
297 * <ul>
298 * <li>{@code ofPolar(}\( -0.0 \){@code , }\( 0 \){@code ) = }\( \text{NaN} + i \text{NaN} \)
299 * <li>{@code ofPolar(}\( 0.0 \){@code , }\( 0 \){@code ) = }\( 0 + i 0 \)
300 * <li>{@code ofPolar(}\( 1 \){@code , }\( 0 \){@code ) = }\( 1 + i 0 \)
301 * <li>{@code ofPolar(}\( 1 \){@code , }\( \pi \){@code ) = }\( -1 + i \sin(\pi) \)
302 * <li>{@code ofPolar(}\( \infty \){@code , }\( \pi \){@code ) = }\( -\infty + i \infty \)
303 * <li>{@code ofPolar(}\( \infty \){@code , }\( 0 \){@code ) = }\( -\infty + i \text{NaN} \)
304 * <li>{@code ofPolar(}\( \infty \){@code , }\( -\frac{\pi}{4} \){@code ) = }\( \infty - i \infty \)
305 * <li>{@code ofPolar(}\( \infty \){@code , }\( 5\frac{\pi}{4} \){@code ) = }\( -\infty - i \infty \)
306 * </ul>
307 *
308 * <p>This method is the functional equivalent of the C++ method {@code std::polar}.
309 *
310 * @param rho The modulus of the complex number.
311 * @param theta The argument of the complex number.
312 * @return {@code Complex} number.
313 * @see <a href="http://mathworld.wolfram.com/PolarCoordinates.html">Polar Coordinates</a>
314 */
315 public static Complex ofPolar(double rho, double theta) {
316 // Require finite theta and non-negative, non-nan rho
317 if (!Double.isFinite(theta) || negative(rho) || Double.isNaN(rho)) {
318 return NAN;
319 }
320 final double x = rho * Math.cos(theta);
321 final double y = rho * Math.sin(theta);
322 return new Complex(x, y);
323 }
324
325 /**
326 * Create a complex cis number. This is also known as the complex exponential:
327 *
328 * \[ \text{cis}(x) = e^{ix} = \cos(x) + i \sin(x) \]
329 *
330 * @param x {@code double} to build the cis number.
331 * @return {@code Complex} cis number.
332 * @see <a href="http://mathworld.wolfram.com/Cis.html">Cis</a>
333 */
334 public static Complex ofCis(double x) {
335 return new Complex(Math.cos(x), Math.sin(x));
336 }
337
338 /**
339 * Returns a {@code Complex} instance representing the specified string {@code s}.
340 *
341 * <p>If {@code s} is {@code null}, then a {@code NullPointerException} is thrown.
342 *
343 * <p>The string must be in a format compatible with that produced by
344 * {@link #toString() Complex.toString()}.
345 * The format expects a start and end parentheses surrounding two numeric parts split
346 * by a separator. Leading and trailing spaces are allowed around each numeric part.
347 * Each numeric part is parsed using {@link Double#parseDouble(String)}. The parts
348 * are interpreted as the real and imaginary parts of the complex number.
349 *
350 * <p>Examples of valid strings and the equivalent {@code Complex} are shown below:
351 *
352 * <pre>
353 * "(0,0)" = Complex.ofCartesian(0, 0)
354 * "(0.0,0.0)" = Complex.ofCartesian(0, 0)
355 * "(-0.0, 0.0)" = Complex.ofCartesian(-0.0, 0)
356 * "(-1.23, 4.56)" = Complex.ofCartesian(-1.23, 4.56)
357 * "(1e300,-1.1e-2)" = Complex.ofCartesian(1e300, -1.1e-2)</pre>
358 *
359 * @param s String representation.
360 * @return {@code Complex} number.
361 * @throws NullPointerException if the string is null.
362 * @throws NumberFormatException if the string does not contain a parsable complex number.
363 * @see Double#parseDouble(String)
364 * @see #toString()
365 */
366 public static Complex parse(String s) {
367 final int len = s.length();
368 if (len < FORMAT_MIN_LEN) {
369 throw new NumberFormatException(
370 parsingExceptionMsg("Input too short, expected format",
371 FORMAT_START + "x" + FORMAT_SEP + "y" + FORMAT_END, s));
372 }
373
374 // Confirm start: '('
375 if (s.charAt(0) != FORMAT_START) {
376 throw new NumberFormatException(
377 parsingExceptionMsg("Expected start delimiter", FORMAT_START, s));
378 }
379
380 // Confirm end: ')'
381 if (s.charAt(len - 1) != FORMAT_END) {
382 throw new NumberFormatException(
383 parsingExceptionMsg("Expected end delimiter", FORMAT_END, s));
384 }
385
386 // Confirm separator ',' is between at least 2 characters from
387 // either end: "(x,x)"
388 // Count back from the end ignoring the last 2 characters.
389 final int sep = s.lastIndexOf(FORMAT_SEP, len - 3);
390 if (sep < BEFORE_SEP) {
391 throw new NumberFormatException(
392 parsingExceptionMsg("Expected separator between two numbers", FORMAT_SEP, s));
393 }
394
395 // Should be no more separators
396 if (s.indexOf(FORMAT_SEP, sep + 1) != -1) {
397 throw new NumberFormatException(
398 parsingExceptionMsg("Incorrect number of parts, expected only 2 using separator",
399 FORMAT_SEP, s));
400 }
401
402 // Try to parse the parts
403
404 final String rePart = s.substring(1, sep);
405 final double re;
406 try {
407 re = Double.parseDouble(rePart);
408 } catch (final NumberFormatException ex) {
409 throw new NumberFormatException(
410 parsingExceptionMsg("Could not parse real part", rePart, s));
411 }
412
413 final String imPart = s.substring(sep + 1, len - 1);
414 final double im;
415 try {
416 im = Double.parseDouble(imPart);
417 } catch (final NumberFormatException ex) {
418 throw new NumberFormatException(
419 parsingExceptionMsg("Could not parse imaginary part", imPart, s));
420 }
421
422 return ofCartesian(re, im);
423 }
424
425 /**
426 * Creates an exception message.
427 *
428 * @param message Message prefix.
429 * @param error Input that caused the error.
430 * @param s String representation.
431 * @return A message.
432 */
433 private static String parsingExceptionMsg(String message,
434 Object error,
435 String s) {
436 final StringBuilder sb = new StringBuilder(100)
437 .append(message)
438 .append(" '").append(error)
439 .append("' for input \"").append(s).append('"');
440 return sb.toString();
441 }
442
443 /**
444 * Gets the real part \( a \) of this complex number \( (a + i b) \).
445 *
446 * @return The real part.
447 */
448 public double getReal() {
449 return real;
450 }
451
452 /**
453 * Gets the real part \( a \) of this complex number \( (a + i b) \).
454 *
455 * <p>This method is the equivalent of the C++ method {@code std::complex::real}.
456 *
457 * @return The real part.
458 * @see #getReal()
459 */
460 public double real() {
461 return getReal();
462 }
463
464 /**
465 * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
466 *
467 * @return The imaginary part.
468 */
469 public double getImaginary() {
470 return imaginary;
471 }
472
473 /**
474 * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
475 *
476 * <p>This method is the equivalent of the C++ method {@code std::complex::imag}.
477 *
478 * @return The imaginary part.
479 * @see #getImaginary()
480 */
481 public double imag() {
482 return getImaginary();
483 }
484
485 /**
486 * Returns the absolute value of this complex number. This is also called complex norm, modulus,
487 * or magnitude.
488 *
489 * <p>\[ \text{abs}(x + i y) = \sqrt{(x^2 + y^2)} \]
490 *
491 * <p>Special cases:
492 *
493 * <ul>
494 * <li>{@code abs(x + iy) == abs(y + ix) == abs(x - iy)}.
495 * <li>If {@code z} is ±∞ + iy for any y, returns +∞.
496 * <li>If {@code z} is x + iNaN for non-infinite x, returns NaN.
497 * <li>If {@code z} is x + i0, returns |x|.
498 * </ul>
499 *
500 * <p>The cases ensure that if either component is infinite then the result is positive
501 * infinity. If either component is NaN and this is not {@link #isInfinite() infinite} then
502 * the result is NaN.
503 *
504 * <p>This method follows the
505 * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
506 * in calculating the returned value without intermediate overflow or underflow.
507 *
508 * <p>The computed result will be within 1 ulp of the exact result.
509 *
510 * @return The absolute value.
511 * @see #isInfinite()
512 * @see #isNaN()
513 * @see <a href="http://mathworld.wolfram.com/ComplexModulus.html">Complex modulus</a>
514 */
515 public double abs() {
516 return abs(real, imaginary);
517 }
518
519 /**
520 * Returns the absolute value of the complex number.
521 * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre>
522 *
523 * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3:
524 * "The hypot functions compute the square root of the sum of the squares of x and y,
525 * without undue overflow or underflow."
526 *
527 * <ul>
528 * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent.
529 * <li>hypot(x, ±0) is equivalent to |x|.
530 * <li>hypot(±∞, y) returns +∞, even if y is a NaN.
531 * </ul>
532 *
533 * <p>This method is called by all methods that require the absolute value of the complex
534 * number, e.g. abs(), sqrt() and log().
535 *
536 * @param real Real part.
537 * @param imaginary Imaginary part.
538 * @return The absolute value.
539 */
540 private static double abs(double real, double imaginary) {
541 // Specialised implementation of hypot.
542 // See NUMBERS-143
543 return hypot(real, imaginary);
544 }
545
546 /**
547 * Returns the argument of this complex number.
548 *
549 * <p>The argument is the angle phi between the positive real axis and
550 * the point representing this number in the complex plane.
551 * The value returned is between \( -\pi \) (not inclusive)
552 * and \( \pi \) (inclusive), with negative values returned for numbers with
553 * negative imaginary parts.
554 *
555 * <p>If either real or imaginary part (or both) is NaN, then the result is NaN.
556 * Infinite parts are handled as {@linkplain Math#atan2} handles them,
557 * essentially treating finite parts as zero in the presence of an
558 * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on
559 * the signs of the infinite parts.
560 *
561 * <p>This code follows the
562 * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
563 * in calculating the returned value using the {@code atan2(y, x)} method for complex
564 * \( x + iy \).
565 *
566 * @return The argument of this complex number.
567 * @see Math#atan2(double, double)
568 */
569 public double arg() {
570 // Delegate
571 return Math.atan2(imaginary, real);
572 }
573
574 /**
575 * Returns the squared norm value of this complex number. This is also called the absolute
576 * square.
577 *
578 * <p>\[ \text{norm}(x + i y) = x^2 + y^2 \]
579 *
580 * <p>If either component is infinite then the result is positive infinity. If either
581 * component is NaN and this is not {@link #isInfinite() infinite} then the result is NaN.
582 *
583 * <p>Note: This method may not return the same value as the square of {@link #abs()} as
584 * that method uses an extended precision computation.
585 *
586 * <p>{@code norm()} can be used as a faster alternative than {@code abs()} for ranking by
587 * magnitude. If used for ranking any overflow to infinity will create an equal ranking for
588 * values that may be still distinguished by {@code abs()}.
589 *
590 * @return The square norm value.
591 * @see #isInfinite()
592 * @see #isNaN()
593 * @see #abs()
594 * @see <a href="http://mathworld.wolfram.com/AbsoluteSquare.html">Absolute square</a>
595 */
596 public double norm() {
597 if (isInfinite()) {
598 return Double.POSITIVE_INFINITY;
599 }
600 return real * real + imaginary * imaginary;
601 }
602
603 /**
604 * Returns {@code true} if either the real <em>or</em> imaginary component of the complex number is NaN
605 * <em>and</em> the complex number is not infinite.
606 *
607 * <p>Note that:
608 * <ul>
609 * <li>There is more than one complex number that can return {@code true}.
610 * <li>Different representations of NaN can be distinguished by the
611 * {@link #equals(Object) Complex.equals(Object)} method.
612 * </ul>
613 *
614 * @return {@code true} if this instance contains NaN and no infinite parts.
615 * @see Double#isNaN(double)
616 * @see #isInfinite()
617 * @see #equals(Object) Complex.equals(Object)
618 */
619 public boolean isNaN() {
620 if (Double.isNaN(real) || Double.isNaN(imaginary)) {
621 return !isInfinite();
622 }
623 return false;
624 }
625
626 /**
627 * Returns {@code true} if either real or imaginary component of the complex number is infinite.
628 *
629 * <p>Note: A complex number with at least one infinite part is regarded
630 * as an infinity (even if its other part is a NaN).
631 *
632 * @return {@code true} if this instance contains an infinite value.
633 * @see Double#isInfinite(double)
634 */
635 public boolean isInfinite() {
636 return Double.isInfinite(real) || Double.isInfinite(imaginary);
637 }
638
639 /**
640 * Returns {@code true} if both real and imaginary component of the complex number are finite.
641 *
642 * @return {@code true} if this instance contains finite values.
643 * @see Double#isFinite(double)
644 */
645 public boolean isFinite() {
646 return Double.isFinite(real) && Double.isFinite(imaginary);
647 }
648
649 /**
650 * Returns the
651 * <a href="http://mathworld.wolfram.com/ComplexConjugate.html">conjugate</a>
652 * \( \overline{z} \) of this complex number \( z \).
653 *
654 * <p>\[ \begin{aligned}
655 * z &= a + i b \\
656 * \overline{z} &= a - i b \end{aligned}\]
657 *
658 * @return The conjugate (\( \overline{z} \)) of this complex number.
659 */
660 public Complex conj() {
661 return new Complex(real, -imaginary);
662 }
663
664 /**
665 * Returns a {@code Complex} whose value is the negation of both the real and imaginary parts
666 * of complex number \( z \).
667 *
668 * <p>\[ \begin{aligned}
669 * z &= a + i b \\
670 * -z &= -a - i b \end{aligned} \]
671 *
672 * @return \( -z \).
673 */
674 public Complex negate() {
675 return new Complex(-real, -imaginary);
676 }
677
678 /**
679 * Returns the projection of this complex number onto the Riemann sphere.
680 *
681 * <p>\( z \) projects to \( z \), except that all complex infinities (even those
682 * with one infinite part and one NaN part) project to positive infinity on the real axis.
683 *
684 * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to:
685 *
686 * <pre>return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());</pre>
687 *
688 * @return \( z \) projected onto the Riemann sphere.
689 * @see #isInfinite()
690 * @see <a href="http://pubs.opengroup.org/onlinepubs/9699919799/functions/cproj.html">
691 * IEEE and ISO C standards: cproj</a>
692 */
693 public Complex proj() {
694 if (isInfinite()) {
695 return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0.0, imaginary));
696 }
697 return this;
698 }
699
700 /**
701 * Returns a {@code Complex} whose value is {@code (this + addend)}.
702 * Implements the formula:
703 *
704 * <p>\[ (a + i b) + (c + i d) = (a + c) + i (b + d) \]
705 *
706 * @param addend Value to be added to this complex number.
707 * @return {@code this + addend}.
708 * @see <a href="http://mathworld.wolfram.com/ComplexAddition.html">Complex Addition</a>
709 */
710 public Complex add(Complex addend) {
711 return new Complex(real + addend.real,
712 imaginary + addend.imaginary);
713 }
714
715 /**
716 * Returns a {@code Complex} whose value is {@code (this + addend)},
717 * with {@code addend} interpreted as a real number.
718 * Implements the formula:
719 *
720 * <p>\[ (a + i b) + c = (a + c) + i b \]
721 *
722 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
723 * real-only and complex numbers.</p>
724 *
725 * <p>Note: This method preserves the sign of the imaginary component \( b \) if it is {@code -0.0}.
726 * The sign would be lost if adding \( (c + i 0) \) using
727 * {@link #add(Complex) add(Complex.ofCartesian(addend, 0))} since
728 * {@code -0.0 + 0.0 = 0.0}.
729 *
730 * @param addend Value to be added to this complex number.
731 * @return {@code this + addend}.
732 * @see #add(Complex)
733 * @see #ofCartesian(double, double)
734 */
735 public Complex add(double addend) {
736 return new Complex(real + addend, imaginary);
737 }
738
739 /**
740 * Returns a {@code Complex} whose value is {@code (this + addend)},
741 * with {@code addend} interpreted as an imaginary number.
742 * Implements the formula:
743 *
744 * <p>\[ (a + i b) + i d = a + i (b + d) \]
745 *
746 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
747 * imaginary-only and complex numbers.</p>
748 *
749 * <p>Note: This method preserves the sign of the real component \( a \) if it is {@code -0.0}.
750 * The sign would be lost if adding \( (0 + i d) \) using
751 * {@link #add(Complex) add(Complex.ofCartesian(0, addend))} since
752 * {@code -0.0 + 0.0 = 0.0}.
753 *
754 * @param addend Value to be added to this complex number.
755 * @return {@code this + addend}.
756 * @see #add(Complex)
757 * @see #ofCartesian(double, double)
758 */
759 public Complex addImaginary(double addend) {
760 return new Complex(real, imaginary + addend);
761 }
762
763 /**
764 * Returns a {@code Complex} whose value is {@code (this - subtrahend)}.
765 * Implements the formula:
766 *
767 * <p>\[ (a + i b) - (c + i d) = (a - c) + i (b - d) \]
768 *
769 * @param subtrahend Value to be subtracted from this complex number.
770 * @return {@code this - subtrahend}.
771 * @see <a href="http://mathworld.wolfram.com/ComplexSubtraction.html">Complex Subtraction</a>
772 */
773 public Complex subtract(Complex subtrahend) {
774 return new Complex(real - subtrahend.real,
775 imaginary - subtrahend.imaginary);
776 }
777
778 /**
779 * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
780 * with {@code subtrahend} interpreted as a real number.
781 * Implements the formula:
782 *
783 * <p>\[ (a + i b) - c = (a - c) + i b \]
784 *
785 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
786 * real-only and complex numbers.</p>
787 *
788 * @param subtrahend Value to be subtracted from this complex number.
789 * @return {@code this - subtrahend}.
790 * @see #subtract(Complex)
791 */
792 public Complex subtract(double subtrahend) {
793 return new Complex(real - subtrahend, imaginary);
794 }
795
796 /**
797 * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
798 * with {@code subtrahend} interpreted as an imaginary number.
799 * Implements the formula:
800 *
801 * <p>\[ (a + i b) - i d = a + i (b - d) \]
802 *
803 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
804 * imaginary-only and complex numbers.</p>
805 *
806 * @param subtrahend Value to be subtracted from this complex number.
807 * @return {@code this - subtrahend}.
808 * @see #subtract(Complex)
809 */
810 public Complex subtractImaginary(double subtrahend) {
811 return new Complex(real, imaginary - subtrahend);
812 }
813
814 /**
815 * Returns a {@code Complex} whose value is {@code (minuend - this)},
816 * with {@code minuend} interpreted as a real number.
817 * Implements the formula:
818 * \[ c - (a + i b) = (c - a) - i b \]
819 *
820 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
821 * real-only and complex numbers.</p>
822 *
823 * <p>Note: This method inverts the sign of the imaginary component \( b \) if it is {@code 0.0}.
824 * The sign would not be inverted if subtracting from \( c + i 0 \) using
825 * {@link #subtract(Complex) Complex.ofCartesian(minuend, 0).subtract(this)} since
826 * {@code 0.0 - 0.0 = 0.0}.
827 *
828 * @param minuend Value this complex number is to be subtracted from.
829 * @return {@code minuend - this}.
830 * @see #subtract(Complex)
831 * @see #ofCartesian(double, double)
832 */
833 public Complex subtractFrom(double minuend) {
834 return new Complex(minuend - real, -imaginary);
835 }
836
837 /**
838 * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
839 * with {@code minuend} interpreted as an imaginary number.
840 * Implements the formula:
841 * \[ i d - (a + i b) = -a + i (d - b) \]
842 *
843 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
844 * imaginary-only and complex numbers.</p>
845 *
846 * <p>Note: This method inverts the sign of the real component \( a \) if it is {@code 0.0}.
847 * The sign would not be inverted if subtracting from \( 0 + i d \) using
848 * {@link #subtract(Complex) Complex.ofCartesian(0, minuend).subtract(this)} since
849 * {@code 0.0 - 0.0 = 0.0}.
850 *
851 * @param minuend Value this complex number is to be subtracted from.
852 * @return {@code this - subtrahend}.
853 * @see #subtract(Complex)
854 * @see #ofCartesian(double, double)
855 */
856 public Complex subtractFromImaginary(double minuend) {
857 return new Complex(-real, minuend - imaginary);
858 }
859
860 /**
861 * Returns a {@code Complex} whose value is {@code this * factor}.
862 * Implements the formula:
863 *
864 * <p>\[ (a + i b)(c + i d) = (ac - bd) + i (ad + bc) \]
865 *
866 * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
867 *
868 * @param factor Value to be multiplied by this complex number.
869 * @return {@code this * factor}.
870 * @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
871 */
872 public Complex multiply(Complex factor) {
873 return multiply(real, imaginary, factor.real, factor.imaginary);
874 }
875
876 /**
877 * Returns a {@code Complex} whose value is:
878 * <pre>
879 * (a + i b)(c + i d) = (ac - bd) + i (ad + bc)</pre>
880 *
881 * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
882 *
883 * @param re1 Real component of first number.
884 * @param im1 Imaginary component of first number.
885 * @param re2 Real component of second number.
886 * @param im2 Imaginary component of second number.
887 * @return (a + b i)(c + d i).
888 */
889 private static Complex multiply(double re1, double im1, double re2, double im2) {
890 double a = re1;
891 double b = im1;
892 double c = re2;
893 double d = im2;
894 final double ac = a * c;
895 final double bd = b * d;
896 final double ad = a * d;
897 final double bc = b * c;
898 double x = ac - bd;
899 double y = ad + bc;
900
901 // --------------
902 // NaN can occur if:
903 // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
904 // - a multiplication of infinity by zero (ac,bd,ad,bc).
905 // - a subtraction of infinity from infinity (e.g. ac - bd)
906 // Note that (ac,bd,ad,bc) can be infinite due to overflow.
907 //
908 // Detect a NaN result and perform correction.
909 //
910 // Modification from the listing in ISO C99 G.5.1 (6)
911 // Do not correct infinity multiplied by zero. This is left as NaN.
912 // --------------
913
914 if (Double.isNaN(x) && Double.isNaN(y)) {
915 // Recover infinities that computed as NaN+iNaN ...
916 boolean recalc = false;
917 if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
918 isNotZero(c, d)) {
919 // This complex is infinite.
920 // "Box" the infinity and change NaNs in the other factor to 0.
921 a = boxInfinity(a);
922 b = boxInfinity(b);
923 c = changeNaNtoZero(c);
924 d = changeNaNtoZero(d);
925 recalc = true;
926 }
927 if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
928 isNotZero(a, b)) {
929 // The other complex is infinite.
930 // "Box" the infinity and change NaNs in the other factor to 0.
931 c = boxInfinity(c);
932 d = boxInfinity(d);
933 a = changeNaNtoZero(a);
934 b = changeNaNtoZero(b);
935 recalc = true;
936 }
937 if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
938 Double.isInfinite(ad) || Double.isInfinite(bc))) {
939 // The result overflowed to infinity.
940 // Recover infinities from overflow by changing NaNs to 0 ...
941 a = changeNaNtoZero(a);
942 b = changeNaNtoZero(b);
943 c = changeNaNtoZero(c);
944 d = changeNaNtoZero(d);
945 recalc = true;
946 }
947 if (recalc) {
948 x = Double.POSITIVE_INFINITY * (a * c - b * d);
949 y = Double.POSITIVE_INFINITY * (a * d + b * c);
950 }
951 }
952 return new Complex(x, y);
953 }
954
955 /**
956 * Box values for the real or imaginary component of an infinite complex number.
957 * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
958 * The sign is maintained.
959 *
960 * <pre>
961 * inf = 1
962 * -inf = -1
963 * x = 0
964 * -x = -0
965 * </pre>
966 *
967 * @param component the component
968 * @return The boxed value
969 */
970 private static double boxInfinity(double component) {
971 return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
972 }
973
974 /**
975 * Checks if the complex number is not zero.
976 *
977 * @param real the real component
978 * @param imaginary the imaginary component
979 * @return true if the complex is not zero
980 */
981 private static boolean isNotZero(double real, double imaginary) {
982 // The use of equals is deliberate.
983 // This method must distinguish NaN from zero thus ruling out:
984 // (real != 0.0 || imaginary != 0.0)
985 return !(real == 0.0 && imaginary == 0.0);
986 }
987
988 /**
989 * Change NaN to zero preserving the sign; otherwise return the value.
990 *
991 * @param value the value
992 * @return The new value
993 */
994 private static double changeNaNtoZero(double value) {
995 return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
996 }
997
998 /**
999 * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
1000 * interpreted as a real number.
1001 * Implements the formula:
1002 *
1003 * <p>\[ (a + i b) c = (ac) + i (bc) \]
1004 *
1005 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1006 * real-only and complex numbers.</p>
1007 *
1008 * <p>Note: This method should be preferred over using
1009 * {@link #multiply(Complex) multiply(Complex.ofCartesian(factor, 0))}. Multiplication
1010 * can generate signed zeros if either {@code this} complex has zeros for the real
1011 * and/or imaginary component, or if the factor is zero. The summation of signed zeros
1012 * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
1013 * from the equivalent call to multiply by a real-only number.
1014 *
1015 * @param factor Value to be multiplied by this complex number.
1016 * @return {@code this * factor}.
1017 * @see #multiply(Complex)
1018 */
1019 public Complex multiply(double factor) {
1020 return new Complex(real * factor, imaginary * factor);
1021 }
1022
1023 /**
1024 * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
1025 * interpreted as an imaginary number.
1026 * Implements the formula:
1027 *
1028 * <p>\[ (a + i b) id = (-bd) + i (ad) \]
1029 *
1030 * <p>This method can be used to compute the multiplication of this complex number \( z \)
1031 * by \( i \) using a factor with magnitude 1.0. This should be used in preference to
1032 * {@link #multiply(Complex) multiply(Complex.I)} with or without {@link #negate() negation}:</p>
1033 *
1034 * \[ \begin{aligned}
1035 * iz &= (-b + i a) \\
1036 * -iz &= (b - i a) \end{aligned} \]
1037 *
1038 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1039 * imaginary-only and complex numbers.</p>
1040 *
1041 * <p>Note: This method should be preferred over using
1042 * {@link #multiply(Complex) multiply(Complex.ofCartesian(0, factor))}. Multiplication
1043 * can generate signed zeros if either {@code this} complex has zeros for the real
1044 * and/or imaginary component, or if the factor is zero. The summation of signed zeros
1045 * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
1046 * from the equivalent call to multiply by an imaginary-only number.
1047 *
1048 * @param factor Value to be multiplied by this complex number.
1049 * @return {@code this * factor}.
1050 * @see #multiply(Complex)
1051 */
1052 public Complex multiplyImaginary(double factor) {
1053 return new Complex(-imaginary * factor, real * factor);
1054 }
1055
1056 /**
1057 * Returns a {@code Complex} whose value is {@code (this / divisor)}.
1058 * Implements the formula:
1059 *
1060 * <p>\[ \frac{a + i b}{c + i d} = \frac{(ac + bd) + i (bc - ad)}{c^2+d^2} \]
1061 *
1062 * <p>Re-calculates NaN result values to recover infinities as specified in C99 standard G.5.1.
1063 *
1064 * @param divisor Value by which this complex number is to be divided.
1065 * @return {@code this / divisor}.
1066 * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
1067 */
1068 public Complex divide(Complex divisor) {
1069 return divide(real, imaginary, divisor.real, divisor.imaginary);
1070 }
1071
1072 /**
1073 * Returns a {@code Complex} whose value is:
1074 * <pre>
1075 * <code>
1076 * a + i b (ac + bd) + i (bc - ad)
1077 * ------- = -----------------------
1078 * c + i d c<sup>2</sup> + d<sup>2</sup>
1079 * </code>
1080 * </pre>
1081 *
1082 * <p>Recalculates to recover infinities as specified in C99
1083 * standard G.5.1. Method is fully in accordance with
1084 * C++11 standards for complex numbers.</p>
1085 *
1086 * <p>Note: In the event of divide by zero this method produces the same result
1087 * as dividing by a real-only zero using {@link #divide(double)}.
1088 *
1089 * @param re1 Real component of first number.
1090 * @param im1 Imaginary component of first number.
1091 * @param re2 Real component of second number.
1092 * @param im2 Imaginary component of second number.
1093 * @return (a + i b) / (c + i d).
1094 * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
1095 * @see #divide(double)
1096 */
1097 private static Complex divide(double re1, double im1, double re2, double im2) {
1098 double a = re1;
1099 double b = im1;
1100 double c = re2;
1101 double d = im2;
1102 int ilogbw = 0;
1103 // Get the exponent to scale the divisor parts to the range [1, 2).
1104 final int exponent = getScale(c, d);
1105 if (exponent <= Double.MAX_EXPONENT) {
1106 ilogbw = exponent;
1107 c = Math.scalb(c, -ilogbw);
1108 d = Math.scalb(d, -ilogbw);
1109 }
1110 final double denom = c * c + d * d;
1111
1112 // Note: Modification from the listing in ISO C99 G.5.1 (8):
1113 // Avoid overflow if a or b are very big.
1114 // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
1115 // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
1116 // (bc - ad) with large negative values.
1117 // Use the maximum exponent as an approximation to the magnitude.
1118 if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
1119 ilogbw -= 2;
1120 a /= 4;
1121 b /= 4;
1122 }
1123
1124 double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
1125 double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
1126 // Recover infinities and zeros that computed as NaN+iNaN
1127 // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
1128 if (Double.isNaN(x) && Double.isNaN(y)) {
1129 if (denom == 0.0 &&
1130 (!Double.isNaN(a) || !Double.isNaN(b))) {
1131 // nonzero/zero
1132 // This case produces the same result as divide by a real-only zero
1133 // using Complex.divide(+/-0.0)
1134 x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
1135 y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
1136 } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
1137 Double.isFinite(c) && Double.isFinite(d)) {
1138 // infinite/finite
1139 a = boxInfinity(a);
1140 b = boxInfinity(b);
1141 x = Double.POSITIVE_INFINITY * (a * c + b * d);
1142 y = Double.POSITIVE_INFINITY * (b * c - a * d);
1143 } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
1144 Double.isFinite(a) && Double.isFinite(b)) {
1145 // finite/infinite
1146 c = boxInfinity(c);
1147 d = boxInfinity(d);
1148 x = 0.0 * (a * c + b * d);
1149 y = 0.0 * (b * c - a * d);
1150 }
1151 }
1152 return new Complex(x, y);
1153 }
1154
1155 /**
1156 * Returns a {@code Complex} whose value is {@code (this / divisor)},
1157 * with {@code divisor} interpreted as a real number.
1158 * Implements the formula:
1159 *
1160 * <p>\[ \frac{a + i b}{c} = \frac{a}{c} + i \frac{b}{c} \]
1161 *
1162 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1163 * real-only and complex numbers.</p>
1164 *
1165 * <p>Note: This method should be preferred over using
1166 * {@link #divide(Complex) divide(Complex.ofCartesian(divisor, 0))}. Division
1167 * can generate signed zeros if {@code this} complex has zeros for the real
1168 * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
1169 * in {@link #divide(Complex)} may create zeros in the result that differ in sign
1170 * from the equivalent call to divide by a real-only number.
1171 *
1172 * @param divisor Value by which this complex number is to be divided.
1173 * @return {@code this / divisor}.
1174 * @see #divide(Complex)
1175 */
1176 public Complex divide(double divisor) {
1177 return new Complex(real / divisor, imaginary / divisor);
1178 }
1179
1180 /**
1181 * Returns a {@code Complex} whose value is {@code (this / divisor)},
1182 * with {@code divisor} interpreted as an imaginary number.
1183 * Implements the formula:
1184 *
1185 * <p>\[ \frac{a + i b}{id} = \frac{b}{d} - i \frac{a}{d} \]
1186 *
1187 * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1188 * imaginary-only and complex numbers.</p>
1189 *
1190 * <p>Note: This method should be preferred over using
1191 * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))}. Division
1192 * can generate signed zeros if {@code this} complex has zeros for the real
1193 * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
1194 * in {@link #divide(Complex)} may create zeros in the result that differ in sign
1195 * from the equivalent call to divide by an imaginary-only number.
1196 *
1197 * <p>Warning: This method will generate a different result from
1198 * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))} if the divisor is zero.
1199 * In this case the divide method using a zero-valued Complex will produce the same result
1200 * as dividing by a real-only zero. The output from dividing by imaginary zero will create
1201 * infinite and NaN values in the same component parts as the output from
1202 * {@code this.divide(Complex.ZERO).multiplyImaginary(1)}, however the sign
1203 * of some infinite values may be negated.
1204 *
1205 * @param divisor Value by which this complex number is to be divided.
1206 * @return {@code this / divisor}.
1207 * @see #divide(Complex)
1208 * @see #divide(double)
1209 */
1210 public Complex divideImaginary(double divisor) {
1211 return new Complex(imaginary / divisor, -real / divisor);
1212 }
1213
1214 /**
1215 * Returns the
1216 * <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
1217 * exponential function</a> of this complex number.
1218 *
1219 * <p>\[ \exp(z) = e^z \]
1220 *
1221 * <p>The exponential function of \( z \) is an entire function in the complex plane.
1222 * Special cases:
1223 *
1224 * <ul>
1225 * <li>{@code z.conj().exp() == z.exp().conj()}.
1226 * <li>If {@code z} is ±0 + i0, returns 1 + i0.
1227 * <li>If {@code z} is x + i∞ for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1228 * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1229 * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
1230 * <li>If {@code z} is −∞ + iy for finite y, returns +0 cis(y) (see {@link #ofCis(double)}).
1231 * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y).
1232 * <li>If {@code z} is −∞ + i∞, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
1233 * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
1234 * <li>If {@code z} is −∞ + iNaN, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
1235 * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
1236 * <li>If {@code z} is NaN + i0, returns NaN + i0.
1237 * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
1238 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1239 * </ul>
1240 *
1241 * <p>Implements the formula:
1242 *
1243 * <p>\[ \exp(x + iy) = e^x (\cos(y) + i \sin(y)) \]
1244 *
1245 * @return The exponential of this complex number.
1246 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
1247 */
1248 public Complex exp() {
1249 if (Double.isInfinite(real)) {
1250 // Set the scale factor applied to cis(y)
1251 double zeroOrInf;
1252 if (real < 0) {
1253 if (!Double.isFinite(imaginary)) {
1254 // (−∞ + i∞) or (−∞ + iNaN) returns (±0 ± i0) (where the signs of the
1255 // real and imaginary parts of the result are unspecified).
1256 // Here we preserve the conjugate equality.
1257 return new Complex(0, Math.copySign(0, imaginary));
1258 }
1259 // (−∞ + iy) returns +0 cis(y), for finite y
1260 zeroOrInf = 0;
1261 } else {
1262 // (+∞ + i0) returns +∞ + i0.
1263 if (imaginary == 0) {
1264 return this;
1265 }
1266 // (+∞ + i∞) or (+∞ + iNaN) returns (±∞ + iNaN) and raises the invalid
1267 // floating-point exception (where the sign of the real part of the
1268 // result is unspecified).
1269 if (!Double.isFinite(imaginary)) {
1270 return new Complex(real, Double.NaN);
1271 }
1272 // (+∞ + iy) returns (+∞ cis(y)), for finite nonzero y.
1273 zeroOrInf = real;
1274 }
1275 return new Complex(zeroOrInf * Math.cos(imaginary),
1276 zeroOrInf * Math.sin(imaginary));
1277 } else if (Double.isNaN(real)) {
1278 // (NaN + i0) returns (NaN + i0)
1279 // (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception
1280 // (NaN + iNaN) returns (NaN + iNaN)
1281 return imaginary == 0 ? this : NAN;
1282 } else if (!Double.isFinite(imaginary)) {
1283 // (x + i∞) or (x + iNaN) returns (NaN + iNaN) and raises the invalid
1284 // floating-point exception, for finite x.
1285 return NAN;
1286 }
1287 // real and imaginary are finite.
1288 // Compute e^a * (cos(b) + i sin(b)).
1289
1290 // Special case:
1291 // (±0 + i0) returns (1 + i0)
1292 final double exp = Math.exp(real);
1293 if (imaginary == 0) {
1294 return new Complex(exp, imaginary);
1295 }
1296 return new Complex(exp * Math.cos(imaginary),
1297 exp * Math.sin(imaginary));
1298 }
1299
1300 /**
1301 * Returns the
1302 * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
1303 * natural logarithm</a> of this complex number.
1304 *
1305 * <p>The natural logarithm of \( z \) is unbounded along the real axis and
1306 * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
1307 * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \).
1308 * Special cases:
1309 *
1310 * <ul>
1311 * <li>{@code z.conj().log() == z.log().conj()}.
1312 * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation).
1313 * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation).
1314 * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
1315 * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1316 * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ.
1317 * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
1318 * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
1319 * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
1320 * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
1321 * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1322 * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
1323 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1324 * </ul>
1325 *
1326 * <p>Implements the formula:
1327 *
1328 * <p>\[ \ln(z) = \ln |z| + i \arg(z) \]
1329 *
1330 * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
1331 *
1332 * <p>The implementation is based on the method described in:</p>
1333 * <blockquote>
1334 * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
1335 * Implementing complex elementary functions using exception handling.
1336 * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
1337 * </blockquote>
1338 *
1339 * @return The natural logarithm of this complex number.
1340 * @see Math#log(double)
1341 * @see #abs()
1342 * @see #arg()
1343 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
1344 */
1345 public Complex log() {
1346 return log(Math::log, HALF, LN_2, Complex::ofCartesian);
1347 }
1348
1349 /**
1350 * Returns the base 10
1351 * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
1352 * common logarithm</a> of this complex number.
1353 *
1354 * <p>The common logarithm of \( z \) is unbounded along the real axis and
1355 * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
1356 * common logarithm has a branch cut along the negative real axis \( (-infty,0] \).
1357 * Special cases are as defined in the {@link #log() natural logarithm}:
1358 *
1359 * <p>Implements the formula:
1360 *
1361 * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
1362 *
1363 * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
1364 *
1365 * @return The base 10 logarithm of this complex number.
1366 * @see Math#log10(double)
1367 * @see #abs()
1368 * @see #arg()
1369 */
1370 public Complex log10() {
1371 return log(Math::log10, LOG_10E_O_2, LOG10_2, Complex::ofCartesian);
1372 }
1373
1374 /**
1375 * Returns the logarithm of this complex number using the provided function.
1376 * Implements the formula:
1377 *
1378 * <pre>
1379 * log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre>
1380 *
1381 * <p>Warning: The argument {@code logOf2} must be equal to {@code log(2)} using the
1382 * provided log function otherwise scaling using powers of 2 in the case of overflow
1383 * will be incorrect. This is provided as an internal optimisation.
1384 *
1385 * @param log Log function.
1386 * @param logOfeOver2 The log function applied to e, then divided by 2.
1387 * @param logOf2 The log function applied to 2.
1388 * @param constructor Constructor for the returned complex.
1389 * @return The logarithm of this complex number.
1390 * @see #abs()
1391 * @see #arg()
1392 */
1393 private Complex log(DoubleUnaryOperator log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
1394 // Handle NaN
1395 if (Double.isNaN(real) || Double.isNaN(imaginary)) {
1396 // Return NaN unless infinite
1397 if (isInfinite()) {
1398 return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
1399 }
1400 // No-use of the input constructor
1401 return NAN;
1402 }
1403
1404 // Returns the real part:
1405 // log(sqrt(x^2 + y^2))
1406 // log(x^2 + y^2) / 2
1407
1408 // Compute with positive values
1409 double x = Math.abs(real);
1410 double y = Math.abs(imaginary);
1411
1412 // Find the larger magnitude.
1413 if (x < y) {
1414 final double tmp = x;
1415 x = y;
1416 y = tmp;
1417 }
1418
1419 if (x == 0) {
1420 // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
1421 return constructor.create(Double.NEGATIVE_INFINITY,
1422 negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
1423 }
1424
1425 double re;
1426
1427 // This alters the implementation of Hull et al (1994) which used a standard
1428 // precision representation of |z|: sqrt(x*x + y*y).
1429 // This formula should use the same definition of the magnitude returned
1430 // by Complex.abs() which is a high precision computation with scaling.
1431 // The checks for overflow thus only require ensuring the output of |z|
1432 // will not overflow or underflow.
1433
1434 if (x > HALF && x < ROOT2) {
1435 // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
1436 re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
1437 } else {
1438 // Check for over/underflow in |z|
1439 // When scaling:
1440 // log(a / b) = log(a) - log(b)
1441 // So initialize the result with the log of the scale factor.
1442 re = 0;
1443 if (x > Double.MAX_VALUE / 2) {
1444 // Potential overflow.
1445 if (isPosInfinite(x)) {
1446 // Handle infinity
1447 return constructor.create(x, arg());
1448 }
1449 // Scale down.
1450 x /= 2;
1451 y /= 2;
1452 // log(2)
1453 re = logOf2;
1454 } else if (y < Double.MIN_NORMAL) {
1455 // Potential underflow.
1456 if (y == 0) {
1457 // Handle real only number
1458 return constructor.create(log.applyAsDouble(x), arg());
1459 }
1460 // Scale up sub-normal numbers to make them normal by scaling by 2^54,
1461 // i.e. more than the mantissa digits.
1462 x *= 0x1.0p54;
1463 y *= 0x1.0p54;
1464 // log(2^-54) = -54 * log(2)
1465 re = -54 * logOf2;
1466 }
1467 re += log.applyAsDouble(abs(x, y));
1468 }
1469
1470 // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
1471 return constructor.create(re, arg());
1472 }
1473
1474 /**
1475 * Returns the complex power of this complex number raised to the power of {@code x}.
1476 * Implements the formula:
1477 *
1478 * <p>\[ z^x = e^{x \ln(z)} \]
1479 *
1480 * <p>If this complex number is zero then this method returns zero if {@code x} is positive
1481 * in the real component and zero in the imaginary component;
1482 * otherwise it returns NaN + iNaN.
1483 *
1484 * @param x The exponent to which this complex number is to be raised.
1485 * @return This complex number raised to the power of {@code x}.
1486 * @see #log()
1487 * @see #multiply(Complex)
1488 * @see #exp()
1489 * @see <a href="http://mathworld.wolfram.com/ComplexExponentiation.html">Complex exponentiation</a>
1490 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
1491 */
1492 public Complex pow(Complex x) {
1493 if (real == 0 &&
1494 imaginary == 0) {
1495 // This value is zero. Test the other.
1496 if (x.real > 0 &&
1497 x.imaginary == 0) {
1498 // 0 raised to positive number is 0
1499 return ZERO;
1500 }
1501 // 0 raised to anything else is NaN
1502 return NAN;
1503 }
1504 return log().multiply(x).exp();
1505 }
1506
1507 /**
1508 * Returns the complex power of this complex number raised to the power of {@code x},
1509 * with {@code x} interpreted as a real number.
1510 * Implements the formula:
1511 *
1512 * <p>\[ z^x = e^{x \ln(z)} \]
1513 *
1514 * <p>If this complex number is zero then this method returns zero if {@code x} is positive;
1515 * otherwise it returns NaN + iNaN.
1516 *
1517 * @param x The exponent to which this complex number is to be raised.
1518 * @return This complex number raised to the power of {@code x}.
1519 * @see #log()
1520 * @see #multiply(double)
1521 * @see #exp()
1522 * @see #pow(Complex)
1523 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
1524 */
1525 public Complex pow(double x) {
1526 if (real == 0 &&
1527 imaginary == 0) {
1528 // This value is zero. Test the other.
1529 if (x > 0) {
1530 // 0 raised to positive number is 0
1531 return ZERO;
1532 }
1533 // 0 raised to anything else is NaN
1534 return NAN;
1535 }
1536 return log().multiply(x).exp();
1537 }
1538
1539 /**
1540 * Returns the
1541 * <a href="http://mathworld.wolfram.com/SquareRoot.html">
1542 * square root</a> of this complex number.
1543 *
1544 * <p>\[ \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) \]
1545 *
1546 * <p>The square root of \( z \) is in the range \( [0, +\infty) \) along the real axis and
1547 * is unbounded along the imaginary axis. The imaginary part of the square root has a
1548 * branch cut along the negative real axis \( (-infty,0) \). Special cases:
1549 *
1550 * <ul>
1551 * <li>{@code z.conj().sqrt() == z.sqrt().conj()}.
1552 * <li>If {@code z} is ±0 + i0, returns +0 + i0.
1553 * <li>If {@code z} is x + i∞ for all x (including NaN), returns +∞ + i∞.
1554 * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1555 * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +0 + i∞.
1556 * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
1557 * <li>If {@code z} is −∞ + iNaN, returns NaN ± i∞ (where the sign of the imaginary part of the result is unspecified).
1558 * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
1559 * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1560 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1561 * </ul>
1562 *
1563 * <p>Implements the following algorithm to compute \( \sqrt{x + iy} \):
1564 * <ol>
1565 * <li>Let \( t = \sqrt{2 (|x| + |x + iy|)} \)
1566 * <li>if \( x \geq 0 \) return \( \frac{t}{2} + i \frac{y}{t} \)
1567 * <li>else return \( \frac{|y|}{t} + i\ \text{sgn}(y) \frac{t}{2} \)
1568 * </ol>
1569 * where:
1570 * <ul>
1571 * <li>\( |x| =\ \){@link Math#abs(double) abs}(x)
1572 * <li>\( |x + y i| =\ \){@link Complex#abs}
1573 * <li>\( \text{sgn}(y) =\ \){@link Math#copySign(double,double) copySign}(1.0, y)
1574 * </ul>
1575 *
1576 * <p>The implementation is overflow and underflow safe based on the method described in:</p>
1577 * <blockquote>
1578 * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
1579 * Implementing complex elementary functions using exception handling.
1580 * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
1581 * </blockquote>
1582 *
1583 * @return The square root of this complex number.
1584 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
1585 */
1586 public Complex sqrt() {
1587 return sqrt(real, imaginary);
1588 }
1589
1590 /**
1591 * Returns the square root of the complex number {@code sqrt(x + i y)}.
1592 *
1593 * @param real Real component.
1594 * @param imaginary Imaginary component.
1595 * @return The square root of the complex number.
1596 */
1597 private static Complex sqrt(double real, double imaginary) {
1598 // Handle NaN
1599 if (Double.isNaN(real) || Double.isNaN(imaginary)) {
1600 // Check for infinite
1601 if (Double.isInfinite(imaginary)) {
1602 return new Complex(Double.POSITIVE_INFINITY, imaginary);
1603 }
1604 if (Double.isInfinite(real)) {
1605 if (real == Double.NEGATIVE_INFINITY) {
1606 return new Complex(Double.NaN, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
1607 }
1608 return new Complex(Double.POSITIVE_INFINITY, Double.NaN);
1609 }
1610 return NAN;
1611 }
1612
1613 // Compute with positive values and determine sign at the end
1614 final double x = Math.abs(real);
1615 final double y = Math.abs(imaginary);
1616
1617 // Compute
1618 double t;
1619
1620 // This alters the implementation of Hull et al (1994) which used a standard
1621 // precision representation of |z|: sqrt(x*x + y*y).
1622 // This formula should use the same definition of the magnitude returned
1623 // by Complex.abs() which is a high precision computation with scaling.
1624 // Worry about overflow if 2 * (|z| + |x|) will overflow.
1625 // Worry about underflow if |z| or |x| are sub-normal components.
1626
1627 if (inRegion(x, y, Double.MIN_NORMAL, SQRT_SAFE_UPPER)) {
1628 // No over/underflow
1629 t = Math.sqrt(2 * (abs(x, y) + x));
1630 } else {
1631 // Potential over/underflow. First check infinites and real/imaginary only.
1632
1633 // Check for infinite
1634 if (isPosInfinite(y)) {
1635 return new Complex(Double.POSITIVE_INFINITY, imaginary);
1636 } else if (isPosInfinite(x)) {
1637 if (real == Double.NEGATIVE_INFINITY) {
1638 return new Complex(0, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
1639 }
1640 return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0, imaginary));
1641 } else if (y == 0) {
1642 // Real only
1643 final double sqrtAbs = Math.sqrt(x);
1644 if (real < 0) {
1645 return new Complex(0, Math.copySign(sqrtAbs, imaginary));
1646 }
1647 return new Complex(sqrtAbs, imaginary);
1648 } else if (x == 0) {
1649 // Imaginary only. This sets the two components to the same magnitude.
1650 // Note: In polar coordinates this does not happen:
1651 // real = sqrt(abs()) * Math.cos(arg() / 2)
1652 // imag = sqrt(abs()) * Math.sin(arg() / 2)
1653 // arg() / 2 = pi/4 and cos and sin should both return sqrt(2)/2 but
1654 // are different by 1 ULP.
1655 final double sqrtAbs = Math.sqrt(y) * ONE_OVER_ROOT2;
1656 return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
1657 } else {
1658 // Over/underflow.
1659 // Full scaling is not required as this is done in the hypotenuse function.
1660 // Keep the number as big as possible for maximum precision in the second sqrt.
1661 // Note if we scale by an even power of 2, we can re-scale by sqrt of the number.
1662 // a = sqrt(b)
1663 // a = sqrt(b/4) * sqrt(4)
1664
1665 double rescale;
1666 double sx;
1667 double sy;
1668 if (Math.max(x, y) > SQRT_SAFE_UPPER) {
1669 // Overflow. Scale down by 16 and rescale by sqrt(16).
1670 sx = x / 16;
1671 sy = y / 16;
1672 rescale = 4;
1673 } else {
1674 // Sub-normal numbers. Make them normal by scaling by 2^54,
1675 // i.e. more than the mantissa digits, and rescale by sqrt(2^54) = 2^27.
1676 sx = x * 0x1.0p54;
1677 sy = y * 0x1.0p54;
1678 rescale = 0x1.0p-27;
1679 }
1680 t = rescale * Math.sqrt(2 * (abs(sx, sy) + sx));
1681 }
1682 }
1683
1684 if (real >= 0) {
1685 return new Complex(t / 2, imaginary / t);
1686 }
1687 return new Complex(y / t, Math.copySign(t / 2, imaginary));
1688 }
1689
1690 /**
1691 * Returns the
1692 * <a href="http://mathworld.wolfram.com/Sine.html">
1693 * sine</a> of this complex number.
1694 *
1695 * <p>\[ \sin(z) = \frac{1}{2} i \left( e^{-iz} - e^{iz} \right) \]
1696 *
1697 * <p>This is an odd function: \( \sin(z) = -\sin(-z) \).
1698 * The sine is an entire function and requires no branch cuts.
1699 *
1700 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1701 *
1702 * <p>\[ \sin(x + iy) = \sin(x)\cosh(y) + i \cos(x)\sinh(y) \]
1703 *
1704 * <p>As per the C99 standard this function is computed using the trigonomic identity:
1705 *
1706 * <p>\[ \sin(z) = -i \sinh(iz) \]
1707 *
1708 * @return The sine of this complex number.
1709 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
1710 */
1711 public Complex sin() {
1712 // Define in terms of sinh
1713 // sin(z) = -i sinh(iz)
1714 // Multiply this number by I, compute sinh, then multiply by back
1715 return sinh(-imaginary, real, Complex::multiplyNegativeI);
1716 }
1717
1718 /**
1719 * Returns the
1720 * <a href="http://mathworld.wolfram.com/Cosine.html">
1721 * cosine</a> of this complex number.
1722 *
1723 * <p>\[ \cos(z) = \frac{1}{2} \left( e^{iz} + e^{-iz} \right) \]
1724 *
1725 * <p>This is an even function: \( \cos(z) = \cos(-z) \).
1726 * The cosine is an entire function and requires no branch cuts.
1727 *
1728 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1729 *
1730 * <p>\[ \cos(x + iy) = \cos(x)\cosh(y) - i \sin(x)\sinh(y) \]
1731 *
1732 * <p>As per the C99 standard this function is computed using the trigonomic identity:
1733 *
1734 * <p>\[ cos(z) = cosh(iz) \]
1735 *
1736 * @return The cosine of this complex number.
1737 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
1738 */
1739 public Complex cos() {
1740 // Define in terms of cosh
1741 // cos(z) = cosh(iz)
1742 // Multiply this number by I and compute cosh.
1743 return cosh(-imaginary, real, Complex::ofCartesian);
1744 }
1745
1746 /**
1747 * Returns the
1748 * <a href="http://mathworld.wolfram.com/Tangent.html">
1749 * tangent</a> of this complex number.
1750 *
1751 * <p>\[ \tan(z) = \frac{i(e^{-iz} - e^{iz})}{e^{-iz} + e^{iz}} \]
1752 *
1753 * <p>This is an odd function: \( \tan(z) = -\tan(-z) \).
1754 * The tangent is an entire function and requires no branch cuts.
1755 *
1756 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:</p>
1757 * \[ \tan(x + iy) = \frac{\sin(2x)}{\cos(2x)+\cosh(2y)} + i \frac{\sinh(2y)}{\cos(2x)+\cosh(2y)} \]
1758 *
1759 * <p>As per the C99 standard this function is computed using the trigonomic identity:</p>
1760 * \[ \tan(z) = -i \tanh(iz) \]
1761 *
1762 * @return The tangent of this complex number.
1763 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
1764 */
1765 public Complex tan() {
1766 // Define in terms of tanh
1767 // tan(z) = -i tanh(iz)
1768 // Multiply this number by I, compute tanh, then multiply by back
1769 return tanh(-imaginary, real, Complex::multiplyNegativeI);
1770 }
1771
1772 /**
1773 * Returns the
1774 * <a href="http://mathworld.wolfram.com/InverseSine.html">
1775 * inverse sine</a> of this complex number.
1776 *
1777 * <p>\[ \sin^{-1}(z) = - i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
1778 *
1779 * <p>The inverse sine of \( z \) is unbounded along the imaginary axis and
1780 * in the range \( [-\pi, \pi] \) along the real axis. Special cases are handled
1781 * as if the operation is implemented using \( \sin^{-1}(z) = -i \sinh^{-1}(iz) \).
1782 *
1783 * <p>The inverse sine is a multivalued function and requires a branch cut in
1784 * the complex plane; the cut is conventionally placed at the line segments
1785 * \( (\infty,-1) \) and \( (1,\infty) \) of the real axis.
1786 *
1787 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1788 *
1789 * <p>\[ \begin{aligned}
1790 * \sin^{-1}(z) &= \sin^{-1}(B) + i\ \text{sgn}(y)\ln \left(A + \sqrt{A^2-1} \right) \\
1791 * A &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
1792 * B &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \]
1793 *
1794 * <p>where \( \text{sgn}(y) \) is the sign function implemented using
1795 * {@link Math#copySign(double,double) copySign(1.0, y)}.
1796 *
1797 * <p>The implementation is based on the method described in:</p>
1798 * <blockquote>
1799 * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
1800 * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
1801 * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
1802 * </blockquote>
1803 *
1804 * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
1805 * {@code c++} implementation {@code <boost/math/complex/asin.hpp>}.
1806 *
1807 * @return The inverse sine of this complex number.
1808 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
1809 */
1810 public Complex asin() {
1811 return asin(real, imaginary, Complex::ofCartesian);
1812 }
1813
1814 /**
1815 * Returns the inverse sine of the complex number.
1816 *
1817 * <p>This function exists to allow implementation of the identity
1818 * {@code asinh(z) = -i asin(iz)}.
1819 *
1820 * <p>Adapted from {@code <boost/math/complex/asin.hpp>}. This method only (and not
1821 * invoked methods within) is distributed under the Boost Software License V1.0.
1822 * The original notice is shown below and the licence is shown in full in LICENSE:
1823 * <pre>
1824 * (C) Copyright John Maddock 2005.
1825 * Distributed under the Boost Software License, Version 1.0. (See accompanying
1826 * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
1827 * </pre>
1828 *
1829 * @param real Real part.
1830 * @param imaginary Imaginary part.
1831 * @param constructor Constructor.
1832 * @return The inverse sine of this complex number.
1833 */
1834 private static Complex asin(final double real, final double imaginary,
1835 final ComplexConstructor constructor) {
1836 // Compute with positive values and determine sign at the end
1837 final double x = Math.abs(real);
1838 final double y = Math.abs(imaginary);
1839 // The result (without sign correction)
1840 double re;
1841 double im;
1842
1843 // Handle C99 special cases
1844 if (Double.isNaN(x)) {
1845 if (isPosInfinite(y)) {
1846 re = x;
1847 im = y;
1848 } else {
1849 // No-use of the input constructor
1850 return NAN;
1851 }
1852 } else if (Double.isNaN(y)) {
1853 if (x == 0) {
1854 re = 0;
1855 im = y;
1856 } else if (isPosInfinite(x)) {
1857 re = y;
1858 im = x;
1859 } else {
1860 // No-use of the input constructor
1861 return NAN;
1862 }
1863 } else if (isPosInfinite(x)) {
1864 re = isPosInfinite(y) ? PI_OVER_4 : PI_OVER_2;
1865 im = x;
1866 } else if (isPosInfinite(y)) {
1867 re = 0;
1868 im = y;
1869 } else {
1870 // Special case for real numbers:
1871 if (y == 0 && x <= 1) {
1872 return constructor.create(Math.asin(real), imaginary);
1873 }
1874
1875 final double xp1 = x + 1;
1876 final double xm1 = x - 1;
1877
1878 if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
1879 final double yy = y * y;
1880 final double r = Math.sqrt(xp1 * xp1 + yy);
1881 final double s = Math.sqrt(xm1 * xm1 + yy);
1882 final double a = 0.5 * (r + s);
1883 final double b = x / a;
1884
1885 if (b <= B_CROSSOVER) {
1886 re = Math.asin(b);
1887 } else {
1888 final double apx = a + x;
1889 if (x <= 1) {
1890 re = Math.atan(x / Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))));
1891 } else {
1892 re = Math.atan(x / (y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))));
1893 }
1894 }
1895
1896 if (a <= A_CROSSOVER) {
1897 double am1;
1898 if (x < 1) {
1899 am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
1900 } else {
1901 am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
1902 }
1903 im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
1904 } else {
1905 im = Math.log(a + Math.sqrt(a * a - 1));
1906 }
1907 } else {
1908 // Hull et al: Exception handling code from figure 4
1909 if (y <= (EPSILON * Math.abs(xm1))) {
1910 if (x < 1) {
1911 re = Math.asin(x);
1912 im = y / Math.sqrt(xp1 * (1 - x));
1913 } else {
1914 re = PI_OVER_2;
1915 if ((Double.MAX_VALUE / xp1) > xm1) {
1916 // xp1 * xm1 won't overflow:
1917 im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
1918 } else {
1919 im = LN_2 + Math.log(x);
1920 }
1921 }
1922 } else if (y <= SAFE_MIN) {
1923 // Hull et al: Assume x == 1.
1924 // True if:
1925 // E^2 > 8*sqrt(u)
1926 //
1927 // E = Machine epsilon: (1 + epsilon) = 1
1928 // u = Double.MIN_NORMAL
1929 re = PI_OVER_2 - Math.sqrt(y);
1930 im = Math.sqrt(y);
1931 } else if (EPSILON * y - 1 >= x) {
1932 // Possible underflow:
1933 re = x / y;
1934 im = LN_2 + Math.log(y);
1935 } else if (x > 1) {
1936 re = Math.atan(x / y);
1937 final double xoy = x / y;
1938 im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
1939 } else {
1940 final double a = Math.sqrt(1 + y * y);
1941 // Possible underflow:
1942 re = x / a;
1943 im = 0.5 * Math.log1p(2 * y * (y + a));
1944 }
1945 }
1946 }
1947
1948 return constructor.create(changeSign(re, real),
1949 changeSign(im, imaginary));
1950 }
1951
1952 /**
1953 * Returns the
1954 * <a href="http://mathworld.wolfram.com/InverseCosine.html">
1955 * inverse cosine</a> of this complex number.
1956 *
1957 * <p>\[ \cos^{-1}(z) = \frac{\pi}{2} + i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
1958 *
1959 * <p>The inverse cosine of \( z \) is in the range \( [0, \pi) \) along the real axis and
1960 * unbounded along the imaginary axis. Special cases:
1961 *
1962 * <ul>
1963 * <li>{@code z.conj().acos() == z.acos().conj()}.
1964 * <li>If {@code z} is ±0 + i0, returns π/2 − i0.
1965 * <li>If {@code z} is ±0 + iNaN, returns π/2 + iNaN.
1966 * <li>If {@code z} is x + i∞ for finite x, returns π/2 − i∞.
1967 * <li>If {@code z} is x + iNaN, returns NaN + iNaN ("invalid" floating-point operation).
1968 * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns π − i∞.
1969 * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +0 − i∞.
1970 * <li>If {@code z} is −∞ + i∞, returns 3π/4 − i∞.
1971 * <li>If {@code z} is +∞ + i∞, returns π/4 − i∞.
1972 * <li>If {@code z} is ±∞ + iNaN, returns NaN ± i∞ where the sign of the imaginary part of the result is unspecified.
1973 * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1974 * <li>If {@code z} is NaN + i∞, returns NaN − i∞.
1975 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1976 * </ul>
1977 *
1978 * <p>The inverse cosine is a multivalued function and requires a branch cut in
1979 * the complex plane; the cut is conventionally placed at the line segments
1980 * \( (-\infty,-1) \) and \( (1,\infty) \) of the real axis.
1981 *
1982 * <p>This function is implemented using real \( x \) and imaginary \( y \) parts:
1983 *
1984 * <p>\[ \begin{aligned}
1985 * \cos^{-1}(z) &= \cos^{-1}(B) - i\ \text{sgn}(y) \ln\left(A + \sqrt{A^2-1}\right) \\
1986 * A &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
1987 * B &= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \]
1988 *
1989 * <p>where \( \text{sgn}(y) \) is the sign function implemented using
1990 * {@link Math#copySign(double,double) copySign(1.0, y)}.
1991 *
1992 * <p>The implementation is based on the method described in:</p>
1993 * <blockquote>
1994 * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
1995 * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
1996 * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
1997 * </blockquote>
1998 *
1999 * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
2000 * {@code c++} implementation {@code <boost/math/complex/acos.hpp>}.
2001 *
2002 * @return The inverse cosine of this complex number.
2003 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
2004 */
2005 public Complex acos() {
2006 return acos(real, imaginary, Complex::ofCartesian);
2007 }
2008
2009 /**
2010 * Returns the inverse cosine of the complex number.
2011 *
2012 * <p>This function exists to allow implementation of the identity
2013 * {@code acosh(z) = +-i acos(z)}.
2014 *
2015 * <p>Adapted from {@code <boost/math/complex/acos.hpp>}. This method only (and not
2016 * invoked methods within) is distributed under the Boost Software License V1.0.
2017 * The original notice is shown below and the licence is shown in full in LICENSE:
2018 * <pre>
2019 * (C) Copyright John Maddock 2005.
2020 * Distributed under the Boost Software License, Version 1.0. (See accompanying
2021 * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
2022 * </pre>
2023 *
2024 * @param real Real part.
2025 * @param imaginary Imaginary part.
2026 * @param constructor Constructor.
2027 * @return The inverse cosine of the complex number.
2028 */
2029 private static Complex acos(final double real, final double imaginary,
2030 final ComplexConstructor constructor) {
2031 // Compute with positive values and determine sign at the end
2032 final double x = Math.abs(real);
2033 final double y = Math.abs(imaginary);
2034 // The result (without sign correction)
2035 double re;
2036 double im;
2037
2038 // Handle C99 special cases
2039 if (isPosInfinite(x)) {
2040 if (isPosInfinite(y)) {
2041 re = PI_OVER_4;
2042 im = y;
2043 } else if (Double.isNaN(y)) {
2044 // sign of the imaginary part of the result is unspecified
2045 return constructor.create(imaginary, real);
2046 } else {
2047 re = 0;
2048 im = Double.POSITIVE_INFINITY;
2049 }
2050 } else if (Double.isNaN(x)) {
2051 if (isPosInfinite(y)) {
2052 return constructor.create(x, -imaginary);
2053 }
2054 // No-use of the input constructor
2055 return NAN;
2056 } else if (isPosInfinite(y)) {
2057 re = PI_OVER_2;
2058 im = y;
2059 } else if (Double.isNaN(y)) {
2060 return constructor.create(x == 0 ? PI_OVER_2 : y, y);
2061 } else {
2062 // Special case for real numbers:
2063 if (y == 0 && x <= 1) {
2064 return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
2065 }
2066
2067 final double xp1 = x + 1;
2068 final double xm1 = x - 1;
2069
2070 if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
2071 final double yy = y * y;
2072 final double r = Math.sqrt(xp1 * xp1 + yy);
2073 final double s = Math.sqrt(xm1 * xm1 + yy);
2074 final double a = 0.5 * (r + s);
2075 final double b = x / a;
2076
2077 if (b <= B_CROSSOVER) {
2078 re = Math.acos(b);
2079 } else {
2080 final double apx = a + x;
2081 if (x <= 1) {
2082 re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
2083 } else {
2084 re = Math.atan((y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))) / x);
2085 }
2086 }
2087
2088 if (a <= A_CROSSOVER) {
2089 double am1;
2090 if (x < 1) {
2091 am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
2092 } else {
2093 am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
2094 }
2095 im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
2096 } else {
2097 im = Math.log(a + Math.sqrt(a * a - 1));
2098 }
2099 } else {
2100 // Hull et al: Exception handling code from figure 6
2101 if (y <= (EPSILON * Math.abs(xm1))) {
2102 if (x < 1) {
2103 re = Math.acos(x);
2104 im = y / Math.sqrt(xp1 * (1 - x));
2105 } else {
2106 // This deviates from Hull et al's paper as per
2107 // https://svn.boost.org/trac/boost/ticket/7290
2108 if ((Double.MAX_VALUE / xp1) > xm1) {
2109 // xp1 * xm1 won't overflow:
2110 re = y / Math.sqrt(xm1 * xp1);
2111 im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
2112 } else {
2113 re = y / x;
2114 im = LN_2 + Math.log(x);
2115 }
2116 }
2117 } else if (y <= SAFE_MIN) {
2118 // Hull et al: Assume x == 1.
2119 // True if:
2120 // E^2 > 8*sqrt(u)
2121 //
2122 // E = Machine epsilon: (1 + epsilon) = 1
2123 // u = Double.MIN_NORMAL
2124 re = Math.sqrt(y);
2125 im = Math.sqrt(y);
2126 } else if (EPSILON * y - 1 >= x) {
2127 re = PI_OVER_2;
2128 im = LN_2 + Math.log(y);
2129 } else if (x > 1) {
2130 re = Math.atan(y / x);
2131 final double xoy = x / y;
2132 im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
2133 } else {
2134 re = PI_OVER_2;
2135 final double a = Math.sqrt(1 + y * y);
2136 im = 0.5 * Math.log1p(2 * y * (y + a));
2137 }
2138 }
2139 }
2140
2141 return constructor.create(negative(real) ? Math.PI - re : re,
2142 negative(imaginary) ? im : -im);
2143 }
2144
2145 /**
2146 * Returns the
2147 * <a href="http://mathworld.wolfram.com/InverseTangent.html">
2148 * inverse tangent</a> of this complex number.
2149 *
2150 * <p>\[ \tan^{-1}(z) = \frac{i}{2} \ln \left( \frac{i + z}{i - z} \right) \]
2151 *
2152 * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the imaginary axis and
2153 * in the range \( [-\pi/2, \pi/2] \) along the real axis.
2154 *
2155 * <p>The inverse tangent is a multivalued function and requires a branch cut in
2156 * the complex plane; the cut is conventionally placed at the line segments
2157 * \( (i \infty,-i] \) and \( [i,i \infty) \) of the imaginary axis.
2158 *
2159 * <p>As per the C99 standard this function is computed using the trigonomic identity:
2160 * \[ \tan^{-1}(z) = -i \tanh^{-1}(iz) \]
2161 *
2162 * @return The inverse tangent of this complex number.
2163 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
2164 */
2165 public Complex atan() {
2166 // Define in terms of atanh
2167 // atan(z) = -i atanh(iz)
2168 // Multiply this number by I, compute atanh, then multiply by back
2169 return atanh(-imaginary, real, Complex::multiplyNegativeI);
2170 }
2171
2172 /**
2173 * Returns the
2174 * <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
2175 * hyperbolic sine</a> of this complex number.
2176 *
2177 * <p>\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \]
2178 *
2179 * <p>The hyperbolic sine of \( z \) is an entire function in the complex plane
2180 * and is periodic with respect to the imaginary component with period \( 2\pi i \).
2181 * Special cases:
2182 *
2183 * <ul>
2184 * <li>{@code z.conj().sinh() == z.sinh().conj()}.
2185 * <li>This is an odd function: \( \sinh(z) = -\sinh(-z) \).
2186 * <li>If {@code z} is +0 + i0, returns +0 + i0.
2187 * <li>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).
2188 * <li>If {@code z} is +0 + iNaN, returns ±0 + iNaN (where the sign of the real part of the result is unspecified).
2189 * <li>If {@code z} is x + i∞ for positive finite x, returns NaN + iNaN ("invalid" floating-point operation).
2190 * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2191 * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
2192 * <li>If {@code z} is +∞ + iy for positive finite y, returns +∞ cis(y) (see {@link #ofCis(double)}.
2193 * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
2194 * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2195 * <li>If {@code z} is NaN + i0, returns NaN + i0.
2196 * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2197 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2198 * </ul>
2199 *
2200 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2201 *
2202 * <p>\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \]
2203 *
2204 * @return The hyperbolic sine of this complex number.
2205 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
2206 */
2207 public Complex sinh() {
2208 return sinh(real, imaginary, Complex::ofCartesian);
2209 }
2210
2211 /**
2212 * Returns the hyperbolic sine of the complex number.
2213 *
2214 * <p>This function exists to allow implementation of the identity
2215 * {@code sin(z) = -i sinh(iz)}.<p>
2216 *
2217 * @param real Real part.
2218 * @param imaginary Imaginary part.
2219 * @param constructor Constructor.
2220 * @return The hyperbolic sine of the complex number.
2221 */
2222 private static Complex sinh(double real, double imaginary, ComplexConstructor constructor) {
2223 if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
2224 return constructor.create(real, Double.NaN);
2225 }
2226 if (real == 0) {
2227 // Imaginary-only sinh(iy) = i sin(y).
2228 if (Double.isFinite(imaginary)) {
2229 // Maintain periodic property with respect to the imaginary component.
2230 // sinh(+/-0.0) * cos(+/-x) = +/-0 * cos(x)
2231 return constructor.create(changeSign(real, Math.cos(imaginary)),
2232 Math.sin(imaginary));
2233 }
2234 // If imaginary is inf/NaN the sign of the real part is unspecified.
2235 // Returning the same real value maintains the conjugate equality.
2236 // It is not possible to also maintain the odd function (hence the unspecified sign).
2237 return constructor.create(real, Double.NaN);
2238 }
2239 if (imaginary == 0) {
2240 // Real-only sinh(x).
2241 return constructor.create(Math.sinh(real), imaginary);
2242 }
2243 final double x = Math.abs(real);
2244 if (x > SAFE_EXP) {
2245 // Approximate sinh/cosh(x) using exp^|x| / 2
2246 return coshsinh(x, real, imaginary, true, constructor);
2247 }
2248 // No overflow of sinh/cosh
2249 return constructor.create(Math.sinh(real) * Math.cos(imaginary),
2250 Math.cosh(real) * Math.sin(imaginary));
2251 }
2252
2253 /**
2254 * Returns the
2255 * <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
2256 * hyperbolic cosine</a> of this complex number.
2257 *
2258 * <p>\[ \cosh(z) = \frac{1}{2} \left( e^{z} + e^{-z} \right) \]
2259 *
2260 * <p>The hyperbolic cosine of \( z \) is an entire function in the complex plane
2261 * and is periodic with respect to the imaginary component with period \( 2\pi i \).
2262 * Special cases:
2263 *
2264 * <ul>
2265 * <li>{@code z.conj().cosh() == z.cosh().conj()}.
2266 * <li>This is an even function: \( \cosh(z) = \cosh(-z) \).
2267 * <li>If {@code z} is +0 + i0, returns 1 + i0.
2268 * <li>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).
2269 * <li>If {@code z} is +0 + iNaN, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
2270 * <li>If {@code z} is x + i∞ for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2271 * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2272 * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
2273 * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y) (see {@link #ofCis(double)}).
2274 * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2275 * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
2276 * <li>If {@code z} is NaN + i0, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
2277 * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2278 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2279 * </ul>
2280 *
2281 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2282 *
2283 * <p>\[ \cosh(x + iy) = \cosh(x)\cos(y) + i \sinh(x)\sin(y) \]
2284 *
2285 * @return The hyperbolic cosine of this complex number.
2286 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cosh/">Cosh</a>
2287 */
2288 public Complex cosh() {
2289 return cosh(real, imaginary, Complex::ofCartesian);
2290 }
2291
2292 /**
2293 * Returns the hyperbolic cosine of the complex number.
2294 *
2295 * <p>This function exists to allow implementation of the identity
2296 * {@code cos(z) = cosh(iz)}.<p>
2297 *
2298 * @param real Real part.
2299 * @param imaginary Imaginary part.
2300 * @param constructor Constructor.
2301 * @return The hyperbolic cosine of the complex number.
2302 */
2303 private static Complex cosh(double real, double imaginary, ComplexConstructor constructor) {
2304 // ISO C99: Preserve the even function by mapping to positive
2305 // f(z) = f(-z)
2306 if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
2307 return constructor.create(Math.abs(real), Double.NaN);
2308 }
2309 if (real == 0) {
2310 // Imaginary-only cosh(iy) = cos(y).
2311 if (Double.isFinite(imaginary)) {
2312 // Maintain periodic property with respect to the imaginary component.
2313 // sinh(+/-0.0) * sin(+/-x) = +/-0 * sin(x)
2314 return constructor.create(Math.cos(imaginary),
2315 changeSign(real, Math.sin(imaginary)));
2316 }
2317 // If imaginary is inf/NaN the sign of the imaginary part is unspecified.
2318 // Although not required by C99 changing the sign maintains the conjugate equality.
2319 // It is not possible to also maintain the even function (hence the unspecified sign).
2320 return constructor.create(Double.NaN, changeSign(real, imaginary));
2321 }
2322 if (imaginary == 0) {
2323 // Real-only cosh(x).
2324 // Change sign to preserve conjugate equality and even function.
2325 // sin(+/-0) * sinh(+/-x) = +/-0 * +/-a (sinh is monotonic and same sign)
2326 // => change the sign of imaginary using real. Handles special case of infinite real.
2327 // If real is NaN the sign of the imaginary part is unspecified.
2328 return constructor.create(Math.cosh(real), changeSign(imaginary, real));
2329 }
2330 final double x = Math.abs(real);
2331 if (x > SAFE_EXP) {
2332 // Approximate sinh/cosh(x) using exp^|x| / 2
2333 return coshsinh(x, real, imaginary, false, constructor);
2334 }
2335 // No overflow of sinh/cosh
2336 return constructor.create(Math.cosh(real) * Math.cos(imaginary),
2337 Math.sinh(real) * Math.sin(imaginary));
2338 }
2339
2340 /**
2341 * Compute cosh or sinh when the absolute real component |x| is large. In this case
2342 * cosh(x) and sinh(x) can be approximated by exp(|x|) / 2:
2343 *
2344 * <pre>
2345 * cosh(x+iy) real = (e^|x| / 2) * cos(y)
2346 * cosh(x+iy) imag = (e^|x| / 2) * sin(y) * sign(x)
2347 * sinh(x+iy) real = (e^|x| / 2) * cos(y) * sign(x)
2348 * sinh(x+iy) imag = (e^|x| / 2) * sin(y)
2349 * </pre>
2350 *
2351 * @param x Absolute real component |x|.
2352 * @param real Real part (x).
2353 * @param imaginary Imaginary part (y).
2354 * @param sinh Set to true to compute sinh, otherwise cosh.
2355 * @param constructor Constructor.
2356 * @return The hyperbolic sine/cosine of the complex number.
2357 */
2358 private static Complex coshsinh(double x, double real, double imaginary, boolean sinh,
2359 ComplexConstructor constructor) {
2360 // Always require the cos and sin.
2361 double re = Math.cos(imaginary);
2362 double im = Math.sin(imaginary);
2363 // Compute the correct function
2364 if (sinh) {
2365 re = changeSign(re, real);
2366 } else {
2367 im = changeSign(im, real);
2368 }
2369 // Multiply by (e^|x| / 2).
2370 // Overflow safe computation since sin/cos can be very small allowing a result
2371 // when e^x overflows: e^x / 2 = (e^m / 2) * e^m * e^(x-2m)
2372 if (x > SAFE_EXP * 3) {
2373 // e^x > e^m * e^m * e^m
2374 // y * (e^m / 2) * e^m * e^m will overflow when starting with Double.MIN_VALUE.
2375 // Note: Do not multiply by +inf to safeguard against sin(y)=0.0 which
2376 // will create 0 * inf = nan.
2377 re *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
2378 im *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
2379 } else {
2380 // Initial part of (e^x / 2) using (e^m / 2)
2381 re *= EXP_M / 2;
2382 im *= EXP_M / 2;
2383 double xm;
2384 if (x > SAFE_EXP * 2) {
2385 // e^x = e^m * e^m * e^(x-2m)
2386 re *= EXP_M;
2387 im *= EXP_M;
2388 xm = x - SAFE_EXP * 2;
2389 } else {
2390 // e^x = e^m * e^(x-m)
2391 xm = x - SAFE_EXP;
2392 }
2393 final double exp = Math.exp(xm);
2394 re *= exp;
2395 im *= exp;
2396 }
2397 return constructor.create(re, im);
2398 }
2399
2400 /**
2401 * Returns the
2402 * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
2403 * hyperbolic tangent</a> of this complex number.
2404 *
2405 * <p>\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \]
2406 *
2407 * <p>The hyperbolic tangent of \( z \) is an entire function in the complex plane
2408 * and is periodic with respect to the imaginary component with period \( \pi i \)
2409 * and has poles of the first order along the imaginary line, at coordinates
2410 * \( (0, \pi(\frac{1}{2} + n)) \).
2411 * Note that the {@code double} floating-point representation is unable to exactly represent
2412 * \( \pi/2 \) and there is no value for which a pole error occurs. Special cases:
2413 *
2414 * <ul>
2415 * <li>{@code z.conj().tanh() == z.tanh().conj()}.
2416 * <li>This is an odd function: \( \tanh(z) = -\tanh(-z) \).
2417 * <li>If {@code z} is +0 + i0, returns +0 + i0.
2418 * <li>If {@code z} is 0 + i∞, returns 0 + iNaN.
2419 * <li>If {@code z} is x + i∞ for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2420 * <li>If {@code z} is 0 + iNaN, returns 0 + iNAN.
2421 * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2422 * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns 1 + i0 sin(2y).
2423 * <li>If {@code z} is +∞ + i∞, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
2424 * <li>If {@code z} is +∞ + iNaN, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
2425 * <li>If {@code z} is NaN + i0, returns NaN + i0.
2426 * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2427 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2428 * </ul>
2429 *
2430 * <p>Special cases include the technical corrigendum
2431 * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
2432 * DR 471: Complex math functions cacosh and ctanh</a>.
2433 *
2434 * <p>This is defined using real \( x \) and imaginary \( y \) parts:
2435 *
2436 * <p>\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \]
2437 *
2438 * <p>The implementation uses double-angle identities to avoid overflow of {@code 2x}
2439 * and {@code 2y}.
2440 *
2441 * @return The hyperbolic tangent of this complex number.
2442 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
2443 */
2444 public Complex tanh() {
2445 return tanh(real, imaginary, Complex::ofCartesian);
2446 }
2447
2448 /**
2449 * Returns the hyperbolic tangent of this complex number.
2450 *
2451 * <p>This function exists to allow implementation of the identity
2452 * {@code tan(z) = -i tanh(iz)}.<p>
2453 *
2454 * @param real Real part.
2455 * @param imaginary Imaginary part.
2456 * @param constructor Constructor.
2457 * @return The hyperbolic tangent of the complex number.
2458 */
2459 private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
2460 // Cache the absolute real value
2461 final double x = Math.abs(real);
2462
2463 // Handle inf or nan.
2464 if (!isPosFinite(x) || !Double.isFinite(imaginary)) {
2465 if (isPosInfinite(x)) {
2466 if (Double.isFinite(imaginary)) {
2467 // The sign is copied from sin(2y)
2468 // The identity sin(2a) = 2 sin(a) cos(a) is used for consistency
2469 // with the computation below. Only the magnitude is important
2470 // so drop the 2. When |y| is small sign(sin(2y)) = sign(y).
2471 final double sign = Math.abs(imaginary) < PI_OVER_2 ?
2472 imaginary :
2473 Math.sin(imaginary) * Math.cos(imaginary);
2474 return constructor.create(Math.copySign(1, real),
2475 Math.copySign(0, sign));
2476 }
2477 // imaginary is infinite or NaN
2478 return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
2479 }
2480 // Remaining cases:
2481 // (0 + i inf), returns (0 + i NaN)
2482 // (0 + i NaN), returns (0 + i NaN)
2483 // (x + i inf), returns (NaN + i NaN) for non-zero x (including infinite)
2484 // (x + i NaN), returns (NaN + i NaN) for non-zero x (including infinite)
2485 // (NaN + i 0), returns (NaN + i 0)
2486 // (NaN + i y), returns (NaN + i NaN) for non-zero y (including infinite)
2487 // (NaN + i NaN), returns (NaN + i NaN)
2488 return constructor.create(real == 0 ? real : Double.NaN,
2489 imaginary == 0 ? imaginary : Double.NaN);
2490 }
2491
2492 // Finite components
2493 // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))
2494
2495 if (real == 0) {
2496 // Imaginary-only tanh(iy) = i tan(y)
2497 // Identity: sin 2y / (1 + cos 2y) = tan(y)
2498 return constructor.create(real, Math.tan(imaginary));
2499 }
2500 if (imaginary == 0) {
2501 // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
2502 return constructor.create(Math.tanh(real), imaginary);
2503 }
2504
2505 // The double angles can be avoided using the identities:
2506 // sinh(2x) = 2 sinh(x) cosh(x)
2507 // sin(2y) = 2 sin(y) cos(y)
2508 // cosh(2x) = 2 sinh^2(x) + 1
2509 // cos(2y) = 2 cos^2(y) - 1
2510 // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
2511 // To avoid a junction when swapping between the double angles and the identities
2512 // the identities are used in all cases.
2513
2514 if (x > SAFE_EXP / 2) {
2515 // Potential overflow in sinh/cosh(2x).
2516 // Approximate sinh/cosh using exp^x.
2517 // Ignore cos^2(y) in the divisor as it is insignificant.
2518 // real = sinh(x)cosh(x) / sinh^2(x) = +/-1
2519 final double re = Math.copySign(1, real);
2520 // imag = sin(2y) / 2 sinh^2(x)
2521 // sinh(x) -> sign(x) * e^|x| / 2 when x is large.
2522 // sinh^2(x) -> e^2|x| / 4 when x is large.
2523 // imag = sin(2y) / 2 (e^2|x| / 4) = 2 sin(2y) / e^2|x|
2524 // = 4 * sin(y) cos(y) / e^2|x|
2525 // Underflow safe divide as e^2|x| may overflow:
2526 // imag = 4 * sin(y) cos(y) / e^m / e^(2|x| - m)
2527 // (|im| is a maximum of 2)
2528 double im = Math.sin(imaginary) * Math.cos(imaginary);
2529 if (x > SAFE_EXP) {
2530 // e^2|x| > e^m * e^m
2531 // This will underflow 2.0 / e^m / e^m
2532 im = Math.copySign(0.0, im);
2533 } else {
2534 // e^2|x| = e^m * e^(2|x| - m)
2535 im = 4 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
2536 }
2537 return constructor.create(re, im);
2538 }
2539
2540 // No overflow of sinh(2x) and cosh(2x)
2541
2542 // Note: This does not use the definitional formula but uses the identity:
2543 // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
2544
2545 final double sinhx = Math.sinh(real);
2546 final double coshx = Math.cosh(real);
2547 final double siny = Math.sin(imaginary);
2548 final double cosy = Math.cos(imaginary);
2549 final double divisor = sinhx * sinhx + cosy * cosy;
2550 return constructor.create(sinhx * coshx / divisor,
2551 siny * cosy / divisor);
2552 }
2553
2554 /**
2555 * Returns the
2556 * <a href="http://mathworld.wolfram.com/InverseHyperbolicSine.html">
2557 * inverse hyperbolic sine</a> of this complex number.
2558 *
2559 * <p>\[ \sinh^{-1}(z) = \ln \left(z + \sqrt{1 + z^2} \right) \]
2560 *
2561 * <p>The inverse hyperbolic sine of \( z \) is unbounded along the real axis and
2562 * in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
2563 *
2564 * <ul>
2565 * <li>{@code z.conj().asinh() == z.asinh().conj()}.
2566 * <li>This is an odd function: \( \sinh^{-1}(z) = -\sinh^{-1}(-z) \).
2567 * <li>If {@code z} is +0 + i0, returns 0 + i0.
2568 * <li>If {@code z} is x + i∞ for positive-signed finite x, returns +∞ + iπ/2.
2569 * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
2570 * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
2571 * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
2572 * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
2573 * <li>If {@code z} is NaN + i0, returns NaN + i0.
2574 * <li>If {@code z} is NaN + iy for finite nonzero y, returns NaN + iNaN ("invalid" floating-point operation).
2575 * <li>If {@code z} is NaN + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2576 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2577 * </ul>
2578 *
2579 * <p>The inverse hyperbolic sine is a multivalued function and requires a branch cut in
2580 * the complex plane; the cut is conventionally placed at the line segments
2581 * \( (-i \infty,-i) \) and \( (i,i \infty) \) of the imaginary axis.
2582 *
2583 * <p>This function is computed using the trigonomic identity:
2584 *
2585 * <p>\[ \sinh^{-1}(z) = -i \sin^{-1}(iz) \]
2586 *
2587 * @return The inverse hyperbolic sine of this complex number.
2588 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
2589 */
2590 public Complex asinh() {
2591 // Define in terms of asin
2592 // asinh(z) = -i asin(iz)
2593 // Note: This is the opposite to the identity defined in the C99 standard:
2594 // asin(z) = -i asinh(iz)
2595 // Multiply this number by I, compute asin, then multiply by back
2596 return asin(-imaginary, real, Complex::multiplyNegativeI);
2597 }
2598
2599 /**
2600 * Returns the
2601 * <a href="http://mathworld.wolfram.com/InverseHyperbolicCosine.html">
2602 * inverse hyperbolic cosine</a> of this complex number.
2603 *
2604 * <p>\[ \cosh^{-1}(z) = \ln \left(z + \sqrt{z + 1} \sqrt{z - 1} \right) \]
2605 *
2606 * <p>The inverse hyperbolic cosine of \( z \) is in the range \( [0, \infty) \) along the
2607 * real axis and in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
2608 *
2609 * <ul>
2610 * <li>{@code z.conj().acosh() == z.acosh().conj()}.
2611 * <li>If {@code z} is ±0 + i0, returns +0 + iπ/2.
2612 * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
2613 * <li>If {@code z} is 0 + iNaN, returns NaN + iπ/2 <sup>[1]</sup>.
2614 * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2615 * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns +∞ + iπ.
2616 * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
2617 * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
2618 * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
2619 * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
2620 * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
2621 * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
2622 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2623 * </ul>
2624 *
2625 * <p>Special cases include the technical corrigendum
2626 * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
2627 * DR 471: Complex math functions cacosh and ctanh</a>.
2628 *
2629 * <p>The inverse hyperbolic cosine is a multivalued function and requires a branch cut in
2630 * the complex plane; the cut is conventionally placed at the line segment
2631 * \( (-\infty,-1) \) of the real axis.
2632 *
2633 * <p>This function is computed using the trigonomic identity:
2634 *
2635 * <p>\[ \cosh^{-1}(z) = \pm i \cos^{-1}(z) \]
2636 *
2637 * <p>The sign of the multiplier is chosen to give {@code z.acosh().real() >= 0}
2638 * and compatibility with the C99 standard.
2639 *
2640 * @return The inverse hyperbolic cosine of this complex number.
2641 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
2642 */
2643 public Complex acosh() {
2644 // Define in terms of acos
2645 // acosh(z) = +-i acos(z)
2646 // Note the special case:
2647 // acos(+-0 + iNaN) = π/2 + iNaN
2648 // acosh(0 + iNaN) = NaN + iπ/2
2649 // will not appropriately multiply by I to maintain positive imaginary if
2650 // acos() imaginary computes as NaN. So do this explicitly.
2651 if (Double.isNaN(imaginary) && real == 0) {
2652 return new Complex(Double.NaN, PI_OVER_2);
2653 }
2654 return acos(real, imaginary, (re, im) ->
2655 // Set the sign appropriately for real >= 0
2656 (negative(im)) ?
2657 // Multiply by I
2658 new Complex(-im, re) :
2659 // Multiply by -I
2660 new Complex(im, -re)
2661 );
2662 }
2663
2664 /**
2665 * Returns the
2666 * <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
2667 * inverse hyperbolic tangent</a> of this complex number.
2668 *
2669 * <p>\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \]
2670 *
2671 * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the real axis and
2672 * in the range \( [-\pi/2, \pi/2] \) along the imaginary axis. Special cases:
2673 *
2674 * <ul>
2675 * <li>{@code z.conj().atanh() == z.atanh().conj()}.
2676 * <li>This is an odd function: \( \tanh^{-1}(z) = -\tanh^{-1}(-z) \).
2677 * <li>If {@code z} is +0 + i0, returns +0 + i0.
2678 * <li>If {@code z} is +0 + iNaN, returns +0 + iNaN.
2679 * <li>If {@code z} is +1 + i0, returns +∞ + i0 ("divide-by-zero" floating-point operation).
2680 * <li>If {@code z} is x + i∞ for finite positive-signed x, returns +0 + iπ/2.
2681 * <li>If {@code z} is x+iNaN for nonzero finite x, returns NaN+iNaN ("invalid" floating-point operation).
2682 * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +0 + iπ/2.
2683 * <li>If {@code z} is +∞ + i∞, returns +0 + iπ/2.
2684 * <li>If {@code z} is +∞ + iNaN, returns +0 + iNaN.
2685 * <li>If {@code z} is NaN+iy for finite y, returns NaN+iNaN ("invalid" floating-point operation).
2686 * <li>If {@code z} is NaN + i∞, returns ±0 + iπ/2 (where the sign of the real part of the result is unspecified).
2687 * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2688 * </ul>
2689 *
2690 * <p>The inverse hyperbolic tangent is a multivalued function and requires a branch cut in
2691 * the complex plane; the cut is conventionally placed at the line segments
2692 * \( (\infty,-1] \) and \( [1,\infty) \) of the real axis.
2693 *
2694 * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2695 *
2696 * <p>\[ \tanh^{-1}(z) = \frac{1}{4} \ln \left(1 + \frac{4x}{(1-x)^2+y^2} \right) + \\
2697 * 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) \]
2698 *
2699 * <p>The imaginary part is computed using {@link Math#atan2(double, double)} to ensure the
2700 * correct quadrant is returned from \( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) \).
2701 *
2702 * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
2703 * {@code c++} implementation {@code <boost/math/complex/atanh.hpp>}.
2704 *
2705 * @return The inverse hyperbolic tangent of this complex number.
2706 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
2707 */
2708 public Complex atanh() {
2709 return atanh(real, imaginary, Complex::ofCartesian);
2710 }
2711
2712 /**
2713 * Returns the inverse hyperbolic tangent of this complex number.
2714 *
2715 * <p>This function exists to allow implementation of the identity
2716 * {@code atan(z) = -i atanh(iz)}.
2717 *
2718 * <p>Adapted from {@code <boost/math/complex/atanh.hpp>}. This method only (and not
2719 * invoked methods within) is distributed under the Boost Software License V1.0.
2720 * The original notice is shown below and the licence is shown in full in LICENSE:
2721 * <pre>
2722 * (C) Copyright John Maddock 2005.
2723 * Distributed under the Boost Software License, Version 1.0. (See accompanying
2724 * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
2725 * </pre>
2726 *
2727 * @param real Real part.
2728 * @param imaginary Imaginary part.
2729 * @param constructor Constructor.
2730 * @return The inverse hyperbolic tangent of the complex number.
2731 */
2732 private static Complex atanh(final double real, final double imaginary,
2733 final ComplexConstructor constructor) {
2734 // Compute with positive values and determine sign at the end
2735 double x = Math.abs(real);
2736 double y = Math.abs(imaginary);
2737 // The result (without sign correction)
2738 double re;
2739 double im;
2740
2741 // Handle C99 special cases
2742 if (Double.isNaN(x)) {
2743 if (isPosInfinite(y)) {
2744 // The sign of the real part of the result is unspecified
2745 return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
2746 }
2747 // No-use of the input constructor.
2748 // Optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
2749 return NAN;
2750 } else if (Double.isNaN(y)) {
2751 if (isPosInfinite(x)) {
2752 return constructor.create(Math.copySign(0, real), Double.NaN);
2753 }
2754 if (x == 0) {
2755 return constructor.create(real, Double.NaN);
2756 }
2757 // No-use of the input constructor
2758 return NAN;
2759 } else {
2760 // x && y are finite or infinite.
2761
2762 // Check the safe region.
2763 // The lower and upper bounds have been copied from boost::math::atanh.
2764 // They are different from the safe region for asin and acos.
2765 // x >= SAFE_UPPER: (1-x) == -x
2766 // x <= SAFE_LOWER: 1 - x^2 = 1
2767
2768 if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
2769 // Normal computation within a safe region.
2770
2771 // minus x plus 1: (-x+1)
2772 final double mxp1 = 1 - x;
2773 final double yy = y * y;
2774 // The definition of real component is:
2775 // real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4
2776 // This simplifies by adding 1 and subtracting 1 as a fraction:
2777 // = log(1 + ((x+1)^2+y^2) / ((1-x)^2+y^2) - ((1-x)^2+y^2)/((1-x)^2+y^2) ) / 4
2778 //
2779 // real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4
2780 // imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2
2781 // imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2
2782 // The division is done at the end of the function.
2783 re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
2784 // Modified from boost which does not switch the magnitude of x and y.
2785 // The denominator for atan2 is 1 - x^2 - y^2.
2786 // This can be made more precise if |x| > |y|.
2787 final double numerator = 2 * y;
2788 double denominator;
2789 if (x < y) {
2790 final double tmp = x;
2791 x = y;
2792 y = tmp;
2793 }
2794 // 1 - x is precise if |x| >= 1
2795 if (x >= 1) {
2796 denominator = (1 - x) * (1 + x) - y * y;
2797 } else {
2798 // |x| < 1: Use high precision if possible:
2799 // 1 - x^2 - y^2 = -(x^2 + y^2 - 1)
2800 // Modified from boost to use the custom high precision method.
2801 denominator = -x2y2m1(x, y);
2802 }
2803 im = Math.atan2(numerator, denominator);
2804 } else {
2805 // This section handles exception cases that would normally cause
2806 // underflow or overflow in the main formulas.
2807
2808 // C99. G.7: Special case for imaginary only numbers
2809 if (x == 0) {
2810 if (imaginary == 0) {
2811 return constructor.create(real, imaginary);
2812 }
2813 // atanh(iy) = i atan(y)
2814 return constructor.create(real, Math.atan(imaginary));
2815 }
2816
2817 // Real part:
2818 // real = Math.log1p(4x / ((1-x)^2 + y^2))
2819 // real = Math.log1p(4x / (1 - 2x + x^2 + y^2))
2820 // real = Math.log1p(4x / (1 + x(x-2) + y^2))
2821 // without either overflow or underflow in the squared terms.
2822 if (x >= SAFE_UPPER) {
2823 // (1-x) = -x to machine precision:
2824 // log1p(4x / (x^2 + y^2))
2825 if (isPosInfinite(x) || isPosInfinite(y)) {
2826 re = 0;
2827 } else if (y >= SAFE_UPPER) {
2828 // Big x and y: divide by x*y
2829 re = Math.log1p((4 / y) / (x / y + y / x));
2830 } else if (y > 1) {
2831 // Big x: divide through by x:
2832 re = Math.log1p(4 / (x + y * y / x));
2833 } else {
2834 // Big x small y, as above but neglect y^2/x:
2835 re = Math.log1p(4 / x);
2836 }
2837 } else if (y >= SAFE_UPPER) {
2838 if (x > 1) {
2839 // Big y, medium x, divide through by y:
2840 final double mxp1 = 1 - x;
2841 re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
2842 } else {
2843 // Big y, small x, as above but neglect (1-x)^2/y:
2844 // Note: log1p(v) == v - v^2/2 + v^3/3 ... Taylor series when v is small.
2845 // Here v is so small only the first term matters.
2846 re = 4 * x / y / y;
2847 }
2848 } else if (x == 1) {
2849 // x = 1, small y:
2850 // Special case when x == 1 as (1-x) is invalid.
2851 // Simplify the following formula:
2852 // real = log( sqrt((x+1)^2+y^2) ) / 2 - log( sqrt((1-x)^2+y^2) ) / 2
2853 // = log( sqrt(4+y^2) ) / 2 - log(y) / 2
2854 // if: 4+y^2 -> 4
2855 // = log( 2 ) / 2 - log(y) / 2
2856 // = (log(2) - log(y)) / 2
2857 // Multiply by 2 as it will be divided by 4 at the end.
2858 // C99: if y=0 raises the ‘‘divide-by-zero’’ floating-point exception.
2859 re = 2 * (LN_2 - Math.log(y));
2860 } else {
2861 // Modified from boost which checks y > SAFE_LOWER.
2862 // if y*y -> 0 it will be ignored so always include it.
2863 final double mxp1 = 1 - x;
2864 re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
2865 }
2866
2867 // Imaginary part:
2868 // imag = atan2(2y, (1-x)(1+x) - y^2)
2869 // if x or y are large, then the formula:
2870 // atan2(2y, (1-x)(1+x) - y^2)
2871 // evaluates to +(PI - theta) where theta is negligible compared to PI.
2872 if (x >= SAFE_UPPER || y >= SAFE_UPPER) {
2873 im = Math.PI;
2874 } else if (x <= SAFE_LOWER) {
2875 // (1-x)^2 -> 1
2876 if (y <= SAFE_LOWER) {
2877 // 1 - y^2 -> 1
2878 im = Math.atan2(2 * y, 1);
2879 } else {
2880 im = Math.atan2(2 * y, 1 - y * y);
2881 }
2882 } else {
2883 // Medium x, small y.
2884 // Modified from boost which checks (y == 0) && (x == 1) and sets re = 0.
2885 // This is same as the result from calling atan2(0, 0) so exclude this case.
2886 // 1 - y^2 = 1 so ignore subtracting y^2
2887 im = Math.atan2(2 * y, (1 - x) * (1 + x));
2888 }
2889 }
2890 }
2891
2892 re /= 4;
2893 im /= 2;
2894 return constructor.create(changeSign(re, real),
2895 changeSign(im, imaginary));
2896 }
2897
2898 /**
2899 * Compute {@code x^2 + y^2 - 1} in high precision.
2900 * Assumes that the values x and y can be multiplied without overflow; that
2901 * {@code x >= y}; and both values are positive.
2902 *
2903 * @param x the x value
2904 * @param y the y value
2905 * @return {@code x^2 + y^2 - 1}.
2906 */
2907 private static double x2y2m1(double x, double y) {
2908 // Hull et al used (x-1)*(x+1)+y*y.
2909 // From the paper on page 236:
2910
2911 // If x == 1 there is no cancellation.
2912
2913 // If x > 1, there is also no cancellation, but the argument is now accurate
2914 // only to within a factor of 1 + 3 EPSILSON (note that x – 1 is exact),
2915 // so that error = 3 EPSILON.
2916
2917 // If x < 1, there can be serious cancellation:
2918
2919 // If 4 y^2 < |x^2 – 1| the cancellation is not serious ... the argument is accurate
2920 // only to within a factor of 1 + 4 EPSILSON so that error = 4 EPSILON.
2921
2922 // Otherwise there can be serious cancellation and the relative error in the real part
2923 // could be enormous.
2924
2925 final double xx = x * x;
2926 final double yy = y * y;
2927 // Modify to use high precision before the threshold set by Hull et al.
2928 // This is to preserve the monotonic output of the computation at the switch.
2929 // Set the threshold when x^2 + y^2 is above 0.5 thus subtracting 1 results in a number
2930 // that can be expressed with a higher precision than any number in the range 0.5-1.0
2931 // due to the variable exponent used below 0.5.
2932 if (x < 1 && xx + yy > 0.5) {
2933 // Large relative error.
2934 // This does not use o.a.c.numbers.LinearCombination.value(x, x, y, y, 1, -1).
2935 // It is optimised knowing that:
2936 // - the products are squares
2937 // - the final term is -1 (which does not require split multiplication and addition)
2938 // - The answer will not be NaN as the terms are not NaN components
2939 // - The order is known to be 1 > |x| >= |y|
2940 // The squares are computed using a split multiply algorithm and
2941 // the summation using an extended precision summation algorithm.
2942
2943 // Split x and y as one 26 bits number and one 27 bits number
2944 final double xHigh = splitHigh(x);
2945 final double xLow = x - xHigh;
2946 final double yHigh = splitHigh(y);
2947 final double yLow = y - yHigh;
2948
2949 // Accurate split multiplication x * x and y * y
2950 final double x2Low = squareLow(xLow, xHigh, xx);
2951 final double y2Low = squareLow(yLow, yHigh, yy);
2952
2953 return sumx2y2m1(xx, x2Low, yy, y2Low);
2954 }
2955 return (x - 1) * (x + 1) + yy;
2956 }
2957
2958 /**
2959 * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) create
2960 * a big value from which to derive the two split parts.
2961 * <pre>
2962 * c = (2^s + 1) * a
2963 * a_big = c - a
2964 * a_hi = c - a_big
2965 * a_lo = a - a_hi
2966 * a = a_hi + a_lo
2967 * </pre>
2968 *
2969 * <p>The multiplicand must be odd allowing a p-bit value to be split into
2970 * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
2971 * Combined they have (p-1) bits of significand but the sign bit of {@code a_lo}
2972 * contains a bit of information.
2973 *
2974 * @param a Value.
2975 * @return the high part of the value.
2976 * @see <a href="https://doi.org/10.1007/BF01397083">
2977 * Dekker (1971) A floating-point technique for extending the available precision</a>
2978 */
2979 private static double splitHigh(double a) {
2980 final double c = MULTIPLIER * a;
2981 return c - (c - a);
2982 }
2983
2984 /**
2985 * Compute the round-off from the square of a split number with {@code low} and {@code high}
2986 * components. Uses Dekker's algorithm for split multiplication modified for a square product.
2987 *
2988 * <p>Note: This is candidate to be replaced with {@code Math.fma(x, x, -x * x)} to compute
2989 * the round-off from the square product {@code x * x}. This would remove the requirement
2990 * to compute the split number and make this method redundant. {@code Math.fma} requires
2991 * JDK 9 and FMA hardware support.
2992 *
2993 * @param low Low part of number.
2994 * @param high High part of number.
2995 * @param square Square of the number.
2996 * @return <code>low * low - (((product - high * high) - low * high) - high * low)</code>
2997 * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
2998 * Shewchuk (1997) Theorum 18</a>
2999 */
3000 private static double squareLow(double low, double high, double square) {
3001 final double lh = low * high;
3002 return low * low - (((square - high * high) - lh) - lh);
3003 }
3004
3005 /**
3006 * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
3007 * Dekker's two-sum algorithm. The values are required to be ordered by magnitude:
3008 * {@code |a| >= |b|}.
3009 *
3010 * @param a First part of sum.
3011 * @param b Second part of sum.
3012 * @param x Sum.
3013 * @return <code>b - (x - a)</code>
3014 * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3015 * Shewchuk (1997) Theorum 6</a>
3016 */
3017 private static double fastSumLow(double a, double b, double x) {
3018 // x = a + b
3019 // bVirtual = x - a
3020 // y = b - bVirtual
3021 return b - (x - a);
3022 }
3023
3024 /**
3025 * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
3026 * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
3027 *
3028 * @param a First part of sum.
3029 * @param b Second part of sum.
3030 * @param x Sum.
3031 * @return <code>(a - (x - (x - a))) + (b - (x - a))</code>
3032 * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3033 * Shewchuk (1997) Theorum 7</a>
3034 */
3035 private static double sumLow(double a, double b, double x) {
3036 // x = a + b
3037 // bVirtual = x - a
3038 // aVirtual = x - bVirtual
3039 // bRoundoff = b - bVirtual
3040 // aRoundoff = a - aVirtual
3041 // y = aRoundoff + bRoundoff
3042 final double bVirtual = x - a;
3043 return (a - (x - bVirtual)) + (b - bVirtual);
3044 }
3045
3046 /**
3047 * Sum x^2 + y^2 - 1. It is assumed that {@code y <= x < 1}.
3048 *
3049 * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High] + [-1] + [y2Low, y2High].
3050 *
3051 * @param x2High High part of x^2.
3052 * @param x2Low Low part of x^2.
3053 * @param y2High High part of y^2.
3054 * @param y2Low Low part of y^2.
3055 * @return x^2 + y^2 - 1
3056 * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3057 * Shewchuk (1997) Theorum 12</a>
3058 */
3059 private static double sumx2y2m1(double x2High, double x2Low, double y2High, double y2Low) {
3060 // Let e and f be non-overlapping expansions of components of length m and n.
3061 // The following algorithm will produce a non-overlapping expansion h where the
3062 // sum h_i = e + f and components of h are in increasing order of magnitude.
3063
3064 // Expansion-sum proceeds by a grow-expansion of the first part from one expansion
3065 // into the other, extending its length by 1. The process repeats for the next part
3066 // but the grow-expansion starts at the previous merge position + 1.
3067 // Thus expansion-sum requires mn two-sum operations to merge length m into length n
3068 // resulting in length m+n-1.
3069
3070 // Variables numbered from 1 as per Figure 7 (p.12). The output expansion h is placed
3071 // into e increasing its length for each grow expansion.
3072
3073 // We have two expansions for x^2 and y^2 and the whole number -1.
3074 // Expecting (x^2 + y^2) close to 1 we generate first the intermediate expansion
3075 // (x^2 - 1) moving the result away from 1 where there are sparse floating point
3076 // representations. This is then added to a similar magnitude y^2. Leaving the -1
3077 // until last suffers from 1 ulp rounding errors more often and the requirement
3078 // for a distillation sum to reduce rounding error frequency.
3079
3080 // Note: Do not use the alternative fast-expansion-sum of the parts sorted by magnitude.
3081 // The parts can be ordered with a single comparison into:
3082 // [y2Low, (y2High|x2Low), x2High, -1]
3083 // The fast-two-sum saves 1 fast-two-sum and 3 two-sum operations (21 additions) and
3084 // adds a penalty of a single branch condition.
3085 // However the order in not "strongly non-overlapping" and the fast-expansion-sum
3086 // output will not be strongly non-overlapping. The sum of the output has 1 ulp error
3087 // on random cis numbers approximately 1 in 160 events. This can be removed by a
3088 // distillation two-sum pass over the final expansion as a cost of 1 fast-two-sum and
3089 // 3 two-sum operations! So we use the expansion sum with the same operations and
3090 // no branches.
3091
3092 // q=running sum
3093 double q = x2Low - 1;
3094 double e1 = fastSumLow(-1, x2Low, q);
3095 double e3 = q + x2High;
3096 double e2 = sumLow(q, x2High, e3);
3097
3098 final double f1 = y2Low;
3099 final double f2 = y2High;
3100
3101 // Grow expansion of f1 into e
3102 q = f1 + e1;
3103 e1 = sumLow(f1, e1, q);
3104 double p = q + e2;
3105 e2 = sumLow(q, e2, p);
3106 double e4 = p + e3;
3107 e3 = sumLow(p, e3, e4);
3108
3109 // Grow expansion of f2 into e (only required to start at e2)
3110 q = f2 + e2;
3111 e2 = sumLow(f2, e2, q);
3112 p = q + e3;
3113 e3 = sumLow(q, e3, p);
3114 final double e5 = p + e4;
3115 e4 = sumLow(p, e4, e5);
3116
3117 // Final summation:
3118 // The sum of the parts is within 1 ulp of the true expansion value e:
3119 // |e - sum| < ulp(sum).
3120 // To achieve the exact result requires iteration of a distillation two-sum through
3121 // the expansion until convergence, i.e. no smaller term changes higher terms.
3122 // This requires (n-1) iterations for length n. Here we neglect this as
3123 // although the method is not ensured to be exact is it robust on random
3124 // cis numbers.
3125 return e1 + e2 + e3 + e4 + e5;
3126 }
3127
3128 /**
3129 * Returns the n-th roots of this complex number.
3130 * The nth roots are defined by the formula:
3131 *
3132 * <p>\[ z_k = |z|^{\frac{1}{n}} \left( \cos \left(\phi + \frac{2\pi k}{n} \right) + i \sin \left(\phi + \frac{2\pi k}{n} \right) \right) \]
3133 *
3134 * <p>for \( k=0, 1, \ldots, n-1 \), where \( |z| \) and \( \phi \)
3135 * are respectively the {@link #abs() modulus} and
3136 * {@link #arg() argument} of this complex number.
3137 *
3138 * <p>If one or both parts of this complex number is NaN, a list with all
3139 * all elements set to {@code NaN + i NaN} is returned.</p>
3140 *
3141 * @param n Degree of root.
3142 * @return A list of all {@code n}-th roots of this complex number.
3143 * @throws IllegalArgumentException if {@code n} is zero.
3144 * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Root/">Root</a>
3145 */
3146 public List<Complex> nthRoot(int n) {
3147 if (n == 0) {
3148 throw new IllegalArgumentException("cannot compute zeroth root");
3149 }
3150
3151 final List<Complex> result = new ArrayList<>();
3152
3153 // nth root of abs -- faster / more accurate to use a solver here?
3154 final double nthRootOfAbs = Math.pow(abs(), 1.0 / n);
3155
3156 // Compute nth roots of complex number with k = 0, 1, ... n-1
3157 final double nthPhi = arg() / n;
3158 final double slice = 2 * Math.PI / n;
3159 double innerPart = nthPhi;
3160 for (int k = 0; k < Math.abs(n); k++) {
3161 // inner part
3162 final double realPart = nthRootOfAbs * Math.cos(innerPart);
3163 final double imaginaryPart = nthRootOfAbs * Math.sin(innerPart);
3164 result.add(ofCartesian(realPart, imaginaryPart));
3165 innerPart += slice;
3166 }
3167
3168 return result;
3169 }
3170
3171 /**
3172 * Test for equality with another object. If the other object is a {@code Complex} then a
3173 * comparison is made of the real and imaginary parts; otherwise {@code false} is returned.
3174 *
3175 * <p>If both the real and imaginary parts of two complex numbers
3176 * are exactly the same the two {@code Complex} objects are considered to be equal.
3177 * For this purpose, two {@code double} values are considered to be
3178 * the same if and only if the method {@link Double #doubleToLongBits(double)}
3179 * returns the identical {@code long} value when applied to each.
3180 *
3181 * <p>Note that in most cases, for two instances of class
3182 * {@code Complex}, {@code c1} and {@code c2}, the
3183 * value of {@code c1.equals(c2)} is {@code true} if and only if
3184 *
3185 * <pre>
3186 * {@code c1.getReal() == c2.getReal() && c1.getImaginary() == c2.getImaginary()}</pre>
3187 *
3188 * <p>also has the value {@code true}. However, there are exceptions:
3189 *
3190 * <ul>
3191 * <li>
3192 * Instances that contain {@code NaN} values in the same part
3193 * are considered to be equal for that part, even though {@code Double.NaN == Double.NaN}
3194 * has the value {@code false}.
3195 * </li>
3196 * <li>
3197 * Instances that share a {@code NaN} value in one part
3198 * but have different values in the other part are <em>not</em> considered equal.
3199 * </li>
3200 * <li>
3201 * Instances that contain different representations of zero in the same part
3202 * are <em>not</em> considered to be equal for that part, even though {@code -0.0 == 0.0}
3203 * has the value {@code true}.
3204 * </li>
3205 * </ul>
3206 *
3207 * <p>The behavior is the same as if the components of the two complex numbers were passed
3208 * to {@link java.util.Arrays#equals(double[], double[]) Arrays.equals(double[], double[])}:
3209 *
3210 * <pre>
3211 * Arrays.equals(new double[]{c1.getReal(), c1.getImaginary()},
3212 * new double[]{c2.getReal(), c2.getImaginary()}); </pre>
3213 *
3214 * @param other Object to test for equality with this instance.
3215 * @return {@code true} if the objects are equal, {@code false} if object
3216 * is {@code null}, not an instance of {@code Complex}, or not equal to
3217 * this instance.
3218 * @see Double#doubleToLongBits(double)
3219 * @see java.util.Arrays#equals(double[], double[])
3220 */
3221 @Override
3222 public boolean equals(Object other) {
3223 if (this == other) {
3224 return true;
3225 }
3226 if (other instanceof Complex) {
3227 final Complex c = (Complex) other;
3228 return equals(real, c.real) &&
3229 equals(imaginary, c.imaginary);
3230 }
3231 return false;
3232 }
3233
3234 /**
3235 * Gets a hash code for the complex number.
3236 *
3237 * <p>The behavior is the same as if the components of the complex number were passed
3238 * to {@link java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])}:
3239 *
3240 * <pre>
3241 * {@code Arrays.hashCode(new double[] {getReal(), getImaginary()})}</pre>
3242 *
3243 * @return A hash code value for this object.
3244 * @see java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])
3245 */
3246 @Override
3247 public int hashCode() {
3248 return 31 * (31 + Double.hashCode(real)) + Double.hashCode(imaginary);
3249 }
3250
3251 /**
3252 * Returns a string representation of the complex number.
3253 *
3254 * <p>The string will represent the numeric values of the real and imaginary parts.
3255 * The values are split by a separator and surrounded by parentheses.
3256 * The string can be {@link #parse(String) parsed} to obtain an instance with the same value.
3257 *
3258 * <p>The format for complex number \( x + i y \) is {@code "(x,y)"}, with \( x \) and
3259 * \( y \) converted as if using {@link Double#toString(double)}.
3260 *
3261 * @return A string representation of the complex number.
3262 * @see #parse(String)
3263 * @see Double#toString(double)
3264 */
3265 @Override
3266 public String toString() {
3267 return new StringBuilder(TO_STRING_SIZE)
3268 .append(FORMAT_START)
3269 .append(real).append(FORMAT_SEP)
3270 .append(imaginary)
3271 .append(FORMAT_END)
3272 .toString();
3273 }
3274
3275 /**
3276 * Returns {@code true} if the values are equal according to semantics of
3277 * {@link Double#equals(Object)}.
3278 *
3279 * @param x Value
3280 * @param y Value
3281 * @return {@code Double.valueof(x).equals(Double.valueOf(y))}.
3282 */
3283 private static boolean equals(double x, double y) {
3284 return Double.doubleToLongBits(x) == Double.doubleToLongBits(y);
3285 }
3286
3287 /**
3288 * Check that a value is negative. It must meet all the following conditions:
3289 * <ul>
3290 * <li>it is not {@code NaN},</li>
3291 * <li>it is negative signed,</li>
3292 * </ul>
3293 *
3294 * <p>Note: This is true for negative zero.</p>
3295 *
3296 * @param d Value.
3297 * @return {@code true} if {@code d} is negative.
3298 */
3299 private static boolean negative(double d) {
3300 return d < 0 || Double.doubleToLongBits(d) == NEGATIVE_ZERO_LONG_BITS;
3301 }
3302
3303 /**
3304 * Check that a value is positive infinity. Used to replace {@link Double#isInfinite()}
3305 * when the input value is known to be positive (i.e. in the case where it has been
3306 * set using {@link Math#abs(double)}).
3307 *
3308 * @param d Value.
3309 * @return {@code true} if {@code d} is +inf.
3310 */
3311 private static boolean isPosInfinite(double d) {
3312 return d == Double.POSITIVE_INFINITY;
3313 }
3314
3315 /**
3316 * Check that an absolute value is finite. Used to replace {@link Double#isFinite(double)}
3317 * when the input value is known to be positive (i.e. in the case where it has been
3318 * set using {@link Math#abs(double)}).
3319 *
3320 * @param d Value.
3321 * @return {@code true} if {@code d} is +finite.
3322 */
3323 private static boolean isPosFinite(double d) {
3324 return d <= Double.MAX_VALUE;
3325 }
3326
3327 /**
3328 * Create a complex number given the real and imaginary parts, then multiply by {@code -i}.
3329 * This is used in functions that implement trigonomic identities. It is the functional
3330 * equivalent of:
3331 *
3332 * <pre>
3333 * z = new Complex(real, imaginary).multiplyImaginary(-1);</pre>
3334 *
3335 * @param real Real part.
3336 * @param imaginary Imaginary part.
3337 * @return {@code Complex} object.
3338 */
3339 private static Complex multiplyNegativeI(double real, double imaginary) {
3340 return new Complex(imaginary, -real);
3341 }
3342
3343 /**
3344 * Change the sign of the magnitude based on the signed value.
3345 *
3346 * <p>If the signed value is negative then the result is {@code -magnitude}; otherwise
3347 * return {@code magnitude}.
3348 *
3349 * <p>A signed value of {@code -0.0} is treated as negative. A signed value of {@code NaN}
3350 * is treated as positive.
3351 *
3352 * <p>This is not the same as {@link Math#copySign(double, double)} as this method
3353 * will change the sign based on the signed value rather than copy the sign.
3354 *
3355 * @param magnitude the magnitude
3356 * @param signedValue the signed value
3357 * @return magnitude or -magnitude.
3358 * @see #negative(double)
3359 */
3360 private static double changeSign(double magnitude, double signedValue) {
3361 return negative(signedValue) ? -magnitude : magnitude;
3362 }
3363
3364 /**
3365 * Returns a scale suitable for use with {@link Math#scalb(double, int)} to normalise
3366 * the number to the interval {@code [1, 2)}.
3367 *
3368 * <p>The scale is typically the largest unbiased exponent used in the representation of the
3369 * two numbers. In contrast to {@link Math#getExponent(double)} this handles
3370 * sub-normal numbers by computing the number of leading zeros in the mantissa
3371 * and shifting the unbiased exponent. The result is that for all finite, non-zero,
3372 * numbers {@code a, b}, the magnitude of {@code scalb(x, -getScale(a, b))} is
3373 * always in the range {@code [1, 2)}, where {@code x = max(|a|, |b|)}.
3374 *
3375 * <p>This method is a functional equivalent of the c function ilogb(double) adapted for
3376 * two input arguments.
3377 *
3378 * <p>The result is to be used to scale a complex number using {@link Math#scalb(double, int)}.
3379 * Hence the special case of both zero arguments is handled using the return value for NaN
3380 * as zero cannot be scaled. This is different from {@link Math#getExponent(double)}
3381 * or {@link #getMaxExponent(double, double)}.
3382 *
3383 * <p>Special cases:
3384 *
3385 * <ul>
3386 * <li>If either argument is NaN or infinite, then the result is
3387 * {@link Double#MAX_EXPONENT} + 1.
3388 * <li>If both arguments are zero, then the result is
3389 * {@link Double#MAX_EXPONENT} + 1.
3390 * </ul>
3391 *
3392 * @param a the first value
3393 * @param b the second value
3394 * @return The maximum unbiased exponent of the values to be used for scaling
3395 * @see Math#getExponent(double)
3396 * @see Math#scalb(double, int)
3397 * @see <a href="http://www.cplusplus.com/reference/cmath/ilogb/">ilogb</a>
3398 */
3399 private static int getScale(double a, double b) {
3400 // Only interested in the exponent and mantissa so remove the sign bit
3401 final long x = Double.doubleToRawLongBits(a) & UNSIGN_MASK;
3402 final long y = Double.doubleToRawLongBits(b) & UNSIGN_MASK;
3403 // Only interested in the maximum
3404 final long bits = Math.max(x, y);
3405 // Get the unbiased exponent
3406 int exp = ((int) (bits >>> 52)) - EXPONENT_OFFSET;
3407
3408 // No case to distinguish nan/inf
3409 // Handle sub-normal numbers
3410 if (exp == Double.MIN_EXPONENT - 1) {
3411 // Special case for zero, return as nan/inf to indicate scaling is not possible
3412 if (bits == 0) {
3413 return Double.MAX_EXPONENT + 1;
3414 }
3415 // A sub-normal number has an exponent below -1022. The amount below
3416 // is defined by the number of shifts of the most significant bit in
3417 // the mantissa that is required to get a 1 at position 53 (i.e. as
3418 // if it were a normal number with assumed leading bit)
3419 final long mantissa = bits & MANTISSA_MASK;
3420 exp -= Long.numberOfLeadingZeros(mantissa << 12);
3421 }
3422 return exp;
3423 }
3424
3425 /**
3426 * Returns the largest unbiased exponent used in the representation of the
3427 * two numbers. Special cases:
3428 *
3429 * <ul>
3430 * <li>If either argument is NaN or infinite, then the result is
3431 * {@link Double#MAX_EXPONENT} + 1.
3432 * <li>If both arguments are zero or subnormal, then the result is
3433 * {@link Double#MIN_EXPONENT} -1.
3434 * </ul>
3435 *
3436 * <p>This is used by {@link #divide(double, double, double, double)} as
3437 * a simple detection that a number may overflow if multiplied
3438 * by a value in the interval [1, 2).
3439 *
3440 * @param a the first value
3441 * @param b the second value
3442 * @return The maximum unbiased exponent of the values.
3443 * @see Math#getExponent(double)
3444 * @see #divide(double, double, double, double)
3445 */
3446 private static int getMaxExponent(double a, double b) {
3447 // This could return:
3448 // Math.getExponent(Math.max(Math.abs(a), Math.abs(b)))
3449 // A speed test is required to determine performance.
3450 return Math.max(Math.getExponent(a), Math.getExponent(b));
3451 }
3452
3453 /**
3454 * Checks if both x and y are in the region defined by the minimum and maximum.
3455 *
3456 * @param x x value.
3457 * @param y y value.
3458 * @param min the minimum (exclusive).
3459 * @param max the maximum (exclusive).
3460 * @return true if inside the region.
3461 */
3462 private static boolean inRegion(double x, double y, double min, double max) {
3463 return x < max && x > min && y < max && y > min;
3464 }
3465
3466 /**
3467 * Returns {@code sqrt(x^2 + y^2)} without intermediate overflow or underflow.
3468 *
3469 * <p>Special cases:
3470 * <ul>
3471 * <li>If either argument is infinite, then the result is positive infinity.
3472 * <li>If either argument is NaN and neither argument is infinite, then the result is NaN.
3473 * </ul>
3474 *
3475 * <p>The computed result is expected to be within 1 ulp of the exact result.
3476 *
3477 * <p>This method is a replacement for {@link Math#hypot(double, double)}. There
3478 * will be differences between this method and {@code Math.hypot(double, double)} due
3479 * to the use of a different algorithm to compute the high precision sum of
3480 * {@code x^2 + y^2}. This method has been tested to have a lower maximum error from
3481 * the exact result; any differences are expected to be 1 ULP indicating a rounding
3482 * change in the sum.
3483 *
3484 * <p>JDK9 ported the hypot function to Java for bug JDK-7130085 due to the slow performance
3485 * of the method as a native function. Benchmarks of the Complex class for functions that
3486 * use hypot confirm this is slow pre-Java 9. This implementation outperforms the new faster
3487 * {@code Math.hypot(double, double)} on JDK 11 (LTS). See the Commons numbers examples JMH
3488 * module for benchmarks. Comparisons with alternative implementations indicate
3489 * performance gains are related to edge case handling and elimination of an unpredictable
3490 * branch in the computation of {@code x^2 + y^2}.
3491 *
3492 * <p>This port was adapted from the "Freely Distributable Math Library" hypot function.
3493 * This method only (and not invoked methods within) is distributed under the terms of the
3494 * original notice as shown below:
3495 * <pre>
3496 * ====================================================
3497 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
3498 *
3499 * Developed at SunSoft, a Sun Microsystems, Inc. business.
3500 * Permission to use, copy, modify, and distribute this
3501 * software is freely granted, provided that this notice
3502 * is preserved.
3503 * ====================================================
3504 * </pre>
3505 *
3506 * <p>Note: The fdlibm c code makes use of the language ability to read and write directly
3507 * to the upper and lower 32-bits of the 64-double. The function performs
3508 * checking on the upper 32-bits for the magnitude of the two numbers by accessing
3509 * the exponent and 20 most significant bits of the mantissa. These upper bits
3510 * are manipulated during scaling and then used to perform extended precision
3511 * computation of the sum {@code x^2 + y^2} where the high part of the number has 20-bit
3512 * precision. Manipulation of direct bits has no equivalent in Java
3513 * other than use of {@link Double#doubleToLongBits(double)} and
3514 * {@link Double#longBitsToDouble(long)}. To avoid conversion to and from long and double
3515 * representations this implementation only scales the double representation. The high
3516 * and low parts of a double for the extended precision computation are extracted
3517 * using the method of Dekker (1971) to create two 26-bit numbers. This works for sub-normal
3518 * numbers and reduces the maximum error in comparison to fdlibm hypot which does not
3519 * use a split number algorithm for sub-normal numbers.
3520 *
3521 * @param x Value x
3522 * @param y Value y
3523 * @return sqrt(x^2 + y^2)
3524 * @see Math#hypot(double, double)
3525 * @see <a href="https://www.netlib.org/fdlibm/e_hypot.c">fdlibm e_hypot.c</a>
3526 * @see <a href="https://bugs.java.com/bugdatabase/view_bug.do?bug_id=7130085">JDK-7130085 : Port fdlibm hypot to Java</a>
3527 */
3528 private static double hypot(double x, double y) {
3529 // Differences to the fdlibm reference:
3530 //
3531 // 1. fdlibm orders the two parts using the magnitude of the upper 32-bits.
3532 // This incorrectly orders numbers which differ only in the lower 32-bits.
3533 // This invalidates hypot(x, y) = hypot(y, x) for small sub-normal numbers and a minority
3534 // of cases of normal numbers. This implementation forces the |x| >= |y| order
3535 // using the entire 63-bits of the unsigned doubles to ensure the function
3536 // is commutative.
3537 //
3538 // 2. fdlibm computed scaling by directly writing changes to the exponent bits
3539 // and maintained the high part (ha) during scaling for use in the high
3540 // precision sum x^2 + y^2. Since exponent scaling cannot be applied to sub-normals
3541 // the original version dropped the split number representation for sub-normals
3542 // and can produce maximum errors above 1 ULP for sub-normal numbers.
3543 // This version uses Dekker's method to split the number. This can be applied to
3544 // sub-normals and allows dropping the condition to check for sub-normal numbers
3545 // since all small numbers are handled with a single scaling factor.
3546 // The effect is increased precision for the majority of sub-normal cases where
3547 // the implementations compute a different result.
3548 //
3549 // 3. An alteration is done here to add an 'else if' instead of a second
3550 // 'if' statement. Thus you cannot scale down and up at the same time.
3551 //
3552 // 4. There is no use of the absolute double value. The magnitude comparison is
3553 // performed using the long bit representation. The computation x^2+y^2 is
3554 // insensitive to the sign bit. Thus use of Math.abs(double) is only in edge-case
3555 // branches.
3556 //
3557 // 5. The exponent different to ignore the smaller component has changed from 60 to 54.
3558 //
3559 // Original comments from fdlibm are in c style: /* */
3560 // Extra comments added for reference.
3561 //
3562 // Note that the high 32-bits are compared to constants.
3563 // The lowest 20-bits are the upper bits of the 52-bit mantissa.
3564 // The next 11-bits are the biased exponent. The sign bit has been cleared.
3565 // Scaling factors are powers of two for exact scaling.
3566 // For clarity the values have been refactored to named constants.
3567
3568 // The mask is used to remove the sign bit.
3569 final long xbits = Double.doubleToRawLongBits(x) & UNSIGN_MASK;
3570 final long ybits = Double.doubleToRawLongBits(y) & UNSIGN_MASK;
3571
3572 // Order by magnitude: |a| >= |b|
3573 double a;
3574 double b;
3575 /* High word of x & y */
3576 int ha;
3577 int hb;
3578 if (ybits > xbits) {
3579 a = y;
3580 b = x;
3581 ha = (int) (ybits >>> 32);
3582 hb = (int) (xbits >>> 32);
3583 } else {
3584 a = x;
3585 b = y;
3586 ha = (int) (xbits >>> 32);
3587 hb = (int) (ybits >>> 32);
3588 }
3589
3590 // Check if the smaller part is significant.
3591 // a^2 is computed in extended precision for an effective mantissa of 106-bits.
3592 // An exponent difference of 54 is where b^2 will not overlap a^2.
3593 if ((ha - hb) > EXP_54) {
3594 /* a/b > 2**54 */
3595 // or a is Inf or NaN.
3596 // No addition of a + b for sNaN.
3597 return Math.abs(a);
3598 }
3599
3600 double rescale = 1.0;
3601 if (ha > EXP_500) {
3602 /* a > 2^500 */
3603 if (ha >= EXP_1024) {
3604 /* Inf or NaN */
3605 // Check b is infinite for the IEEE754 result.
3606 // No addition of a + b for sNaN.
3607 return Math.abs(b) == Double.POSITIVE_INFINITY ?
3608 Double.POSITIVE_INFINITY :
3609 Math.abs(a);
3610 }
3611 /* scale a and b by 2^-600 */
3612 // Before scaling: a in [2^500, 2^1023].
3613 // After scaling: a in [2^-100, 2^423].
3614 // After scaling: b in [2^-154, 2^423].
3615 a *= TWO_POW_NEG_600;
3616 b *= TWO_POW_NEG_600;
3617 rescale = TWO_POW_600;
3618 } else if (hb < EXP_NEG_500) {
3619 // No special handling of sub-normals.
3620 // These do not matter when we do not manipulate the exponent bits
3621 // for scaling the split representation.
3622
3623 // Intentional comparison with zero.
3624 if (b == 0) {
3625 return Math.abs(a);
3626 }
3627
3628 /* scale a and b by 2^600 */
3629 // Effective min exponent of a sub-normal = -1022 - 52 = -1074.
3630 // Before scaling: b in [2^-1074, 2^-501].
3631 // After scaling: b in [2^-474, 2^99].
3632 // After scaling: a in [2^-474, 2^153].
3633 a *= TWO_POW_600;
3634 b *= TWO_POW_600;
3635 rescale = TWO_POW_NEG_600;
3636 }
3637
3638 // High precision x^2 + y^2
3639 return Math.sqrt(x2y2(a, b)) * rescale;
3640 }
3641
3642 /**
3643 * Return {@code x^2 + y^2} with high accuracy.
3644 *
3645 * <p>It is assumed that {@code 2^500 > |x| >= |y| > 2^-500}. Thus there will be no
3646 * overflow or underflow of the result. The inputs are not assumed to be unsigned.
3647 *
3648 * <p>The computation is performed using Dekker's method for extended precision
3649 * multiplication of x and y and then summation of the extended precision squares.
3650 *
3651 * @param x Value x.
3652 * @param y Value y
3653 * @return x^2 + y^2
3654 * @see <a href="https://doi.org/10.1007/BF01397083">
3655 * Dekker (1971) A floating-point technique for extending the available precision</a>
3656 */
3657 private static double x2y2(double x, double y) {
3658 // Note:
3659 // This method is different from the high-accuracy summation used in fdlibm for hypot.
3660 // The summation could be any valid computation of x^2+y^2. However since this follows
3661 // the re-scaling logic in hypot(x, y) the use of high precision has relatively
3662 // less performance overhead than if used without scaling.
3663 // The Dekker algorithm is branchless for better performance
3664 // than the fdlibm method with a maximum ULP error of approximately 0.86.
3665 //
3666 // See NUMBERS-143 for analysis.
3667
3668 // Do a Dekker summation of double length products x*x and y*y
3669 // (10 multiply and 20 additions).
3670 final double xx = x * x;
3671 final double yy = y * y;
3672 // Compute the round-off from the products.
3673 // With FMA hardware support in JDK 9+ this can be replaced with the much faster:
3674 // xxLow = Math.fma(x, x, -xx)
3675 // yyLow = Math.fma(y, y, -yy)
3676 // Dekker mul12
3677 final double xHigh = splitHigh(x);
3678 final double xLow = x - xHigh;
3679 final double xxLow = squareLow(xLow, xHigh, xx);
3680 // Dekker mul12
3681 final double yHigh = splitHigh(y);
3682 final double yLow = y - yHigh;
3683 final double yyLow = squareLow(yLow, yHigh, yy);
3684 // Dekker add2
3685 final double r = xx + yy;
3686 // Note: The order is important. Assume xx > yy and drop Dekker's conditional
3687 // check for which is the greater magnitude.
3688 // s = xx - r + yy + yyLow + xxLow
3689 // z = r + s
3690 // zz = r - z + s
3691 // Here we compute z inline and ignore computing the round-off zz.
3692 // Note: The round-off could be used with Dekker's sqrt2 method.
3693 // That adds 7 multiply, 1 division and 19 additions doubling the cost
3694 // and reducing error to < 0.5 ulp for the final sqrt.
3695 return xx - r + yy + yyLow + xxLow + r;
3696 }
3697 }