1 package org.djutils.complex;
2
3 /**
4 * Complex.java. Immutable complex numbers. This implementation stores the real and imaginary component as a double value. These
5 * fields are directly accessible, or through a getter. There are also getters for the norm and phi, but this are a CPU
6 * intensive.
7 * <p>
8 * Copyright (c) 2021-2021 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
9 * BSD-style license. See <a href="https://djutils.org/docs/current/djutils/licenses.html">DJUTILS License</a>.
10 * </p>
11 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
12 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
13 */
14 public class Complex
15 {
16 /** Real component. */
17 @SuppressWarnings("checkstyle:visibilitymodifier")
18 public final double re;
19
20 /** Imaginary component. */
21 @SuppressWarnings("checkstyle:visibilitymodifier")
22 public final double im;
23
24 /** The zero value in the complex number space. */
25 public static final Complexml#Complex">Complex ZERO = new Complex(0, 0);
26
27 /** The (real) one value in the complex number space. */
28 public static final Complextml#Complex">Complex ONE = new Complex(1, 0);
29
30 /** The (real) minus one value in the complex number space. */
31 public static final Complexmplex">Complex MINUS_ONE = new Complex(-1, 0);
32
33 /** The imaginary unit value (i). */
34 public static final Complex.html#Complex">Complex I = new Complex(0, 1);
35
36 /** The negative imaginary unit value (i). */
37 public static final ComplexComplex">Complex MINUS_I = new Complex(0, -1);
38
39 /**
40 * Construct a new complex number.
41 * @param re double; real component of the new complex number
42 * @param im double; imaginary component of the new complex number
43 */
44 public Complex(final double re, final double im)
45 {
46 this.re = re;
47 this.im = im;
48 }
49
50 /**
51 * Construct a new complex number with specified real component and zero imaginary component.
52 * @param re double; the real component of the new complex number
53 */
54 public Complex(final double re)
55 {
56 this.re = re;
57 this.im = 0;
58 }
59
60 /**
61 * Retrieve the real component of this complex number.
62 * @return double; the real component of this complex number
63 */
64 public double getRe()
65 {
66 return this.re;
67 }
68
69 /**
70 * Retrieve the imaginary component of this complex number.
71 * @return double; the imaginary component of this complex number
72 */
73 public double getIm()
74 {
75 return this.im;
76 }
77
78 /**
79 * Compute and return the norm, or radius, or absolute value of this complex number.
80 * @return double; the norm, or radius, or absolute value of this complex number. Due to the fact that in this
81 * implementation of complex numbers, all values are stored as a real an imaginary double value; the result may
82 * overflow to Double.POSITIVE_INFINITY, even though both the real and imaginary part of the complex number can be
83 * represented.
84 */
85 public double norm()
86 {
87 return hypot(this.re, this.im);
88 }
89
90 /** Precision limit. */
91 private static final double EPSILONSQRT = Math.sqrt(Math.ulp(1.0) / 2);
92 /** Square root of smallest floating point value. */
93 private static final double SQRT_OF_MIN_VALUE = Math.sqrt(Double.MIN_VALUE);
94 /** Square root of biggest floating point value. */
95 private static final double SQRT_OF_MAX_VALUE = Math.sqrt(Double.MAX_VALUE);
96
97 /**
98 * Better implementation of the hypotenuse function (faster and more accurate than the one in the java Math library). <br>
99 * Derived from <a href="https://arxiv.org/abs/1904.09481">An improved algorithm for hypot(a, b) by Carlos F. Borges</a>.
100 * @param x double; the x value
101 * @param y double; the y value
102 * @return double; hypot(x, y)
103 */
104 public static double hypot(final double x, final double y)
105 {
106 if (x != x || y != y)
107 {
108 return Double.NaN;
109 }
110 double absX = Math.abs(x);
111 double absY = Math.abs(y);
112 if (absX == Double.POSITIVE_INFINITY || absY == Double.POSITIVE_INFINITY)
113 {
114 return Double.POSITIVE_INFINITY;
115 }
116 if (absX < absY)
117 {
118 double swap = absX;
119 absX = absY;
120 absY = swap;
121 }
122 if (absY <= absX * EPSILONSQRT)
123 {
124 return absX;
125 }
126 double scale = SQRT_OF_MIN_VALUE;
127 if (absX > SQRT_OF_MAX_VALUE)
128 {
129 absX *= scale;
130 absY *= scale;
131 scale = 1.0 / scale;
132 }
133 else if (absY < SQRT_OF_MIN_VALUE)
134 {
135 absX /= scale;
136 absY /= scale;
137 }
138 else
139 {
140 scale = 1.0;
141 }
142 double h = Math.sqrt(Math.fma(absX, absX, absY * absY));
143 double hsq = h * h;
144 double xsq = absX * absX;
145 double a = Math.fma(-absY, absY, hsq - xsq) + Math.fma(h, h, -hsq) - Math.fma(absX, absX, -xsq);
146 return scale * (h - a / (2.0 * h));
147 }
148
149 /**
150 * Compute and return the phase or phi of this complex number.
151 * @return double; the phase or phi of this complex number in Radians. Due to the fact that in this implementation of
152 * complex numbers, all values are stored as a real and imaginary value; the result of this method is always
153 * normalized to the interval (-π,π].
154 */
155 public double phi()
156 {
157 return Math.atan2(this.im, this.re);
158 }
159
160 /**
161 * Determine if this Complex has an imaginary component of zero.
162 * @return boolean; true if the imaginary component of this Complex number is 0.0; false if the imaginary component of this
163 * Complex number is not 0.0.
164 */
165 public boolean isReal()
166 {
167 return this.im == 0.0;
168 }
169
170 /**
171 * Determine if this Complex has a real component of zero.
172 * @return boolean; true if the real component of this Complex number is 0.0; false if the real component of this Complex
173 * number is not 0.0
174 */
175 public boolean isImaginary()
176 {
177 return this.re == 0.0;
178 }
179
180 /**
181 * Construct the complex conjugate of this Complex.
182 * @return Complex; the complex conjugate of this
183 */
184 public Complex conjugate()
185 {
186 return new Complex(this.re, -this.im);
187 }
188
189 /**
190 * Rotate this Complex by an angle.
191 * @param angle double; the angle (in Radians)
192 * @return complex; the result of the rotation
193 */
194 public Complex rotate(final double angle)
195 {
196 double sin = Math.sin(angle);
197 double cos = Math.cos(angle);
198 return new Complex(this.re * cos - this.im * sin, this.im * cos + this.re * sin);
199 }
200
201 /**
202 * Add this Complex and another Complex.
203 * @param rightOperand Complex; the other Complex
204 * @return Complex; the sum of this Complex and the other Complex
205 */
206 public Complexml#Complex">Complex plus(final Complex rightOperand)
207 {
208 return new Complex(this.re + rightOperand.re, this.im + rightOperand.im);
209 }
210
211 /**
212 * Add a scalar to this Complex.
213 * @param rightOperand double; the scalar
214 * @return Complex; the sum of this Complex and the scalar
215 */
216 public Complex plus(final double rightOperand)
217 {
218 return new Complex(this.re + rightOperand, this.im);
219 }
220
221 /**
222 * Subtract another Complex from this Complex.
223 * @param rightOperand Complex; the other Complex
224 * @return Complex; the difference of this Complex and the other Complex
225 */
226 public Complexl#Complex">Complex minus(final Complex rightOperand)
227 {
228 return new Complex(this.re - rightOperand.re, this.im - rightOperand.im);
229 }
230
231 /**
232 * Subtract a scalar from this Complex.
233 * @param rightOperand double; the scalar
234 * @return Complex; the difference of this Complex and the scalar
235 */
236 public Complex minus(final double rightOperand)
237 {
238 return new Complex(this.re - rightOperand, this.im);
239 }
240
241 /**
242 * Multiply this Complex with another Complex.
243 * @param rightOperand Complex; the right hand side operand
244 * @return Complex; the product of this Complex and the other Complex
245 */
246 public Complexl#Complex">Complex times(final Complex rightOperand)
247 {
248 return new Complex(this.re * rightOperand.re - this.im * rightOperand.im,
249 this.im * rightOperand.re + this.re * rightOperand.im);
250 }
251
252 /**
253 * Multiply this Complex with a scalar.
254 * @param rightOperand double; the right hand side operand
255 * @return Complex; the product of this Complex and the scalar
256 */
257 public Complex times(final double rightOperand)
258 {
259 return new Complex(this.re * rightOperand, this.im * rightOperand);
260 }
261
262 /**
263 * Compute the reciprocal of this Complex. If this is zero, the result will have the re field set to
264 * Double.POSITIVE_INFINITY and the imaginary part set to NEGATIVE_INFINITY.
265 * @return Complex; the reciprocal of this Complex (1 / this)
266 */
267 public Complex reciprocal()
268 {
269 double divisor = this.re * this.re + this.im * this.im;
270 if (0.0 == divisor)
271 {
272 return new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
273 }
274 return new Complex(this.re / divisor, -this.im / divisor);
275 }
276
277 /**
278 * Divide this Complex by another Complex. Division by ZERO yields a Complex with re set to Infinity if this.re != 0 and NaN
279 * if this.re == 0 and im set to Infinity if this.im != 0 and NaN if this.im == 0.
280 * @param rightOperand Complex; the right hand side operand
281 * @return Complex; the ratio of this Complex and the right hand side operand
282 */
283 public Complexomplex">Complex divideBy(final Complex rightOperand)
284 {
285 if (rightOperand.re == 0 && rightOperand.im == 0)
286 {
287 return new Complex(this.re / 0.0, this.im / 0.0);
288 }
289 return this.times(rightOperand.reciprocal());
290 }
291
292 /**
293 * Divide this Complex by a scalar. Division by 0.0 yields a Complex with re set to Infinity if this.re != 0 and NaN if
294 * this.re == 0 and im set to Infinity if this.im != 0 and NaN if this.im == 0.
295 * @param rightOperand double; the scalar right hand side operand
296 * @return Complex; the ratio of this Complex and the right hand side operand
297 */
298 public Complex divideBy(final double rightOperand)
299 {
300 return new Complex(this.re / rightOperand, this.im / rightOperand);
301 }
302
303 /** {@inheritDoc} */
304 @Override
305 public String toString()
306 {
307 return "Complex [re=" + this.re + ", im=" + this.im + "]";
308 }
309
310 /** {@inheritDoc} */
311 @Override
312 public int hashCode()
313 {
314 final int prime = 31;
315 int result = 1;
316 long temp;
317 temp = Double.doubleToLongBits(this.im);
318 result = prime * result + (int) (temp ^ (temp >>> 32));
319 temp = Double.doubleToLongBits(this.re);
320 result = prime * result + (int) (temp ^ (temp >>> 32));
321 return result;
322 }
323
324 /** {@inheritDoc} */
325 @SuppressWarnings("checkstyle:needbraces")
326 @Override
327 public boolean equals(final Object obj)
328 {
329 if (this == obj)
330 return true;
331 if (obj == null)
332 return false;
333 if (getClass() != obj.getClass())
334 return false;
335 Complex../../../org/djutils/complex/Complex.html#Complex">Complex other = (Complex) obj;
336 if (Double.doubleToLongBits(this.im) != Double.doubleToLongBits(other.im))
337 return false;
338 if (Double.doubleToLongBits(this.re) != Double.doubleToLongBits(other.re))
339 return false;
340 return true;
341 }
342
343 }