Skip to content

Commit

Permalink
Fix bugs in ArcTan(x,y) in calling atan2()
Browse files Browse the repository at this point in the history 
- Hipparchus-Math/hipparchus#140
- Note: the first argument of atan2 is y and the second is x.
  • Loading branch information
@axkr
axkr committed Jun 18, 2021
1 parent 441ed19 commit 88f844b
Show file tree
Hide file tree
Showing 4 changed files with 138 additions and 15 deletions.
78 changes: 67 additions & 11 deletions ...ibrary/matheclipse-core/src/main/java/org/matheclipse/core/builtin/ExpTrigsFunctions.java
View file
Expand Up @@ -40,7 +40,9 @@
import org.apfloat.ApfloatMath;
import org.apfloat.Apint;
import org.hipparchus.complex.Complex;
import org.hipparchus.special.elliptic.carlson.CarlsonEllipticIntegral;
import org.hipparchus.util.FastMath;
import org.matheclipse.core.builtin.functions.EllipticIntegralsJS;
import org.matheclipse.core.eval.EvalEngine;
import org.matheclipse.core.eval.exception.ValidateException;
import org.matheclipse.core.eval.interfaces.AbstractArg1;
Expand Down Expand Up @@ -994,7 +996,7 @@ public IExpr e1ObjArg(final IExpr arg1) {
}

@Override
public IExpr e2ObjArg(final IExpr x, final IExpr y) {
public IExpr e2ObjArg(IExpr x, IExpr y) {
if (x.isZero() && y.isRealResult()) {
if (y.isZero()) {
return S.Indeterminate;
Expand Down Expand Up @@ -1029,6 +1031,43 @@ public IExpr e2ObjArg(final IExpr x, final IExpr y) {
if (x.isNegativeInfinity()) {
return F.Times(F.Subtract(F.Times(F.C2, F.UnitStep(F.Re(y))), F.C1), S.Pi);
}
EvalEngine engine = EvalEngine.get();
if (engine.isNumericMode()) {
if (engine.isArbitraryMode()) {
x = engine.evalN(x);
y = engine.evalN(y);

if (x.isReal() && y.isReal()) {
long precision = engine.getNumericPrecision();
Apfloat xa = ((ISignedNumber) x).apfloatValue(precision);
Apfloat ya = ((ISignedNumber) y).apfloatValue(precision);
return F.num(ApfloatMath.atan2(ya, xa));
}
if (x.isNumber() && y.isNumber()) {
long precision = engine.getNumericPrecision();
x = ((INumber) x).apcomplexNumValue(precision);
y = ((INumber) y).apcomplexNumValue(precision);
return y.atan2(x);
}
return F.NIL;
}
double xd = Double.NaN;
double yd = Double.NaN;
try {
xd = x.evalDouble();
yd = y.evalDouble();
} catch (ValidateException ve) {
}
if (Double.isNaN(xd) || Double.isNaN(yd)) {
Complex xc = x.evalComplex();
Complex yc = y.evalComplex();
// the first argument of atan2 is y and the second is x.
return F.complexNum(yc.atan2(xc));
} else {
// the first argument of atan2 is y and the second is x.
return F.num(Math.atan2(yd, xd));
}
}
return F.NIL;
}

Expand All @@ -1053,8 +1092,9 @@ public IExpr e1DblComArg(final IComplexNum val) {
}

@Override
public IExpr e2DblArg(final INum d0, final INum d1) {
return F.num(Math.atan2(d0.getRealPart(), d1.getRealPart()));
public IExpr e2DblArg(final INum x, final INum y) {
// the first argument of atan2 is y and the second is x.
return F.num(Math.atan2(y.getRealPart(), x.getRealPart()));
}

@Override
Expand All @@ -1073,7 +1113,8 @@ public double evalReal(final double[] stack, final int top, final int size) {
throw new UnsupportedOperationException();
}
if (size == 2) {
return Math.atan2(stack[top - 1], stack[top]);
// the first argument of atan2 is y and the second is x.
return Math.atan2(stack[top], stack[top - 1]);
}

return Math.atan(stack[top]);
Expand Down Expand Up @@ -2090,17 +2131,29 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
if (engine.isNumericMode()) {
try {
if (z.isNumber()) {
// see https://github.com/paulmasson/math/issues/24

// Re(z)>0 || (Re(z)==0&&Im(z)>=0)
if (((INumber) z).complexSign() >= 0) {
// (1/2)*(Pi - 4*ArcCot(E^z))
return F.Times.of(
engine, F.C1D2, F.Subtract(S.Pi, F.Times(F.C4, F.ArcCot(F.Power(S.E, z)))));
}
// (1/2)*(-Pi + 4*ArcTan(E^z))

// 2*ArcTan( Tan( z/2 )
return F.Times.of(
engine,
F.C1D2,
F.Plus(F.Times(F.CN1, S.Pi), F.Times(F.C4, F.ArcTan(F.Power(S.E, z)))));
F.C2,
F.ArcTan( //
F.Tanh(F.Times(F.C1D2, z)) //
));

// // (1/2)*(-Pi + 4*ArcTan(E^z))
// return F.Times.of(
// engine,
// F.C1D2,
// F.Plus(F.Times(F.CN1, S.Pi), F.Times(F.C4, F.ArcTan(F.Power(S.E,
// z)))));
}
} catch (ValidateException ve) {
if (FEConfig.SHOW_STACKTRACE) {
Expand Down Expand Up @@ -2167,11 +2220,14 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
if (engine.isNumericMode()) {
try {
if (z.isNumber()) {
// Log(Tan(Pi/4 + z/2))
return F.Log.of(
// see https://github.com/paulmasson/math/issues/24

// 2*ArcTanh( Tan( z/2 )
return F.Times.of(
engine,
F.Tan( //
F.Plus(F.Times(F.C1D4, S.Pi), F.Times(F.C1D2, z)) //
F.C2,
F.ArcTanh( //
F.Tan(F.Times(F.C1D2, z)) //
));
}
} catch (ValidateException ve) {
Expand Down
26 changes: 25 additions & 1 deletion ..._library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ApcomplexNum.java
View file
Expand Up @@ -24,7 +24,7 @@
import org.matheclipse.core.visit.IVisitorBoolean;
import org.matheclipse.core.visit.IVisitorInt;
import org.matheclipse.core.visit.IVisitorLong;
import org.matheclipse.parser.client.FEConfig;
import org.matheclipse.parser.client.FEConfig;

/**
* <code>IComplexNum</code> implementation which wraps a <code>
Expand Down Expand Up @@ -696,6 +696,30 @@ public IExpr atan() {
return valueOf(ApcomplexMath.atan(fApcomplex));
}

@Override
public IExpr atan2(IExpr value) {
if (value instanceof ApcomplexNum) {
Apcomplex th = fApcomplex;
final long precision = th.precision();
Apcomplex x = ((ApcomplexNum)value).fApcomplex;

// compute r = sqrt(x^2+y^2)
final Apcomplex r = ApcomplexMath.sqrt(x.multiply(x).add(th.multiply(th)));

if (x.real().compareTo(Apfloat.ZERO) >= 0) {
// compute atan2(y, x) = 2 atan(y / (r + x))
return valueOf(ApcomplexMath.atan(th.divide(r.add(x))).multiply(new Apfloat(2, precision)));
} else {
// compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
return valueOf(
ApcomplexMath.atan(th.divide(r.subtract(x)))
.multiply(new Apfloat(-2, precision))
.add(ApfloatMath.pi(precision)));
}
}
return IComplexNum.super.atan2(value);
}

@Override
public IExpr atanh() {
return valueOf(ApcomplexMath.atanh(fApcomplex));
Expand Down
5 changes: 3 additions & 2 deletions ...atheclipse-frontend/src/main/java/org/matheclipse/parser/client/eval/DoubleEvaluator.java
View file
Expand Up @@ -75,8 +75,9 @@ public double applyAsDouble(double arg1) {
}

@Override
public double applyAsDouble(double arg1, double arg2) {
return Math.atan2(arg1, arg2);
public double applyAsDouble(double x, double y) {
// the first argument of atan2 is y and the second is x.
return Math.atan2(y, x);
}
}

Expand Down
44 changes: 43 additions & 1 deletion ...oid_library/matheclipse-io/src/test/java/org/matheclipse/io/system/LowercaseTestCase.java
View file
Expand Up @@ -1176,6 +1176,33 @@ public void testArcTan() {
// github #110 avoid infinite recursion
// check("ArcTan(Re(Sin(3+I*2)),Im(Sin(3+I*2)))", //
// "ArcTan(Im(Sin(3+I*2))/Re(Sin(3+I*2)))");

check(
"N(ArcTan(2, 1), 50)", //
"0.4636476090008061162142562314612144020285370542861");
check(
"N(ArcTan(2, 4), 50)", //
"1.1071487177940905030170654601785370400700476454014");
check(
"N(ArcTan(2, I*4), 50)", //
"1.5707963267948966192313216916397514420985846996875+I*0.5493061443340548456976226184612628523237452789113");
check(
"N(ArcTan(I*4, 2), 50)", //
"I*(-0.54930614433405484569762261846126285232374527891137)");

check(
"N(ArcTan(2, 1))", //
"0.463648");
check(
"N(ArcTan(2, 4))", //
"1.10715");
check(
"N(ArcTan(I*4, 2))", //
"I*(-0.549306)");
check(
"N(ArcTan(2,I*4 ))", //
"1.5708+I*0.549306");

check(
"ArcTan(y_,x^2)", //
"ArcTan(y_,x^2)");
Expand Down Expand Up @@ -14839,6 +14866,9 @@ public void testGudermannian() {
check(
"N(Gudermannian(4/3), 50)", //
"1.055327395759396716725120121987678867953438481189");
check(
"N(Gudermannian(I*Pi+5), 50)", //
"1.584272016863742216280221584803507365335847712497");
check(
"Gudermannian(6/4*Pi*I)", //
"-I*Infinity");
Expand Down Expand Up @@ -14943,6 +14973,12 @@ public void testHarmonicNumber() {
}

public void testHaversine() {
check(
"N(Haversine(1/7), 50)", //
"0.00509336977669245865213693149328312706298364955021");
check(
"N(Haversine(798+I), 50)", //
"-0.27105513255154801719412701505423122282395758411+I*0.020835467253162356620662405472920178540829827673");
checkNumeric(
"Haversine(1.5)", //
"0.46463139916614854");
Expand Down Expand Up @@ -18057,7 +18093,7 @@ public void testInverseGammaRegularized() {
public void testInverseGudermannian() {
check(
"N(InverseGudermannian(4/3), 50)", //
"2.12617609078226939193623627117215618545032512233803");
"2.126176090782269391936236271172156185450325122338");
check(
"N(InverseGudermannian(4/3-2/3*I), 50)", //
"1.0709631602353002334250501606064505306184199829854+I*(-1.253841344203559414899366482437511230389585771795)");
Expand All @@ -18081,6 +18117,12 @@ public void testInverseGudermannian() {

public void testInverseHaversine() {

check(
"N(InverseHaversine(1/7), 50)", //
"0.775193373310361307204093711182473425799817434023");
check(
"N(InverseHaversine(798+I), 50)", //
"3.140338735505235220711482104291646256237738908735+I*8.067776883665043131189388013520931297785865736445");
checkNumeric(
"InverseHaversine(0.5)", //
"1.5707963267948968");
Expand Down

0 comments on commit 88f844b

Please sign in to comment.