Complex.java

package org.djutils.complex;

/**
 * Complex.java. Immutable complex numbers. This implementation stores the real and imaginary component as a double value. These
 * fields are directly accessible, or through a getter. There are also getters for the norm and phi, but this are a CPU
 * intensive.
 * <p>
 * Copyright (c) 2021-2021 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
 * BSD-style license. See <a href="https://djutils.org/docs/current/djutils/licenses.html">DJUTILS License</a>.
 * </p>
 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
 */
public class Complex
{
    /** Real component. */
    @SuppressWarnings("checkstyle:visibilitymodifier")
    public final double re;

    /** Imaginary component. */
    @SuppressWarnings("checkstyle:visibilitymodifier")
    public final double im;

    /** The zero value in the complex number space. */
    public static final Complex ZERO = new Complex(0, 0);

    /** The (real) one value in the complex number space. */
    public static final Complex ONE = new Complex(1, 0);

    /** The (real) minus one value in the complex number space. */
    public static final Complex MINUS_ONE = new Complex(-1, 0);

    /** The imaginary unit value (i). */
    public static final Complex I = new Complex(0, 1);

    /** The negative imaginary unit value (i). */
    public static final Complex MINUS_I = new Complex(0, -1);

    /**
     * Construct a new complex number.
     * @param re double; real component of the new complex number
     * @param im double; imaginary component of the new complex number
     */
    public Complex(final double re, final double im)
    {
        this.re = re;
        this.im = im;
    }

    /**
     * Construct a new complex number with specified real component and zero imaginary component.
     * @param re double; the real component of the new complex number
     */
    public Complex(final double re)
    {
        this.re = re;
        this.im = 0;
    }

    /**
     * Retrieve the real component of this complex number.
     * @return double; the real component of this complex number
     */
    public double getRe()
    {
        return this.re;
    }

    /**
     * Retrieve the imaginary component of this complex number.
     * @return double; the imaginary component of this complex number
     */
    public double getIm()
    {
        return this.im;
    }

    /**
     * Compute and return the norm, or radius, or absolute value of this complex number.
     * @return double; the norm, or radius, or absolute value of this complex number. Due to the fact that in this
     *         implementation of complex numbers, all values are stored as a real an imaginary double value; the result may
     *         overflow to Double.POSITIVE_INFINITY, even though both the real and imaginary part of the complex number can be
     *         represented.
     */
    public double norm()
    {
        return hypot(this.re, this.im);
    }

    /** Precision limit. */
    private static final double EPSILONSQRT = Math.sqrt(Math.ulp(1.0) / 2);
    /** Square root of smallest floating point value. */
    private static final double SQRT_OF_MIN_VALUE = Math.sqrt(Double.MIN_VALUE);
    /** Square root of biggest floating point value. */
    private static final double SQRT_OF_MAX_VALUE = Math.sqrt(Double.MAX_VALUE);
    
    /**
     * Better implementation of the hypotenuse function (faster and more accurate than the one in the java Math library). <br>
     * Derived from <a href="https://arxiv.org/abs/1904.09481">An improved algorithm for hypot(a, b) by Carlos F. Borges</a>.
     * @param x double; the x value
     * @param y double; the y value
     * @return double; hypot(x, y)
     */
    public static double hypot(final double x, final double y)
    {
        if (x != x || y != y)
        {
            return Double.NaN;
        }
        double absX = Math.abs(x);
        double absY = Math.abs(y);
        if (absX == Double.POSITIVE_INFINITY || absY == Double.POSITIVE_INFINITY)
        {
            return Double.POSITIVE_INFINITY;
        }
        if (absX < absY)
        {
            double swap = absX;
            absX = absY;
            absY = swap;
        }
        if (absY <= absX * EPSILONSQRT)
        {
            return absX;
        }
        double scale = SQRT_OF_MIN_VALUE;
        if (absX > SQRT_OF_MAX_VALUE)
        {
            absX *= scale;
            absY *= scale;
            scale = 1.0 / scale;
        }
        else if (absY < SQRT_OF_MIN_VALUE)
        {
            absX /= scale;
            absY /= scale;
        }
        else
        {
            scale = 1.0;
        }
        double h = Math.sqrt(Math.fma(absX, absX, absY * absY));
        double hsq = h * h;
        double xsq = absX * absX;
        double a = Math.fma(-absY, absY, hsq - xsq) + Math.fma(h, h, -hsq) - Math.fma(absX, absX, -xsq);
        return scale * (h - a / (2.0 * h));
    }

    /**
     * Compute and return the phase or phi of this complex number.
     * @return double; the phase or phi of this complex number in Radians. Due to the fact that in this implementation of
     *         complex numbers, all values are stored as a real and imaginary value; the result of this method is always
     *         normalized to the interval (-&pi;,&pi;].
     */
    public double phi()
    {
        return Math.atan2(this.im, this.re);
    }

    /**
     * Determine if this Complex has an imaginary component of zero.
     * @return boolean; true if the imaginary component of this Complex number is 0.0; false if the imaginary component of this
     *         Complex number is not 0.0.
     */
    public boolean isReal()
    {
        return this.im == 0.0;
    }

    /**
     * Determine if this Complex has a real component of zero.
     * @return boolean; true if the real component of this Complex number is 0.0; false if the real component of this Complex
     *         number is not 0.0
     */
    public boolean isImaginary()
    {
        return this.re == 0.0;
    }

    /**
     * Construct the complex conjugate of this Complex.
     * @return Complex; the complex conjugate of this
     */
    public Complex conjugate()
    {
        return new Complex(this.re, -this.im);
    }

    /**
     * Rotate this Complex by an angle.
     * @param angle double; the angle (in Radians)
     * @return complex; the result of the rotation
     */
    public Complex rotate(final double angle)
    {
        double sin = Math.sin(angle);
        double cos = Math.cos(angle);
        return new Complex(this.re * cos - this.im * sin, this.im * cos + this.re * sin);
    }

    /**
     * Add this Complex and another Complex.
     * @param rightOperand Complex; the other Complex
     * @return Complex; the sum of this Complex and the other Complex
     */
    public Complex plus(final Complex rightOperand)
    {
        return new Complex(this.re + rightOperand.re, this.im + rightOperand.im);
    }

    /**
     * Add a scalar to this Complex.
     * @param rightOperand double; the scalar
     * @return Complex; the sum of this Complex and the scalar
     */
    public Complex plus(final double rightOperand)
    {
        return new Complex(this.re + rightOperand, this.im);
    }

    /**
     * Subtract another Complex from this Complex.
     * @param rightOperand Complex; the other Complex
     * @return Complex; the difference of this Complex and the other Complex
     */
    public Complex minus(final Complex rightOperand)
    {
        return new Complex(this.re - rightOperand.re, this.im - rightOperand.im);
    }

    /**
     * Subtract a scalar from this Complex.
     * @param rightOperand double; the scalar
     * @return Complex; the difference of this Complex and the scalar
     */
    public Complex minus(final double rightOperand)
    {
        return new Complex(this.re - rightOperand, this.im);
    }

    /**
     * Multiply this Complex with another Complex.
     * @param rightOperand Complex; the right hand side operand
     * @return Complex; the product of this Complex and the other Complex
     */
    public Complex times(final Complex rightOperand)
    {
        return new Complex(this.re * rightOperand.re - this.im * rightOperand.im,
                this.im * rightOperand.re + this.re * rightOperand.im);
    }

    /**
     * Multiply this Complex with a scalar.
     * @param rightOperand double; the right hand side operand
     * @return Complex; the product of this Complex and the scalar
     */
    public Complex times(final double rightOperand)
    {
        return new Complex(this.re * rightOperand, this.im * rightOperand);
    }

    /**
     * Compute the reciprocal of this Complex. If this is zero, the result will have the re field set to
     * Double.POSITIVE_INFINITY and the imaginary part set to NEGATIVE_INFINITY.
     * @return Complex; the reciprocal of this Complex (1 / this)
     */
    public Complex reciprocal()
    {
        double divisor = this.re * this.re + this.im * this.im;
        if (0.0 == divisor)
        {
            return new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
        }
        return new Complex(this.re / divisor, -this.im / divisor);
    }

    /**
     * Divide this Complex by another Complex. Division by ZERO yields a Complex with re set to Infinity if this.re != 0 and NaN
     * if this.re == 0 and im set to Infinity if this.im != 0 and NaN if this.im == 0.
     * @param rightOperand Complex; the right hand side operand
     * @return Complex; the ratio of this Complex and the right hand side operand
     */
    public Complex divideBy(final Complex rightOperand)
    {
        if (rightOperand.re == 0 && rightOperand.im == 0)
        {
            return new Complex(this.re / 0.0, this.im / 0.0);
        }
        return this.times(rightOperand.reciprocal());
    }

    /**
     * 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
     * this.re == 0 and im set to Infinity if this.im != 0 and NaN if this.im == 0.
     * @param rightOperand double; the scalar right hand side operand
     * @return Complex; the ratio of this Complex and the right hand side operand
     */
    public Complex divideBy(final double rightOperand)
    {
        return new Complex(this.re / rightOperand, this.im / rightOperand);
    }

    /** {@inheritDoc} */
    @Override
    public String toString()
    {
        return "Complex [re=" + this.re + ", im=" + this.im + "]";
    }

    /** {@inheritDoc} */
    @Override
    public int hashCode()
    {
        final int prime = 31;
        int result = 1;
        long temp;
        temp = Double.doubleToLongBits(this.im);
        result = prime * result + (int) (temp ^ (temp >>> 32));
        temp = Double.doubleToLongBits(this.re);
        result = prime * result + (int) (temp ^ (temp >>> 32));
        return result;
    }

    /** {@inheritDoc} */
    @SuppressWarnings("checkstyle:needbraces")
    @Override
    public boolean equals(final Object obj)
    {
        if (this == obj)
            return true;
        if (obj == null)
            return false;
        if (getClass() != obj.getClass())
            return false;
        Complex other = (Complex) obj;
        if (Double.doubleToLongBits(this.im) != Double.doubleToLongBits(other.im))
            return false;
        if (Double.doubleToLongBits(this.re) != Double.doubleToLongBits(other.re))
            return false;
        return true;
    }

}