View Javadoc
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 org.junit.jupiter.api.Assertions;
21  import org.junit.jupiter.api.Test;
22  
23  import java.math.BigDecimal;
24  import java.util.function.BiFunction;
25  import java.util.function.UnaryOperator;
26  
27  /**
28   * Edge case tests for the functions defined by the C.99 standard for complex numbers
29   * defined in ISO/IEC 9899, Annex G.
30   *
31   * <p>The test contained here are specifically written to target edge cases of finite valued
32   * input values that cause overflow/underflow during the computation.
33   *
34   * <p>The test data is generated from a known implementation of the standard.
35   *
36   * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards">
37   *    ISO/IEC 9899 - Programming languages - C</a>
38   */
39  class ComplexEdgeCaseTest {
40      private static final double inf = Double.POSITIVE_INFINITY;
41      private static final double nan = Double.NaN;
42  
43      /**
44       * Assert the operation on the complex number is equal to the expected value.
45       *
46       * <p>The results are are considered equal if there are no floating-point values between them.
47       *
48       * @param a Real part.
49       * @param b Imaginary part.
50       * @param name The operation name.
51       * @param operation The operation.
52       * @param x Expected real part.
53       * @param y Expected imaginary part.
54       */
55      private static void assertComplex(double a, double b,
56              String name, UnaryOperator<Complex> operation,
57              double x, double y) {
58          assertComplex(a, b, name, operation, x, y, 1);
59      }
60  
61      /**
62       * Assert the operation on the complex number is equal to the expected value.
63       *
64       * <p>The results are considered equal within the provided units of least
65       * precision. The maximum count of numbers allowed between the two values is
66       * {@code maxUlps - 1}.
67       *
68       * @param a Real part.
69       * @param b Imaginary part.
70       * @param name The operation name.
71       * @param operation The operation.
72       * @param x Expected real part.
73       * @param y Expected imaginary part.
74       * @param maxUlps the maximum units of least precision between the two values
75       */
76      private static void assertComplex(double a, double b,
77              String name, UnaryOperator<Complex> operation,
78              double x, double y, long maxUlps) {
79          final Complex c = Complex.ofCartesian(a, b);
80          final Complex e = Complex.ofCartesian(x, y);
81          CReferenceTest.assertComplex(c, name, operation, e, maxUlps);
82      }
83  
84      /**
85       * Assert the operation on the complex numbers is equal to the expected value.
86       *
87       * <p>The results are considered equal if there are no floating-point values between them.
88       *
89       * @param a Real part of first number.
90       * @param b Imaginary part of first number.
91       * @param c Real part of second number.
92       * @param d Imaginary part of second number.
93       * @param name The operation name.
94       * @param operation The operation.
95       * @param x Expected real part.
96       * @param y Expected imaginary part.
97       */
98      // CHECKSTYLE: stop ParameterNumberCheck
99      private static void assertComplex(double a, double b, double c, double d,
100             String name, BiFunction<Complex, Complex, Complex> operation,
101             double x, double y) {
102         assertComplex(a, b, c, d, name, operation, x, y, 1);
103     }
104 
105     /**
106      * Assert the operation on the complex numbers is equal to the expected value.
107      *
108      * <p>The results are considered equal within the provided units of least
109      * precision. The maximum count of numbers allowed between the two values is
110      * {@code maxUlps - 1}.
111      *
112      * @param a Real part of first number.
113      * @param b Imaginary part of first number.
114      * @param c Real part of second number.
115      * @param d Imaginary part of second number.
116      * @param name The operation name
117      * @param operation the operation
118      * @param x Expected real part.
119      * @param y Expected imaginary part.
120      * @param maxUlps the maximum units of least precision between the two values
121      */
122     private static void assertComplex(double a, double b, double c, double d,
123             String name, BiFunction<Complex, Complex, Complex> operation,
124             double x, double y, long maxUlps) {
125         final Complex c1 = Complex.ofCartesian(a, b);
126         final Complex c2 = Complex.ofCartesian(c, d);
127         final Complex e = Complex.ofCartesian(x, y);
128         CReferenceTest.assertComplex(c1, c2, name, operation, e, maxUlps);
129     }
130 
131     @Test
132     void testAcos() {
133         // acos(z) = (pi / 2) + i ln(iz + sqrt(1 - z^2))
134         final String name = "acos";
135         final UnaryOperator<Complex> operation = Complex::acos;
136 
137         // Edge cases are when values are big but not infinite and small but not zero.
138         // Big and small are set using the limits in atanh.
139         // A medium value is used to test outside the range of the CReferenceTest.
140         // The results have been generated using g++ -std=c++11 acos.
141         // xp1 * xm1 will overflow:
142         final double huge = Math.sqrt(Double.MAX_VALUE) * 2;
143         final double big = Math.sqrt(Double.MAX_VALUE) / 8;
144         final double medium = 100;
145         final double small = Math.sqrt(Double.MIN_NORMAL) * 4;
146         assertComplex(huge, big, name, operation, 0.06241880999595735, -356.27960012801969);
147         assertComplex(huge, medium, name, operation, 3.7291703656001039e-153, -356.27765080781188);
148         assertComplex(huge, small, name, operation, 2.2250738585072019e-308, -356.27765080781188);
149         assertComplex(big, big, name, operation, 0.78539816339744828, -353.85163567585209);
150         assertComplex(big, medium, name, operation, 5.9666725849601662e-152, -353.50506208557209);
151         assertComplex(big, small, name, operation, 3.560118173611523e-307, -353.50506208557209);
152         assertComplex(medium, big, name, operation, 1.5707963267948966, -353.50506208557209);
153         assertComplex(medium, medium, name, operation, 0.78541066339744181, -5.6448909570623842);
154         assertComplex(medium, small, name, operation, 5.9669709409662999e-156, -5.298292365610485);
155         assertComplex(small, big, name, operation, 1.5707963267948966, -353.50506208557209);
156         assertComplex(small, medium, name, operation, 1.5707963267948966, -5.2983423656105888);
157         assertComplex(small, small, name, operation, 1.5707963267948966, -5.9666725849601654e-154);
158         // Additional cases to achieve full coverage
159         // xm1 = 0
160         assertComplex(1, small, name, operation, 2.4426773395109241e-77, -2.4426773395109241e-77);
161         // https://svn.boost.org/trac10/ticket/7290
162         assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 2.4252018043912224e-196, -0.00023604834149293664);
163     }
164 
165     // acosh is defined by acos so is not tested
166 
167     @Test
168     void testAsin() {
169         // asin(z) = -i (ln(iz + sqrt(1 - z^2)))
170         final String name = "asin";
171         final UnaryOperator<Complex> operation = Complex::asin;
172 
173         // This method is essentially the same as acos and the edge cases are the same.
174         // The results have been generated using g++ -std=c++11 asin.
175         final double huge = Math.sqrt(Double.MAX_VALUE) * 2;
176         final double big = Math.sqrt(Double.MAX_VALUE) / 8;
177         final double medium = 100;
178         final double small = Math.sqrt(Double.MIN_NORMAL) * 4;
179         assertComplex(huge, big, name, operation, 1.5083775167989393, 356.27960012801969);
180         assertComplex(huge, medium, name, operation, 1.5707963267948966, 356.27765080781188);
181         assertComplex(huge, small, name, operation, 1.5707963267948966, 356.27765080781188);
182         assertComplex(big, big, name, operation, 0.78539816339744828, 353.85163567585209);
183         assertComplex(big, medium, name, operation, 1.5707963267948966, 353.50506208557209);
184         assertComplex(big, small, name, operation, 1.5707963267948966, 353.50506208557209);
185         assertComplex(medium, big, name, operation, 5.9666725849601662e-152, 353.50506208557209);
186         assertComplex(medium, medium, name, operation, 0.78538566339745486, 5.6448909570623842);
187         assertComplex(medium, small, name, operation, 1.5707963267948966, 5.298292365610485);
188         assertComplex(small, big, name, operation, 3.560118173611523e-307, 353.50506208557209);
189         assertComplex(small, medium, name, operation, 5.9663742737040751e-156, 5.2983423656105888);
190         assertComplex(small, small, name, operation, 5.9666725849601654e-154, 5.9666725849601654e-154);
191         // Additional cases to achieve full coverage
192         // xm1 = 0
193         assertComplex(1, small, name, operation, 1.5707963267948966, 2.4426773395109241e-77);
194         // https://svn.boost.org/trac10/ticket/7290
195         assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 1.5707963267948966, 0.00023604834149293664);
196     }
197 
198     // asinh is defined by asin so is not tested
199 
200     @Test
201     void testAtanh() {
202         // atanh(z) = (1/2) ln((1 + z) / (1 - z))
203         // Odd function: negative real cases defined by positive real cases
204         final String name = "atanh";
205         final UnaryOperator<Complex> operation = Complex::atanh;
206 
207         // Edge cases are when values are big but not infinite and small but not zero.
208         // Big and small are set using the limits in atanh.
209         // A medium value is used to test outside the range of the CReferenceTest.
210         // It hits an edge case when x is big and y > 1.
211         // The results have been generated using g++ -std=c++11 atanh.
212         final double big = Math.sqrt(Double.MAX_VALUE) / 2;
213         final double medium = 100;
214         final double small = Math.sqrt(Double.MIN_NORMAL) * 2;
215         assertComplex(big, big, name, operation, 7.4583407312002067e-155, 1.5707963267948966);
216         assertComplex(big, medium, name, operation, 1.4916681462400417e-154, 1.5707963267948966);
217         assertComplex(big, small, name, operation, 1.4916681462400417e-154, 1.5707963267948966);
218         assertComplex(medium, big, name, operation, 2.225073858507202e-306, 1.5707963267948966);
219         assertComplex(medium, medium, name, operation, 0.0049999166641667555, 1.5657962434640633);
220         assertComplex(medium, small, name, operation, 0.010000333353334761, 1.5707963267948966);
221         assertComplex(small, big, name, operation, 0, 1.5707963267948966);
222         assertComplex(small, medium, name, operation, 2.9830379886812147e-158, 1.5607966601082315);
223         assertComplex(small, small, name, operation, 2.9833362924800827e-154, 2.9833362924800827e-154);
224         // Additional cases to achieve full coverage
225         assertComplex(inf, big, name, operation, 0, 1.5707963267948966);
226         assertComplex(big, inf, name, operation, 0, 1.5707963267948966);
227     }
228 
229     @Test
230     void testCosh() {
231         // cosh(a + b i) = cosh(a)cos(b) + i sinh(a)sin(b)
232         // Even function: negative real cases defined by positive real cases
233         final String name = "cosh";
234         final UnaryOperator<Complex> operation = Complex::cosh;
235 
236         // Implementation defers to java.util.Math.
237         // Hit edge cases for extreme values.
238         final double big = Double.MAX_VALUE;
239         final double medium = 2;
240         final double small = Double.MIN_NORMAL;
241         assertComplex(big, big, name, operation, -inf, inf);
242         assertComplex(big, medium, name, operation, -inf, inf);
243         assertComplex(big, small, name, operation, inf, inf);
244         assertComplex(medium, big, name, operation, -3.7621493762972804, 0.017996317370418576);
245         assertComplex(medium, medium, name, operation, -1.5656258353157435, 3.297894836311237);
246         assertComplex(medium, small, name, operation, 3.7621956910836314, 8.0700322819551687e-308);
247         assertComplex(small, big, name, operation, -0.99998768942655991, 1.1040715888508271e-310);
248         assertComplex(small, medium, name, operation, -0.41614683654714241, 2.0232539340376892e-308);
249         assertComplex(small, small, name, operation, 1, 0);
250 
251         // Overflow test.
252         // Based on MATH-901 discussion of FastMath functionality.
253         // https://issues.apache.org/jira/browse/MATH-901#comment-13500669
254         // sinh(x)/cosh(x) can be approximated by exp(x) but must be overflow safe.
255 
256         // sinh(x) = sign(x) * e^|x| / 2 when x is large.
257         // cosh(x) = e^|x| / 2 when x is large.
258         // Thus e^|x| can overflow but e^|x| / 2 may not.
259         // (e^|x| / 2) * sin/cos will always be smaller.
260         final double tiny = Double.MIN_VALUE;
261         final double x = 709.783;
262         Assertions.assertEquals(inf, Math.exp(x));
263         // As computed by GNU g++
264         assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
265         assertComplex(-x, 0, name, operation, 8.9910466927705402e+307, -0.0);
266         // sub-normal number x:
267         // cos(x) = 1 => real = (e^|x| / 2)
268         // sin(x) = x => imaginary = x * (e^|x| / 2)
269         assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
270         assertComplex(-x, small, name, operation, 8.9910466927705402e+307, -2.0005742956701358);
271         assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
272         assertComplex(-x, tiny, name, operation, 8.9910466927705402e+307, -4.4421672910524807e-16);
273         // Should not overflow imaginary.
274         assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
275         assertComplex(-2 * x, tiny, name, operation, inf, -7.9879467061901743e+292);
276         // Test when large enough to overflow any non-zero value to infinity. Result should be
277         // as if x was infinite and y was finite.
278         assertComplex(3 * x, tiny, name, operation, inf, inf);
279         assertComplex(-3 * x, tiny, name, operation, inf, -inf);
280         // pi / 2 x:
281         // cos(x) = ~0 => real = x * (e^|x| / 2)
282         // sin(x) = ~1 => imaginary = (e^|x| / 2)
283         final double pi2 = Math.PI / 2;
284         assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
285         assertComplex(-x, pi2, name, operation, 5.5054282766429199e+291, -8.9910466927705402e+307);
286         assertComplex(2 * x, pi2, name, operation, inf, inf);
287         assertComplex(-2 * x, pi2, name, operation, inf, -inf);
288         // Test when large enough to overflow any non-zero value to infinity. Result should be
289         // as if x was infinite and y was finite.
290         assertComplex(3 * x, pi2, name, operation, inf, inf);
291         assertComplex(-3 * x, pi2, name, operation, inf, -inf);
292     }
293 
294     @Test
295     void testSinh() {
296         // sinh(a + b i) = sinh(a)cos(b) + i cosh(a)sin(b)
297         // Odd function: negative real cases defined by positive real cases
298         final String name = "sinh";
299         final UnaryOperator<Complex> operation = Complex::sinh;
300 
301         // Implementation defers to java.util.Math.
302         // Hit edge cases for extreme values.
303         final double big = Double.MAX_VALUE;
304         final double medium = 2;
305         final double small = Double.MIN_NORMAL;
306         assertComplex(big, big, name, operation, -inf, inf);
307         assertComplex(big, medium, name, operation, -inf, inf);
308         assertComplex(big, small, name, operation, inf, inf);
309         assertComplex(medium, big, name, operation, -3.6268157591156114, 0.018667844927220067);
310         assertComplex(medium, medium, name, operation, -1.5093064853236158, 3.4209548611170133);
311         assertComplex(medium, small, name, operation, 3.626860407847019, 8.3711632828186228e-308);
312         assertComplex(small, big, name, operation, -2.2250464665720564e-308, 0.004961954789184062);
313         assertComplex(small, medium, name, operation, -9.2595744730151568e-309, 0.90929742682568171);
314         assertComplex(small, small, name, operation, 2.2250738585072014e-308, 2.2250738585072014e-308);
315 
316         // Overflow test.
317         // As per cosh with sign changes to real and imaginary
318 
319         // sinh(x) = sign(x) * e^|x| / 2 when x is large.
320         // cosh(x) = e^|x| / 2 when x is large.
321         // Thus e^|x| can overflow but e^|x| / 2 may not.
322         // sinh(x) * sin/cos will always be smaller.
323         final double tiny = Double.MIN_VALUE;
324         final double x = 709.783;
325         Assertions.assertEquals(inf, Math.exp(x));
326         // As computed by GNU g++
327         assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
328         assertComplex(-x, 0, name, operation, -8.9910466927705402e+307, 0.0);
329         // sub-normal number:
330         // cos(x) = 1 => real = (e^|x| / 2)
331         // sin(x) = x => imaginary = x * (e^|x| / 2)
332         assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
333         assertComplex(-x, small, name, operation, -8.9910466927705402e+307, 2.0005742956701358);
334         assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
335         assertComplex(-x, tiny, name, operation, -8.9910466927705402e+307, 4.4421672910524807e-16);
336         // Should not overflow imaginary.
337         assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
338         assertComplex(-2 * x, tiny, name, operation, -inf, 7.9879467061901743e+292);
339         // Test when large enough to overflow any non-zero value to infinity. Result should be
340         // as if x was infinite and y was finite.
341         assertComplex(3 * x, tiny, name, operation, inf, inf);
342         assertComplex(-3 * x, tiny, name, operation, -inf, inf);
343         // pi / 2 x:
344         // cos(x) = ~0 => real = x * (e^|x| / 2)
345         // sin(x) = ~1 => imaginary = (e^|x| / 2)
346         final double pi2 = Math.PI / 2;
347         assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
348         assertComplex(-x, pi2, name, operation, -5.5054282766429199e+291, 8.9910466927705402e+307);
349         assertComplex(2 * x, pi2, name, operation, inf, inf);
350         assertComplex(-2 * x, pi2, name, operation, -inf, inf);
351         // Test when large enough to overflow any non-zero value to infinity. Result should be
352         // as if x was infinite and y was finite.
353         assertComplex(3 * x, pi2, name, operation, inf, inf);
354         assertComplex(-3 * x, pi2, name, operation, -inf, inf);
355     }
356 
357     @Test
358     void testTanh() {
359         // tan(a + b i) = sinh(2a)/(cosh(2a)+cos(2b)) + i [sin(2b)/(cosh(2a)+cos(2b))]
360         // Odd function: negative real cases defined by positive real cases
361         final String name = "tanh";
362         final UnaryOperator<Complex> operation = Complex::tanh;
363 
364         // Overflow on 2b:
365         // cos(2b) = cos(inf) = NaN
366         // sin(2b) = sin(inf) = NaN
367         assertComplex(1, Double.MAX_VALUE, name, operation, 0.76160203106265523, -0.0020838895895863505);
368 
369         // Underflow on 2b:
370         // cos(2b) -> 1
371         // sin(2b) -> 0
372         assertComplex(1, Double.MIN_NORMAL, name, operation, 0.76159415595576485, 9.344739287691424e-309);
373         assertComplex(1, Double.MIN_VALUE, name, operation, 0.76159415595576485, 0);
374 
375         // Overflow on 2a:
376         // sinh(2a) = sinh(inf) = inf
377         // cosh(2a) = cosh(inf) = inf
378         // Test all sign variants as this execution path to treat real as infinite
379         // is not tested else where.
380         assertComplex(Double.MAX_VALUE, 1, name, operation, 1, 0.0);
381         assertComplex(Double.MAX_VALUE, -1, name, operation, 1, -0.0);
382         assertComplex(-Double.MAX_VALUE, 1, name, operation, -1, 0.0);
383         assertComplex(-Double.MAX_VALUE, -1, name, operation, -1, -0.0);
384 
385         // Underflow on 2a:
386         // sinh(2a) -> 0
387         // cosh(2a) -> 0
388         assertComplex(Double.MIN_NORMAL, 1, name, operation, 7.6220323800193346e-308, 1.5574077246549021);
389         assertComplex(Double.MIN_VALUE, 1, name, operation, 1.4821969375237396e-323, 1.5574077246549021);
390 
391         // Underflow test.
392         // sinh(x) can be approximated by exp(x) but must be overflow safe.
393         // im = 2 sin(2y) / e^2|x|
394         // This can be computed when e^2|x| only just overflows.
395         // Set a case where e^2|x| overflows but the imaginary can be computed
396         double x = 709.783 / 2;
397         double y = Math.PI / 4;
398         Assertions.assertEquals(1.0, Math.sin(2 * y), 1e-16);
399         Assertions.assertEquals(Double.POSITIVE_INFINITY, Math.exp(2 * x));
400         // As computed by GNU g++
401         assertComplex(x, y, name, operation, 1, 1.1122175583895849e-308);
402     }
403 
404     @Test
405     void testExp() {
406         final String name = "exp";
407         final UnaryOperator<Complex> operation = Complex::exp;
408 
409         // exp(a + b i) = exp(a) (cos(b) + i sin(b))
410 
411         // Overflow if exp(a) == inf
412         assertComplex(1000, 0, name, operation, inf, 0.0);
413         assertComplex(1000, 1, name, operation, inf, inf);
414         assertComplex(1000, 2, name, operation, -inf, inf);
415         assertComplex(1000, 3, name, operation, -inf, inf);
416         assertComplex(1000, 4, name, operation, -inf, -inf);
417 
418         // Underflow if exp(a) == 0
419         assertComplex(-1000, 0, name, operation, 0.0, 0.0);
420         assertComplex(-1000, 1, name, operation, 0.0, 0.0);
421         assertComplex(-1000, 2, name, operation, -0.0, 0.0);
422         assertComplex(-1000, 3, name, operation, -0.0, 0.0);
423         assertComplex(-1000, 4, name, operation, -0.0, -0.0);
424     }
425 
426     @Test
427     void testLog() {
428         final String name = "log";
429         final UnaryOperator<Complex> operation = Complex::log;
430 
431         // ln(a + b i) = ln(|a + b i|) + i arg(a + b i)
432         // |a + b i| = sqrt(a^2 + b^2)
433         // arg(a + b i) = Math.atan2(imaginary, real)
434 
435         // Overflow if sqrt(a^2 + b^2) == inf.
436         // Matlab computes this.
437         assertComplex(-Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 7.101292864836639e2, Math.PI * 3 / 4);
438         assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 7.101292864836639e2, Math.PI / 4);
439         assertComplex(-Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.896613990462929);
440         assertComplex(Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.449786631268641e-1, 2);
441 
442         // Underflow if sqrt(a^2 + b^2) -> 0
443         assertComplex(-Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, -708.04984494198413, 2.3561944901923448);
444         assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, -708.04984494198413, 0.78539816339744828);
445         // Math.hypot(min, min) = min.
446         // To compute the expected result do scaling of the actual hypot = sqrt(2).
447         // log(a/n) = log(a) - log(n)
448         // n = 2^1074 => log(a) - log(2) * 1074
449         double expected = Math.log(Math.sqrt(2)) - Math.log(2) * 1074;
450         assertComplex(-Double.MIN_VALUE, Double.MIN_VALUE, name, operation, expected, Math.atan2(1, -1));
451         expected = Math.log(Math.sqrt(5)) - Math.log(2) * 1074;
452         assertComplex(-Double.MIN_VALUE, 2 * Double.MIN_VALUE, name, operation, expected, Math.atan2(2, -1));
453 
454         // Imprecision if sqrt(a^2 + b^2) == 1 as log(1) is 0.
455         // Method should switch to using log1p(x^2 + x^2 - 1) * 0.5.
456 
457         // In the following:
458         // max = max(real, imaginary)
459         // min = min(real, imaginary)
460 
461         // No cancellation error when max > 1
462 
463         assertLog(1.0001, Math.sqrt(1.2 - 1.0001 * 1.0001), 1);
464         assertLog(1.0001, Math.sqrt(1.1 - 1.0001 * 1.0001), 1);
465         assertLog(1.0001, Math.sqrt(1.02 - 1.0001 * 1.0001), 0);
466         assertLog(1.0001, Math.sqrt(1.01 - 1.0001 * 1.0001), 0);
467 
468         // Cancellation error when max < 1.
469 
470         // Hard: 4 * min^2 < |max^2 - 1|
471         // Gets harder as max is further from 1
472         assertLog(0.99, 0.00001, 0);
473         assertLog(0.95, 0.00001, 0);
474         assertLog(0.9, 0.00001, 0);
475         assertLog(0.85, 0.00001, 0);
476         assertLog(0.8, 0.00001, 0);
477         assertLog(0.75, 0.00001, 0);
478         // At this point the log function does not use high precision computation
479         assertLog(0.7, 0.00001, 2);
480 
481         // Very hard: 4 * min^2 > |max^2 - 1|
482 
483         // Radius 0.99
484         assertLog(0.97, Math.sqrt(0.99 - 0.97 * 0.97), 0);
485         // Radius 1.01
486         assertLog(0.97, Math.sqrt(1.01 - 0.97 * 0.97), 0);
487 
488         // Massive relative error
489         // Radius 0.9999
490         assertLog(0.97, Math.sqrt(0.9999 - 0.97 * 0.97), 0);
491 
492         // polar numbers on a 1/8 circle with a magnitude close to 1.
493         final int steps = 20;
494         final double[] magnitude = {0.999, 1.0, 1.001};
495         final int[] ulps = {0, 0, 1};
496         for (int j = 0; j < magnitude.length; j++) {
497             for (int i = 1; i <= steps; i++) {
498                 final double theta = i * Math.PI / (4 * steps);
499                 assertLog(magnitude[j] * Math.sin(theta), magnitude[j] * Math.cos(theta), ulps[j]);
500             }
501         }
502 
503         // cis numbers using an increasingly smaller angle
504         double theta = Math.PI / (4 * steps);
505         while (theta > 0) {
506             theta /= 2;
507             assertLog(Math.sin(theta), Math.cos(theta), 0);
508         }
509 
510         // Extreme cases.
511         final double up1 = Math.nextUp(1.0);
512         final double down1 = Math.nextDown(1.0);
513         assertLog(down1, Double.MIN_NORMAL, 0);
514         assertLog(down1, Double.MIN_VALUE, 0);
515         // No use of high-precision computation
516         assertLog(up1, Double.MIN_NORMAL, 2);
517         assertLog(up1, Double.MIN_VALUE, 2);
518 
519         // Add some cases known to fail without very high precision computation.
520         // These have been found using randomly generated cis numbers and the
521         // previous Dekker split-summation algorithm:
522         // theta = rng.nextDouble()
523         // x = Math.sin(theta)
524         // y = Math.cos(theta)
525         // Easy: <16 ulps with the Dekker summation
526         assertLog(0.007640392270319105, 0.9999708117770016, 0);
527         assertLog(0.40158433204881533, 0.9158220483548684, 0);
528         assertLog(0.13258789214774552, 0.9911712520325727, 0);
529         assertLog(0.2552206803398717, 0.9668828286441191, 0);
530         // Hard: >1024 ulps with the Dekker summation
531         assertLog(0.4650816500945186, 0.8852677892848919, 0);
532         assertLog(0.06548693057069123, 0.9978534270745526, 0);
533         assertLog(0.08223027214657339, 0.9966133564942327, 0);
534         assertLog(0.06548693057069123, 0.9978534270745526, 0);
535         assertLog(0.04590800199633988, 0.9989456718724518, 0);
536         assertLog(0.3019636508581243, 0.9533194394118022, 0);
537     }
538 
539     /**
540      * Assert the Complex log function using BigDecimal to compute the field norm
541      * {@code x*x + y*y} and then {@link Math#log1p(double)} to compute the log of
542      * the modulus \ using {@code 0.5 * log1p(x*x + y*y - 1)}. This test is for the
543      * extreme case for performance around {@code sqrt(x*x + y*y) = 1} where using
544      * {@link Math#log(double)} will fail dramatically.
545      *
546      * @param x the real value of the complex
547      * @param y the imaginary value of the complex
548      * @param maxUlps the maximum units of least precision between the two values
549      */
550     private static void assertLog(double x, double y, long maxUlps) {
551         // Compute the best value we can
552         final BigDecimal bx = new BigDecimal(x);
553         final BigDecimal by = new BigDecimal(y);
554         final BigDecimal exact = bx.multiply(bx).add(by.multiply(by)).subtract(BigDecimal.ONE);
555         final double real = 0.5 * Math.log1p(exact.doubleValue());
556         final double imag = Math.atan2(y, x);
557         assertComplex(x, y, "log", Complex::log, real, imag, maxUlps);
558     }
559 
560     @Test
561     void testSqrt() {
562         final String name = "sqrt";
563         final UnaryOperator<Complex> operation = Complex::sqrt;
564 
565         // Test real/imaginary only numbers satisfy the definition using polar coordinates:
566         // real = sqrt(abs()) * Math.cos(arg() / 2)
567         // imag = sqrt(abs()) * Math.sin(arg() / 2)
568         // However direct use of sin/cos will result in incorrect results due floating-point error.
569         // This test asserts the on the closest result to the exact answer which is possible
570         // if not using a simple polar computation.
571         // Note: If this test fails in the set-up assertions it is due to a change in the
572         // precision of java.util.Math.
573 
574         // For positive real-only the argument is +/-0.
575         // For negative real-only the argument is +/-pi.
576         Assertions.assertEquals(0, Math.atan2(0, 1));
577         Assertions.assertEquals(Math.PI, Math.atan2(0, -1));
578         // In both cases the trigonomic functions should be exact but
579         // cos(pi/2) cannot be as pi/2 is not exact.
580         final double cosArgRe = 1.0;
581         final double sinArgRe = 0.0;
582         Assertions.assertNotEquals(0.0, Math.cos(Math.PI / 2), "Expected cos(pi/2) to be non-zero");
583         Assertions.assertEquals(0.0, Math.cos(Math.PI / 2), 6.123233995736766e-17);
584         // For imaginary-only the argument is Math.atan2(y, 0) = +/- pi/2.
585         Assertions.assertEquals(Math.PI / 2, Math.atan2(1, 0));
586         Assertions.assertEquals(-Math.PI / 2, Math.atan2(-1, 0));
587         // There is 1 ULP difference in the result of cos/sin of pi/4.
588         // It should be sqrt(2) / 2 for both.
589         final double cosArgIm = Math.cos(Math.PI / 4);
590         final double sinArgIm = Math.sin(Math.PI / 4);
591         final double root2over2 = Math.sqrt(2) / 2;
592         final double ulp = Math.ulp(cosArgIm);
593         Assertions.assertNotEquals(cosArgIm, sinArgIm, "Expected cos(pi/4) to not exactly equal sin(pi/4)");
594         Assertions.assertEquals(root2over2, cosArgIm, 0, "Expected cos(pi/4) to be sqrt(2) / 2");
595         Assertions.assertEquals(root2over2, sinArgIm, ulp, "Expected sin(pi/4) to be 1 ulp from sqrt(2) / 2");
596         for (final double a : new double[] {0.5, 1.0, 1.2322, 345345.234523}) {
597             final double rootA = Math.sqrt(a);
598             assertComplex(a, 0, name, operation, rootA * cosArgRe, rootA * sinArgRe, 0);
599             // This should be exact. It will fail if using the polar computation
600             // real = sqrt(abs()) * Math.cos(arg() / 2) as cos(pi/2) is not 0.0 but 6.123233995736766e-17
601             assertComplex(-a, 0, name, operation, rootA * sinArgRe, rootA * cosArgRe, 0);
602             // This should be exact. It won't be if Complex is using polar computation
603             // with sin/cos which does not output the same result for angle pi/4.
604             assertComplex(0, a, name, operation, rootA * root2over2, rootA * root2over2, 0);
605             assertComplex(0, -a, name, operation, rootA * root2over2, -rootA * root2over2, 0);
606         }
607 
608         // Check overflow safe.
609         double a = Double.MAX_VALUE;
610         final double b = a / 4;
611         Assertions.assertEquals(inf, Complex.ofCartesian(a, b).abs(), "Expected overflow");
612         // The expected absolute value has been computed using BigDecimal on Java 9
613         //final double newAbs = new BigDecimal(a).multiply(new BigDecimal(a)).add(
614         //                      new BigDecimal(b).multiply(new BigDecimal(b)))
615         //                     .sqrt(MathContext.DECIMAL128).sqrt(MathContext.DECIMAL128).doubleValue()
616         final double newAbs = 1.3612566508088272E154;
617         assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
618                                              newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
619         assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
620                                              newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
621 
622         // Note that the computation is possible in polar coords if abs() does not overflow.
623         a = Double.MAX_VALUE / 2;
624         assertComplex(-a, a, name, operation, 4.3145940638864765e+153, 1.0416351505169177e+154, 2);
625         assertComplex(a, a, name, operation, 1.0416351505169177e+154, 4.3145940638864758e+153);
626         assertComplex(-a, -a, name, operation, 4.3145940638864765e+153,  -1.0416351505169177e+154, 2);
627         assertComplex(a, -a, name, operation, 1.0416351505169177e+154, -4.3145940638864758e+153);
628 
629         // Check minimum normal value conditions
630         // Computing in polar coords produces a very different result with
631         // MIN_VALUE so use MIN_NORMAL
632         a = Double.MIN_NORMAL;
633         assertComplex(-a, a, name, operation, 6.7884304867749663e-155, 1.6388720948399111e-154);
634         assertComplex(a, a, name, operation, 1.6388720948399111e-154, 6.7884304867749655e-155);
635         assertComplex(-a, -a, name, operation, 6.7884304867749663e-155, -1.6388720948399111e-154);
636         assertComplex(a, -a, name, operation, 1.6388720948399111e-154, -6.7884304867749655e-155);
637     }
638 
639     // Note: inf/nan edge cases for
640     // multiply/divide are tested in CStandardTest
641 
642     @Test
643     void testDivide() {
644         final String name = "divide";
645         final BiFunction<Complex, Complex, Complex> operation = Complex::divide;
646 
647         // Should be able to divide by a complex whose absolute (c*c+d*d)
648         // overflows or underflows including all sub-normal numbers.
649 
650         // Worst case is using Double.MIN_VALUE
651         // Should normalise c and d to range [1, 2) resulting in:
652         // c = d = 1
653         // c * c + d * d = 2
654         // scaled x = (a * c + b * d) / denom = Double.MIN_VALUE
655         // scaled y = (b * c - a * d) / denom = 0
656         // The values are rescaled by 1023 + 51 (shift the last bit of the 52 bit mantissa)
657         double x = Math.scalb(Double.MIN_VALUE, 1023 + 51);
658         Assertions.assertEquals(1.0, x);
659         // In other words the result is (x+iy) / (x+iy) = (1+i0)
660         // The result is the same if imaginary is zero (i.e. a real only divide)
661 
662         assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 1.0, 0.0);
663         assertComplex(Double.MAX_VALUE, 0.0, Double.MAX_VALUE, 0.0, name, operation, 1.0, 0.0);
664 
665         assertComplex(1.0, 1.0, 1.0, 1.0, name, operation, 1.0, 0.0);
666         assertComplex(1.0, 0.0, 1.0, 0.0, name, operation, 1.0, 0.0);
667         // Should work for all small values
668         x = Double.MIN_NORMAL;
669         while (x != 0) {
670             assertComplex(x, x, x, x, name, operation, 1.0, 0.0);
671             assertComplex(x, 0, x, 0, name, operation, 1.0, 0.0);
672             x /= 2;
673         }
674 
675         // Some cases of not self-divide
676         assertComplex(1, 1, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, inf, 0);
677         // As computed by GNU g++
678         assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, 4503599627370496L, 0);
679         assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, 2.2204460492503131e-16, 0);
680     }
681 
682     @Test
683     void testPow() {
684         final String name = "pow";
685         final BiFunction<Complex, Complex, Complex> operation = Complex::pow;
686 
687         // pow(Complex) is log().multiply(Complex).exp()
688         // All are overflow safe and handle infinities as defined in the C99 standard.
689         // TODO: Test edge cases with:
690         // Double.MAX_VALUE, Double.MIN_NORMAL, Inf
691         // using other library implementations.
692 
693         // Test NaN
694         assertComplex(1, 1, nan, nan, name, operation, nan, nan);
695         assertComplex(nan, nan, 1, 1, name, operation, nan, nan);
696         assertComplex(nan, 1, 1, 1, name, operation, nan, nan);
697         assertComplex(1, nan, 1, 1, name, operation, nan, nan);
698         assertComplex(1, 1, nan, 1, name, operation, nan, nan);
699         assertComplex(1, 1, 1, nan, name, operation, nan, nan);
700 
701         // Test overflow.
702         assertComplex(Double.MAX_VALUE, 1, 2, 2, name, operation, inf, -inf);
703         assertComplex(1, Double.MAX_VALUE, 2, 2, name, operation, -inf, inf);
704     }
705 }