plus,
* minus, times or over.
* * public static VNComplex op(VNComplex x, VNComplex y) // x op y * public static VNComplex op(VNComplex x, double y) // x op y * public static VNComplex op(double x, VNComplex y) // x op y * public VNComplex op(VNComplex y) // this op y * public VNComplex op(double y) // this op y * public VNComplex opReverse(double x) // x op this ** * The functions in this class follow the rules for VNComplex arithmetic * as defined C9x Annex G:"IEC 559-compatible VNComplex arithmetic." * The API is not the same, but handling of infinities, NaNs, and positive * and negative zeros is intended to follow the same rules. * * This class depends on the standard java.lang.Math class following * certain rules, as defined in the C9x Annex F, for the handling of * infinities, NaNs, and positive and negative zeros. Sun's specification * is that java.lang.Math should reproduce the results in the Sun's fdlibm * C library. This library appears to follow the Annex F specification. * At least on Windows, Sun's JDK 1.0 and 1.1 do NOT follow this specification. * Sun's JDK 1.2(RC2) does follow the Annex F specification. Thesefore, * this class will not give the expected results for edge cases with * JDK 1.0 and 1.1. */ public class VNComplex implements java.io.Serializable, Cloneable { /** * @serial Real part of the VNComplex. */ private double re; /** * @serial Imaginary part of the VNComplex. */ private double im; /** * Serialization ID */ static final long serialVersionUID = -633126172485117692L; /** * String used in converting VNComplex to String. * Default is "i", but sometimes "j" is desired. * Note that this is set for the class, not for * a particular instance of a VNComplex. */ public static String suffix = "i"; private final static long negZeroBits = Double.doubleToLongBits(1.0/Double.NEGATIVE_INFINITY); /** * Constructs a VNComplex equal to zero. */ public VNComplex() { re = 0.0; im = 0.0; } /** * Constructs a VNComplex with a zero imaginary part. * @param re A double value equal to the real part of the VNComplex object. */ public VNComplex(double re) { this.re = re; this.im = 0.0; } /** * Constructs a VNComplex with real and imaginary parts given * by the input arguments. * @param re A double value equal to the real part of the VNComplex object. * @param im A double value equal to the imaginary part of the VNComplex object. */ public VNComplex(double re, double im) { this.re = re; this.im = im; } /** * Constructs a VNComplex equal to the argument. * @param z A VNComplex object * If z is null then a NullPointerException is thrown. */ public VNComplex(VNComplex z) { re = z.re; im = z.im; } /** * Returns the absolute value (modulus) of a VNComplex, |z|. * @param z A VNComplex object. * @return A double value equal to the absolute value of the argument. */ public static double abs(VNComplex z) { double x = Math.abs(z.re); double y = Math.abs(z.im); if (Double.isInfinite(x) || Double.isInfinite(y)) return Double.POSITIVE_INFINITY; if (x + y == 0.0) { return 0.0; } else if (x > y) { y /= x; return x*Math.sqrt(1.0+y*y); } else { x /= y; return y*Math.sqrt(x*x+1.0); } } /** * Returns the inverse cosine (arc cosine) of a VNComplex, * with branch cuts outside the interval [-1,1] along the * real axis. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * inverse (arc) cosine of the argument. * The real part of the result is in the interval [0,pi]. */ public static VNComplex acos(VNComplex z) { VNComplex result = new VNComplex(); double r = abs(z); if (Double.isInfinite(z.re) && Double.isNaN(z.im)) { result.re = Double.NaN; result.im = Double.NEGATIVE_INFINITY; } else if (Double.isInfinite(r)) { result.re = Math.atan2(Math.abs(z.im),z.re); result.im = z.im*Double.NEGATIVE_INFINITY; } else if (r == 0) { result.re = Math.PI/2; result.im = -z.im; } else { result = minus(Math.PI/2,asin(z)); } return result; } /** * Returns the inverse hyperbolic cosine (arc cosh) of a VNComplex, * with a branch cut at values less than one along the real axis. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * inverse (arc) hyperbolic cosine of the argument. * The real part of the result is non-negative and its * imaginary part is in the interval [-i*pi,i*pi]. */ public static VNComplex acosh(VNComplex z) { VNComplex result = acos(z); double rx = -result.im; result.im = result.re; result.re = rx; if (result.re < 0 || isNegZero(result.re)) { result.re = -result.re; result.im = -result.im; } return result; } /** * Returns the argument (phase) of a VNComplex, in radians, * with a branch cut along the negative real axis. * @param z A VNComplex object. * @return A double value equal to the argument (or phase) of a VNComplex. * It is in the interval [-pi,pi]. */ public static double argument(VNComplex z) { return Math.atan2(z.im, z.re); } /** * Returns the inverse sine (arc sine) of a VNComplex, * with branch cuts outside the interval [-1,1] along the * real axis. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to inverse * (arc) sine of the argument. The real part of the * result is in the interval [-pi/2,+pi/2]. */ public static VNComplex asin(VNComplex z) { VNComplex result = new VNComplex(); double r = abs(z); if (Double.isInfinite(r)) { boolean infiniteX = Double.isInfinite(z.re); boolean infiniteY = Double.isInfinite(z.im); if (infiniteX) { double pi2 = 0.5*Math.PI; result.re = (z.re>0 ? pi2 : -pi2); if (infiniteY) result.re /= 2; } else if (infiniteY) { result.re = z.re/Double.POSITIVE_INFINITY; } if (Double.isNaN(z.im)) { result.im = -z.re; result.re = z.im; } else { result.im = z.im*Double.POSITIVE_INFINITY; } return result; } else if (Double.isNaN(r)) { result.re = result.im = Double.NaN; if (z.re == 0) result.re = z.re; } else if (r < 2.58095e-08) { // sqrt(6.0*dmach(3)) = 2.58095e-08 result.re = z.re; result.im = z.im; } else if (z.re == 0) { result.re = 0; result.im = Sfun.asinh(z.im); } else if (r <= 0.1) { VNComplex z2 = times(z,z); //log(eps)/log(rmax) = 8 where rmax = 0.1 for (int i = 1; i <= 8; i++) { double twoi = 2*(8-i) + 1; result = times(times(result,z2),twoi/(twoi+1.0)); result.re += 1.0/twoi; } result = result.times(z); } else { // A&S 4.4.26 // asin(z) = -i*log(z+sqrt(1-z)*sqrt(1+z)) // or, since log(iz) = log(z) +i*pi/2, // asin(z) = pi/2 - i*log(z+sqrt(z+1)*sqrt(z-1)) VNComplex w = ((z.im < 0) ? negative(z) : z); VNComplex sqzp1 = sqrt(plus(w,1.0)); if (sqzp1.im < 0.0) sqzp1 = negative(sqzp1); VNComplex sqzm1 = sqrt(minus(w,1.0)); result = log(plus(w,times(sqzp1,sqzm1))); double rx = result.re; result.re = 0.5*Math.PI + result.im; result.im = -rx; } if (result.re > 0.5*Math.PI) { result.re = Math.PI - result.re; result.im = -result.im; } if (result.re < -0.5*Math.PI) { result.re = -Math.PI - result.re; result.im = -result.im; } if (z.im < 0) { result.re = -result.re; result.im = -result.im; } return result; } /** * Returns the inverse hyperbolic sine (arc sinh) of a VNComplex, * with a branch cuts outside the interval [-i,i]. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * inverse (arc) hyperbolic sine of the argument. * Its imaginary part is in the interval [-i*pi/2,i*pi/2]. */ public static VNComplex asinh(VNComplex z) { // asinh(z) = i*asin(-i*z) VNComplex miz = new VNComplex(z.im,-z.re); VNComplex result = asin(miz); double rx = result.im; result.im = result.re; result.re = -rx; return result; } /** * Returns the inverse tangent (arc tangent) of a VNComplex, * with branch cuts outside the interval [-i,i] along the * imaginary axis. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * inverse (arc) tangent of the argument. * Its real part is in the interval [-pi/2,pi/2]. */ public static VNComplex atan(VNComplex z) { VNComplex result = new VNComplex(); double r = abs(z); if (Double.isInfinite(r)) { double pi2 = 0.5*Math.PI; double im = (Double.isNaN(z.im) ? 0 : z.im); result.re = (z.re<0 ? -pi2 : pi2); result.im = (im<0 ? -1 : 1)/Double.POSITIVE_INFINITY; if (Double.isNaN(z.re)) result.re = z.re; } else if (Double.isNaN(r)) { result.re = result.im = Double.NaN; if (z.im == 0) result.im = z.im; } else if (r < 1.82501e-08) { // sqrt(3.0*dmach(3)) = 1.82501e-08 result.re = z.re; result.im = z.im; } else if (r < 0.1) { VNComplex z2 = times(z,z); // -0.4343*log(dmach(3))+1 = 17 for (int k = 0; k < 17; k++) { VNComplex temp = times(z2,result); int twoi = 2*(17-k) - 1; result.re = 1.0/twoi - temp.re; result.im = -temp.im; } result = result.times(z); } else if (r < 9.0072e+15) { // 1.0/dmach(3) = 9.0072e+15 double r2 = r*r; result.re = 0.5*Math.atan2(2*z.re,1.0-r2); result.im = 0.25*Math.log((r2+2*z.im+1)/(r2-2*z.im+1)); } else { result.re = ((z.re < 0.0) ? -0.5*Math.PI : 0.5*Math.PI); } return result; } /** * Returns the inverse hyperbolic tangent (arc tanh) of a VNComplex, * with a branch cuts outside the interval [-1,1] on the real axis. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * inverse (arc) hyperbolic tangent of the argument. * The imaginary part of the result is in the interval * [-i*pi/2,i*pi/2]. */ public static VNComplex atanh(VNComplex z) { // atanh(z) = i*atan(-i*z) VNComplex miz = new VNComplex(z.im,-z.re); VNComplex result = atan(miz); double rx = result.im; result.im = result.re; result.re = -rx; return result; } /** * Returns the VNComplex conjugate of a VNComplex object. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to VNComplex conjugate of z. */ public static VNComplex conjugate(VNComplex z) { return new VNComplex(z.re, -z.im); } /* * Returns sign(b)*|a|. */ private static double copysign(double a, double b) { double abs = Math.abs(a); return ((b < 0) ? -abs : abs); } /** * Returns the cosine of a VNComplex. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to cosine of the argument. */ public static VNComplex cos(VNComplex z) { // cos(z) = cosh(i*z) return cosh(new VNComplex(-z.im,z.re)); } /** * Returns the hyperbolic cosh of a VNComplex. * @param z A VNComplex object. * @return A newly constructed VNComplex initialized to * the hyperbolic cosine of the argument. */ public static VNComplex cosh(VNComplex z) { if (z.im == 0) { return new VNComplex(Sfun.cosh(z.re)); } double coshx = Sfun.cosh(z.re); double sinhx = Sfun.sinh(z.re); double cosy = Math.cos(z.im); double siny = Math.sin(z.im); boolean infiniteX = Double.isInfinite(coshx); boolean infiniteY = Double.isInfinite(z.im); // A&S 4.5.50 VNComplex result = new VNComplex(coshx*cosy, sinhx*siny); if (infiniteY) result.re = Double.NaN; if (z.re == 0) { result.im = 0; } else if (infiniteX) { result.re = z.re*cosy; result.im = z.re*siny; if (z.im == 0) result.im = 0; if (Double.isNaN(z.im)) { result.re = z.re; } else if (infiniteY) { result.re = z.im; } } return result; } /** * Compares with another VNComplex. *
Note: To be useful in hashtables this method * considers two NaN double values to be equal. This * is not according to IEEE specification. * @param z A VNComplex object. * @return True if the real and imaginary parts of this object * are equal to their counterparts in the argument; false, otherwise. */ public boolean equals(VNComplex z) { if (isNaN() && z.isNaN()) { return true; } else { return (re == z.re && im == z.im); } } /** * Compares this object against the specified object. *
Note: To be useful in hashtables this method
* considers two NaN double values to be equal. This
* is not according to IEEE specification
* @param obj The object to compare with.
* @return True if the objects are the same; false otherwise.
*/
public boolean equals(Object obj)
{
if (obj == null) {
return false;
} else if (obj instanceof VNComplex) {
return equals((VNComplex)obj);
} else {
return false;
}
}
/**
* Returns the exponential of a VNComplex z, exp(z).
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to exponential
* of the argument.
*/
public static VNComplex exp(VNComplex z)
{
VNComplex result = new VNComplex();
double r = Math.exp(z.re);
double cosa = Math.cos(z.im);
double sina = Math.sin(z.im);
if (Double.isInfinite(z.im) || Double.isNaN(z.im) || Math.abs(cosa)>1) {
cosa = sina = Double.NaN;
}
if (Double.isInfinite(z.re) || Double.isInfinite(r)) {
if (z.re < 0) {
r = 0;
if (Double.isInfinite(z.im) || Double.isNaN(z.im)) {
cosa = sina = 0;
} else {
cosa /= Double.POSITIVE_INFINITY;
sina /= Double.POSITIVE_INFINITY;
}
} else {
r = z.re;
if (Double.isNaN(z.im)) cosa = 1;
}
}
if (z.im == 0.0) {
result.re = r;
result.im = z.im;
} else {
result.re = r*cosa;
result.im = r*sina;
}
return result;
}
/**
* Returns a hashcode for this VNComplex.
* @return A hash code value for this object.
*/
public int hashCode()
{
long re_bits = Double.doubleToLongBits(re);
long im_bits = Double.doubleToLongBits(im);
return (int)((re_bits^im_bits)^((re_bits^im_bits)>>32));
}
/**
* Returns the imaginary part of a VNComplex object.
* @param z A VNComplex object.
* @return The imaginary part of z.
*/
public double imag()
{
return im;
}
/**
* Returns the imaginary part of a VNComplex object.
* @param z A VNComplex object.
* @return The imaginary part of z.
*/
public static double imag(VNComplex z)
{
return z.im;
}
private static boolean isFinite(double x)
{
return !(Double.isInfinite(x) || Double.isNaN(x));
}
/**
* Tests if this is a VNComplex Not-a-Number (NaN) value.
* @return True if either component of the VNComplex object is NaN;
* false, otherwise.
*/
private boolean isNaN()
{
return (Double.isNaN(re) || Double.isNaN(im));
}
/**
* Returns true is x is a negative zero.
*/
private static boolean isNegZero(double x)
{
return (Double.doubleToLongBits(x) == negZeroBits);
}
/**
* Returns the logarithm of a VNComplex z,
* with a branch cut along the negative real axis.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to logarithm
* of the argument. Its imaginary part is in the
* interval [-i*pi,i*pi].
*/
public static VNComplex log(VNComplex z)
{
VNComplex result = new VNComplex();
if (Double.isNaN(z.re)) {
result.re = result.im = z.re;
if (Double.isInfinite(z.im))
result.re = Double.POSITIVE_INFINITY;
} else if (Double.isNaN(z.im)) {
result.re = result.im = z.im;
if (Double.isInfinite(z.re))
result.re = Double.POSITIVE_INFINITY;
} else {
result.re = Math.log(abs(z));
result.im = argument(z);
}
return result;
}
/**
* Subtracts a double from this VNComplex and returns the difference, this-y.
* @param y A double value.
* @return A newly constructed VNComplex initialized to this-y.
*/
public VNComplex minus(double y)
{
return new VNComplex(re-y, im);
}
/**
* Returns the difference of a double and a VNComplex object, x-y.
* @param x A double value.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x-y..
*/
public static VNComplex minus(double x, VNComplex y)
{
return new VNComplex(x-y.re, -y.im);
}
/**
* Returns the difference of this VNComplex object and
* another VNComplex object, this-y.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to this-y.
*/
public VNComplex minus(VNComplex y)
{
return new VNComplex(re-y.re, im-y.im);
}
/**
* Returns the difference of a VNComplex object and a double, x-y.
* @param x A VNComplex object.
* @param y A double value.
* @return A newly constructed VNComplex initialized to x-y.
*/
public static VNComplex minus(VNComplex x, double y)
{
return new VNComplex(x.re-y, x.im);
}
/**
* Returns the difference of two VNComplex objects, x-y.
* @param x A VNComplex object.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x-y.
*/
public static VNComplex minus(VNComplex x, VNComplex y)
{
return new VNComplex(x.re-y.re, x.im-y.im);
}
/**
* Returns the difference of this VNComplex object and a double, this-y.
* @param y A double value.
* @return A newly constructed VNComplex initialized to x-this.
*/
public VNComplex minusReverse(double x)
{
return new VNComplex(x-re, -im);
}
/**
* Returns the negative of a VNComplex object, -z.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to
* the negative of the argument.
*/
public static VNComplex negative(VNComplex z)
{
return new VNComplex(-z.re, -z.im);
}
/**
* Returns this VNComplex object divided by double, this/y.
* @param y The denominator, a double.
* @return A newly constructed VNComplex initialized to x/y.
*/
public VNComplex over(double y)
{
return over(this, y);
}
/**
* Returns a double divided by a VNComplex object, x/y.
* @param x A double value.
* @param y The denominator, a VNComplex object.
* @return A newly constructed VNComplex initialized to x/y.
*/
public static VNComplex over(double x, VNComplex y)
{
return y.overReverse(x);
}
/**
* Returns this VNComplex object divided by another VNComplex object, this/y.
* @param y The denominator, a VNComplex object.
* @return A newly constructed VNComplex initialized to x/y.
*/
public VNComplex over(VNComplex y)
{
return over(this, y);
}
/**
* Returns VNComplex object divided by a double, x/y.
* @param x The numerator, a VNComplex object.
* @param y The denominator, a double.
* @return A newly constructed VNComplex initialized to x/y.
*/
public static VNComplex over(VNComplex x, double y)
{
return new VNComplex(x.re/y, x.im/y);
}
/**
* Returns VNComplex object divided by a VNComplex object, x/y.
* @param x The numerator, a VNComplex object.
* @param y The denominator, a VNComplex object.
* @return A newly constructed VNComplex initialized to x/y.
*/
public static VNComplex over(VNComplex x, VNComplex y)
{
double a = x.re;
double b = x.im;
double c = y.re;
double d = y.im;
double scale = Math.max(Math.abs(c), Math.abs(d));
boolean isScaleFinite = isFinite(scale);
if (isScaleFinite) {
c /= scale;
d /= scale;
}
double den = c*c + d*d;
VNComplex z = new VNComplex((a*c+b*d)/den, (b*c-a*d)/den);
if (isScaleFinite) {
z.re /= scale;
z.im /= scale;
}
// Recover infinities and zeros computed as NaN+iNaN.
if (Double.isNaN(z.re) && Double.isNaN(z.im)) {
if (den == 0.0 && (!Double.isNaN(a) || !Double.isNaN(b))) {
double s = copysign(Double.POSITIVE_INFINITY, c);
z.re = s * a;
z.im = s * b;
} else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
isFinite(c) && isFinite(d)) {
a = copysign(Double.isInfinite(a)?1.0:0.0, a);
b = copysign(Double.isInfinite(b)?1.0:0.0, b);
z.re = Double.POSITIVE_INFINITY * (a*c + b*d);
z.im = Double.POSITIVE_INFINITY * (b*c - a*d);
} else if (Double.isInfinite(scale) &&
isFinite(a) && isFinite(b)) {
c = copysign(Double.isInfinite(c)?1.0:0.0, c);
d = copysign(Double.isInfinite(d)?1.0:0.0, d);
z.re = 0.0 * (a*c + b*d);
z.im = 0.0 * (b*c - a*d);
}
}
return z;
}
/**
* Returns a double dividied by this VNComplex object, x/this.
* @param x The numerator, a double.
* @return A newly constructed VNComplex initialized to x/this.
*/
public VNComplex overReverse(double x)
{
double den, t;
VNComplex z;
if (Math.abs(re) > Math.abs(im)) {
t = im / re;
den = re + im*t;
z = new VNComplex(x/den, -x*t/den);
} else {
t = re / im;
den = im + re*t;
z = new VNComplex(x*t/den, -x/den);
}
return z;
}
/**
* Returns the sum of this VNComplex a double, this+y.
* @param y A double value.
* @return A newly constructed VNComplex initialized to this+y.
*/
public VNComplex plus(double y)
{
return new VNComplex(re+y, im);
}
/**
* Returns the sum of a double and a VNComplex, x+y.
* @param x A double value.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x+y.
*/
public static VNComplex plus(double x, VNComplex y)
{
return new VNComplex(x+y.re, y.im);
}
/**
* Returns the sum of this VNComplex and another VNComplex, this+y.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to this+y.
*/
public VNComplex plus(VNComplex y)
{
return new VNComplex(re+y.re, im+y.im);
}
/**
* Returns the sum of a VNComplex and a double, x+y.
* @param x A VNComplex object.
* @param y A double value.
* @return A newly constructed VNComplex initialized to x+y.
*/
public static VNComplex plus(VNComplex x, double y)
{
return new VNComplex(x.re+y, x.im);
}
/**
* Returns the sum of two VNComplex objects, x+y.
* @param x A VNComplex object.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x+y.
*/
public static VNComplex plus(VNComplex x, VNComplex y)
{
return new VNComplex(x.re+y.re, x.im+y.im);
}
/**
* Returns the sum of this VNComplex and a double, x+this.
* @param x A double value.
* @return A newly constructed VNComplex initialized to x+this.
*/
public VNComplex plusReverse(double x)
{
return new VNComplex(re+x, im);
}
/**
* Returns the VNComplex z raised to the x power,
* with a branch cut for the first parameter (z) along the
* negative real axis.
* @param z A VNComplex object.
* @param x A double value.
* @return A newly constructed VNComplex initialized to z to the power x.
*/
public static VNComplex pow(VNComplex z, double x)
{
double absz = abs(z);
VNComplex result = new VNComplex();
if (absz == 0.0) {
result = z;
} else {
double a = argument(z);
double e = Math.pow(absz, x);
result.re = e*Math.cos(x*a);
result.im = e*Math.sin(x*a);
}
return result;
}
/**
* Returns the VNComplex x raised to the VNComplex y power.
* @param x A VNComplex object.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized
* to xy.
*/
public static VNComplex pow(VNComplex x, VNComplex y)
{
return exp(times(y,log(x)));
}
/**
* Returns the real part of a VNComplex object.
* @return The real part of z.
*/
public double real()
{
return re;
}
/**
* Returns the real part of a VNComplex object.
* @param z A VNComplex object.
* @return The real part of z.
*/
public static double real(VNComplex z)
{
return z.re;
}
/**
* Returns the sine of a VNComplex.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to sine of the argument.
*/
public static VNComplex sin(VNComplex z)
{
// sin(z) = -i*sinh(i*z)
VNComplex iz = new VNComplex(-z.im,z.re);
VNComplex s = sinh(iz);
double re = s.im;
s.im = -s.re;
s.re = re;
return s;
}
/**
* Returns the hyperbolic sine of a VNComplex.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to hyperbolic
* sine of the argument.
*/
public static VNComplex sinh(VNComplex z)
{
double coshx = Sfun.cosh(z.re);
double sinhx = Sfun.sinh(z.re);
double cosy = Math.cos(z.im);
double siny = Math.sin(z.im);
boolean infiniteX = Double.isInfinite(coshx);
boolean infiniteY = Double.isInfinite(z.im);
VNComplex result;
if (z.im == 0) {
result = new VNComplex(Sfun.sinh(z.re));
} else {
// A&S 4.5.49
result = new VNComplex(sinhx*cosy, coshx*siny);
if (infiniteY) {
result.im = Double.NaN;
if (z.re == 0) result.re = 0;
}
if (infiniteX) {
result.re = z.re*cosy;
result.im = z.re*siny;
if (z.im == 0) result.im = 0;
if (infiniteY) result.re = z.im;
}
}
return result;
}
/**
* Returns the square root of a VNComplex,
* with a branch cut along the negative real axis.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized
* to square root of z. Its real part is
* non-negative.
*/
public static VNComplex sqrt(VNComplex z)
{
VNComplex result = new VNComplex();
if (Double.isInfinite(z.im)) {
result.re = Double.POSITIVE_INFINITY;
result.im = z.im;
} else if (Double.isNaN(z.re)) {
result.re = result.im = Double.NaN;
} else if (Double.isNaN(z.im)) {
if (Double.isInfinite(z.re)) {
if (z.re > 0) {
result.re = z.re;
result.im = z.im;
} else {
result.re = z.im;
result.im = Double.POSITIVE_INFINITY;
}
} else {
result.re = result.im = Double.NaN;
}
} else {
// Numerically correct version of formula 3.7.27
// in the NBS Hanbook, as suggested by Pete Stewart.
double t = abs(z);
if (Math.abs(z.re) <= Math.abs(z.im)) {
// No cancellation in these formulas
result.re = Math.sqrt(0.5*(t+z.re));
result.im = Math.sqrt(0.5*(t-z.re));
} else {
// Stable computation of the above formulas
if (z.re > 0) {
result.re = t + z.re;
result.im = Math.abs(z.im)*Math.sqrt(0.5/result.re);
result.re = Math.sqrt(0.5*result.re);
} else {
result.im = t - z.re;
result.re = Math.abs(z.im)*Math.sqrt(0.5/result.im);
result.im = Math.sqrt(0.5*result.im);
}
}
if (z.im < 0)
result.im = -result.im;
}
return result;
}
/**
* Returns the tangent of a VNComplex.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized
* to tangent of the argument.
*/
public static VNComplex tan(VNComplex z)
{
// tan = -i*tanh(i*z)
VNComplex iz = new VNComplex(-z.im,z.re);
VNComplex s = tanh(iz);
double re = s.im;
s.im = -s.re;
s.re = re;
return s;
}
/**
* Returns the hyperbolic tanh of a VNComplex.
* @param z A VNComplex object.
* @return A newly constructed VNComplex initialized to
* the hyperbolic tangent of the argument.
*/
public static VNComplex tanh(VNComplex z)
{
double sinh2x = Sfun.sinh(2*z.re);
if (z.im == 0) {
return new VNComplex(Sfun.tanh(z.re));
} else if (sinh2x == 0) {
return new VNComplex(0,Math.tan(z.im));
}
double cosh2x = Sfun.cosh(2*z.re);
double cos2y = Math.cos(2*z.im);
double sin2y = Math.sin(2*z.im);
boolean infiniteX = Double.isInfinite(cosh2x);
// Workaround for bug in JDK 1.2beta4
if (Double.isInfinite(z.im) || Double.isNaN(z.im)) {
cos2y = sin2y = Double.NaN;
}
if (infiniteX)
return new VNComplex(z.re > 0 ? 1 : -1);
// A&S 4.5.51
double den = (cosh2x + cos2y);
return new VNComplex(sinh2x/den, sin2y/den);
}
/**
* Returns the product of this VNComplex object and a double, this*y.
* @param y A double value.
* @return A newly constructed VNComplex initialized to this*y.
*/
public VNComplex times(double y)
{
return new VNComplex(re*y, im*y);
}
/**
* Returns the product of a double and a VNComplex object, x*y.
* @param x A double value.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x*y.
*/
public static VNComplex times(double x, VNComplex y)
{
return new VNComplex(x*y.re, x*y.im);
}
/**
* Returns the product of this VNComplex object and another VNComplex object, this*y.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to this*y.
*/
public VNComplex times(VNComplex y)
{
return times(this,y);
}
/**
* Returns the product of a VNComplex object and a double, x*y.
* @param x A VNComplex object.
* @param y A double value.
* @return A newly constructed VNComplex initialized to x*y.
*/
public static VNComplex times(VNComplex x, double y)
{
return new VNComplex(x.re*y, x.im*y);
}
/**
* Returns the product of two VNComplex objects, x*y.
* @param x A VNComplex object.
* @param y A VNComplex object.
* @return A newly constructed VNComplex initialized to x*y.
*/
public static VNComplex times(VNComplex x, VNComplex y)
{
VNComplex t = new VNComplex(x.re*y.re-x.im*y.im, x.re*y.im+x.im*y.re);
if (Double.isNaN(t.re) && Double.isNaN(t.im))
timesNaN(x, y, t);
return t;
}
/**
* Recovers infinities when computed x*y = NaN+i*NaN.
* This code is not part of times(), so that times
* could be inlined by an optimizing compiler.
*
* This algorithm is adapted from the C9x Annex G:
* "IEC 559-compatible VNComplex arithmetic."
* @param x First VNComplex operand.
* @param y Second VNComplex operand.
* @param t The product x*y, computed without regard to NaN.
* The real and/or the imaginary part of t is
* expected to be NaN.
* @return The corrected product of x*y.
*/
private static void timesNaN(VNComplex x, VNComplex y, VNComplex t)
{
boolean recalc = false;
double a = x.re;
double b = x.im;
double c = y.re;
double d = y.im;
if (Double.isInfinite(a) || Double.isInfinite(b)) {
// x is infinite
a = copysign(Double.isInfinite(a)?1.0:0.0, a);
b = copysign(Double.isInfinite(b)?1.0:0.0, b);
if (Double.isNaN(c)) c = copysign(0.0, c);
if (Double.isNaN(d)) d = copysign(0.0, d);
recalc = true;
}
if (Double.isInfinite(c) || Double.isInfinite(d)) {
// x is infinite
a = copysign(Double.isInfinite(c)?1.0:0.0, c);
b = copysign(Double.isInfinite(d)?1.0:0.0, d);
if (Double.isNaN(a)) a = copysign(0.0, a);
if (Double.isNaN(b)) b = copysign(0.0, b);
recalc = true;
}
if (!recalc) {
if (Double.isInfinite(a*c) || Double.isInfinite(b*d) ||
Double.isInfinite(a*d) || Double.isInfinite(b*c)) {
// Change all NaNs to 0
if (Double.isNaN(a)) a = copysign(0.0, a);
if (Double.isNaN(b)) b = copysign(0.0, b);
if (Double.isNaN(c)) c = copysign(0.0, c);
if (Double.isNaN(d)) d = copysign(0.0, d);
recalc = true;
}
}
if (recalc) {
t.re = Double.POSITIVE_INFINITY * (a*c - b*d);
t.im = Double.POSITIVE_INFINITY * (a*d + b*c);
}
}
/**
* Returns the product of a double and this VNComplex, x*this.
* @param y A double value.
* @return A newly constructed VNComplex initialized to x*this.
*/
public VNComplex timesReverse(double x)
{
return new VNComplex(x*re, x*im);
}
/**
* Returns a String representation for the specified VNComplex.
* @return A String representation for this object.
*/
public String toString()
{
if (im == 0.0)
return String.valueOf(re);
if (re == 0.0)
return String.valueOf(im) + suffix;
String sign = (im < 0.0) ? "" : "+";
return (String.valueOf(re) + sign + String.valueOf(im) + suffix);
}
/**
* Parses a string into a VNComplex.
* @param s The string to be parsed.
* @return A newly constructed VNComplex initialized to the value represented
* by the string argument.
* @exception NumberFormatException If the string does not contain a parsable VNComplex number.
* @exception NullPointerException If the input argument is null.
*/
public static VNComplex valueOf(String s) throws NumberFormatException
{
String input = s.trim();
int iBeginNumber = 0;
VNComplex z = new VNComplex();
int state = 0;
int sign = 1;
boolean haveRealPart = false;
/*
* state values
* 0 Initial State
* 1 After Initial Sign
* 2 In integer part
* 3 In fractional part
* 4 In exponential part (after 'e' but fore sign or digits)
* 5 In exponential digits
*/
for (int k = 0; k < input.length(); k++) {
char ch = input.charAt(k);
switch (ch) {
case '0': case '1': case '2': case '3': case '4':
case '5': case '6': case '7': case '8': case '9':
if (state == 0 || state == 1) {
state = 2;
} else if (state == 4) {
state = 5;
}
break;
case '-':
case '+':
sign = ((ch=='+') ? 1 : -1);
if (state == 0) {
state = 1;
} else if (state == 4) {
state = 5;
} else {
if (!haveRealPart) {
// have the real part of the number
z.re = Double.valueOf(input.substring(iBeginNumber,k)).doubleValue();
haveRealPart = true;
// perpare to part the imaginary part
iBeginNumber = k;
state = 1;
} else {
throw new NumberFormatException(input);
}
}
break;
case '.':
if (state == 0 || state == 1 || state == 2)
state = 3;
else
throw new NumberFormatException(input);
break;
case 'i': case 'I':
case 'j': case 'J':
if (k+1 != input.length()) {
throw new NumberFormatException(input);
} else if (state == 0 || state == 1) {
z.im = sign;
return z;
} else if (state == 2 || state == 3 || state == 5) {
z.im = Double.valueOf(input.substring(iBeginNumber,k)).doubleValue();
return z;
} else {
throw new NumberFormatException(input);
}
case 'e': case 'E': case 'd': case 'D':
if (state == 2 || state == 3) {
state = 4;
} else {
throw new NumberFormatException(input);
}
break;
default:
throw new NumberFormatException(input);
}
}
if (!haveRealPart) {
z.re = Double.valueOf(input).doubleValue();
return z;
} else {
throw new NumberFormatException(input);
}
}
}