Skip to content

Commit

Permalink
use apfloat hypergeometricU for evaluations
Browse files Browse the repository at this point in the history
  • Loading branch information
@axkr
axkr committed Jan 14, 2024
1 parent c94a4d0 commit dc20fb2
Show file tree
Hide file tree
Showing 6 changed files with 103 additions and 52 deletions.
88 changes: 59 additions & 29 deletions .../matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1389,26 +1389,28 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
return F.Times(F.Exp(z), F.Gamma(F.Subtract(F.C1, a), z));
}
try {
IExpr n = engine.evaluate(F.Subtract(b, a));
if (n.isInteger()) {
if (n.isOne()) {
// b==a+1 ==> z^(-a)
return F.Power(z, a.negate());
}
int nInt = n.toIntDefault();
if (nInt > 0) {
int nMinus1 = nInt - 1;
// Sum((Binomial(-1+n, -1-k+n)*Pochhammer(a, k))/z^k, {k, 0, n-1}) / z^a
return F.Times(F.Power(z, a.negate()), //
F.intSum( //
k -> F.Times(F.Binomial(nMinus1, nMinus1 - k), F.Pochhammer(a, F.ZZ(k)),
F.Power(z, -k)), //
0, nMinus1));
{
IExpr n = engine.evaluate(F.Subtract(b, a));
if (n.isInteger()) {
if (n.isOne()) {
// b==a+1 ==> z^(-a)
return F.Power(z, a.negate());
}
int nInt = n.toIntDefault();
if (nInt > 0) {
int nMinus1 = nInt - 1;
// Sum((Binomial(-1+n, -1-k+n)*Pochhammer(a, k))/z^k, {k, 0, n-1}) / z^a
return F.Times(F.Power(z, a.negate()), //
F.intSum( //
k -> F.Times(F.Binomial(nMinus1, nMinus1 - k), F.Pochhammer(a, F.ZZ(k)),
F.Power(z, -k)), //
0, nMinus1));
}
}
}
if (engine.isArbitraryMode()) {
try {
IExpr res = a.hypergeometric1F1(b, z);
IExpr res = a.hypergeometricU(b, z);
if (res.isNumber()) {
return res;
}
Expand All @@ -1418,22 +1420,50 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
} else if (engine.isDoubleMode()) {
double aDouble = Double.NaN;
double bDouble = Double.NaN;
double zDouble = Double.NaN;
try {
aDouble = a.evalf();
bDouble = b.evalf();
zDouble = z.evalf();
return F.complexNum(HypergeometricJS.hypergeometricU(new Complex(aDouble),
new Complex(bDouble), new Complex(zDouble)));
IExpr res = a.hypergeometricU(b, z);
if (res.isNumber()) {
return res;
}
} catch (ValidateException ve) {
Errors.printMessage(ast.topHead(), ve, engine);
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
Complex ac = a.evalfc();
Complex bc = b.evalfc();
Complex zc = z.evalfc();
return F.complexNum(HypergeometricJS.hypergeometricU(ac, bc, zc));

// double aDouble = Double.NaN;
// double bDouble = Double.NaN;
// double zDouble = Double.NaN;
// try {
// aDouble = a.evalf();
// bDouble = b.evalf();
// zDouble = z.evalf();
// return F.complexNum(HypergeometricJS.hypergeometricU(new Complex(aDouble),
// new Complex(bDouble), new Complex(zDouble)));
// } catch (ValidateException ve) {
// Errors.printMessage(ast.topHead(), ve, engine);
// }
// Complex ac = a.evalfc();
// Complex bc = b.evalfc();
// Complex zc = z.evalfc();
// return F.complexNum(HypergeometricJS.hypergeometricU(ac, bc, zc));
}

if (a.isInteger() && a.isPositive() && (!b.isNumber() || !z.isNumber())) {
IInteger n = (IInteger) a;
ISymbol k = F.Dummy("k");
// (Gamma(-1+b,z)*z^(1-b)*E^z*LaguerreL(-1+n,1-b,-z)-Sum(LaguerreL(-1-k+n,-b+k+1,-z)/k*LaguerreL(-
// 1+k,-1+b-k,z),{k,1,-1+n}))/Pochhammer(2-b,-1+n)
return F.Times(F.Power(F.Pochhammer(F.Subtract(F.C2, b), F.Plus(F.CN1, n)), F.CN1),
F.Subtract(
F.Times(F.Gamma(F.Plus(F.CN1, b), z), F.Power(z, F.Subtract(F.C1, b)), F.Exp(z),
F.LaguerreL(F.Plus(F.CN1, n), F.Subtract(F.C1, b), F.Negate(z))),
F.Sum(
F.Times(F.Power(k, F.CN1),
F.LaguerreL(F.Plus(F.CN1, F.Negate(k), n), F.Plus(F.Negate(b), k, F.C1),
F.Negate(z)),
F.LaguerreL(F.Plus(F.CN1, k), F.Plus(F.CN1, b, F.Negate(k)), z)),
F.list(k, F.C1, F.Plus(F.CN1, n)))));
}
} catch (ThrowException te) {
LOGGER.debug("HypergeometricU.evaluate() failed", te);
Expand Down
20 changes: 14 additions & 6 deletions ...id_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -722,14 +722,22 @@ public IExpr hypergeometric2F1(IExpr arg2, IExpr arg3, IExpr arg4) {

@Override
public IExpr hypergeometricU(IExpr arg2, IExpr arg3) {
try {
return F.complexNum(HypergeometricJS.hypergeometricU(fComplex, //
arg2.evalfc(), //
arg3.evalfc()));
} catch (RuntimeException e) {
// try as computation with complex numbers
if (arg2 instanceof INumber && arg3 instanceof INumber) {
Apcomplex hypergeometricU = EvalEngine.getApfloatDouble().hypergeometricU(apcomplexValue(),
((INumber) arg2).apcomplexValue(), ((INumber) arg3).apcomplexValue());
return F.complexNum(hypergeometricU.real().doubleValue(),
hypergeometricU.imag().doubleValue());
}
return IComplexNum.super.hypergeometricU(arg2, arg3);

// try {
// return F.complexNum(HypergeometricJS.hypergeometricU(fComplex, //
// arg2.evalfc(), //
// arg3.evalfc()));
// } catch (RuntimeException e) {
// // try as computation with complex numbers
// }
// return IComplexNum.super.hypergeometricU(arg2, arg3);
}

/** {@inheritDoc} */
Expand Down
24 changes: 19 additions & 5 deletions ...a_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -564,13 +564,27 @@ public IExpr hypergeometric2F1(IExpr arg2, IExpr arg3, IExpr arg4) {
@Override
public IExpr hypergeometricU(IExpr arg2, IExpr arg3) {
try {
return F.complexNum(HypergeometricJS.hypergeometricU(new Complex(value), //
arg2.evalfc(), //
arg3.evalfc()));
} catch (RuntimeException e) {
// try as computation with complex numbers
Apcomplex hypergeometricU = EvalEngine.getApfloatDouble().hypergeometricU(apcomplexValue(),
((INumber) arg2).apcomplexValue(), ((INumber) arg3).apcomplexValue());
return F.complexNum(hypergeometricU.real().doubleValue(),
hypergeometricU.imag().doubleValue());
} catch (ArithmeticException aex) {
if (aex.getMessage().equals("Division by zero")) {
return F.ComplexInfinity;
} else {
// aex.printStackTrace();
}
}
return INum.super.hypergeometricU(arg2, arg3);

// try {
// return F.complexNum(HypergeometricJS.hypergeometricU(new Complex(value), //
// arg2.evalfc(), //
// arg3.evalfc()));
// } catch (RuntimeException e) {
// // try as computation with complex numbers
// }
// return INum.super.hypergeometricU(arg2, arg3);
}

/** {@inheritDoc} */
Expand Down
5 changes: 1 addition & 4 deletions ...core/src/main/java/org/matheclipse/core/reflection/system/rules/HypergeometricURules.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,6 @@ public class HypergeometricURules {
Times(Exp(z),Power(z,Subtract(C1,m)),Gamma(Plus(CN1,m),z))),
// HypergeometricU(1/2,1,z_):=(E^(z/2)*BesselK(0,z/2))/Sqrt(Pi)
ISetDelayed(HypergeometricU(C1D2,C1,z_),
Times(Exp(Times(C1D2,z)),Power(Pi,CN1D2),BesselK(C0,Times(C1D2,z)))),
// HypergeometricU(n_Integer,b_,z_):=-1/Pochhammer(2-b,-1+n)*(-Gamma(-1+b,z)*z^(1-b)*E^z*LaguerreL(-1+n,1-b,-z)+Sum((LaguerreL(-1-k+n,-b+k+1,-z)*LaguerreL(-1+k,-1+b-k,z))/k,{k,1,-1+n}))/;n>0
ISetDelayed(HypergeometricU($p(n, Integer),b_,z_),
Condition(Times(CN1,Power(Pochhammer(Subtract(C2,b),Plus(CN1,n)),CN1),Plus(Times(CN1,Gamma(Plus(CN1,b),z),Power(z,Subtract(C1,b)),Exp(z),LaguerreL(Plus(CN1,n),Subtract(C1,b),Negate(z))),Sum(Times(Power(k,CN1),LaguerreL(Plus(CN1,Negate(k),n),Plus(Negate(b),k,C1),Negate(z)),LaguerreL(Plus(CN1,k),Plus(CN1,b,Negate(k)),z)),list(k,C1,Plus(CN1,n))))),Greater(n,C0)))
Times(Exp(Times(C1D2,z)),Power(Pi,CN1D2),BesselK(C0,Times(C1D2,z))))
);
}
14 changes: 9 additions & 5 deletions ...library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11655,8 +11655,14 @@ public void testHypergeometricPFQ() {
@Test
public void testHypergeometricU() {
// TODO throws hypergeometric function pole
// check("HypergeometricU(3, 2, 1.0)", //
// "0.105479");
// https://github.com/mtommila/apfloat/issues/36
check("HypergeometricU(3.0, 1.0, 0.0)", //
"HypergeometricU(3.0,1.0,0.0)");
check("HypergeometricU(3, 2, 1.0)", //
"0.105479");
check("N(HypergeometricU(3, 2, 1),50)", //
"0.10547895651520888848838225094608093588873320977117");


check("D(HypergeometricU(a,b,x), x)", //
"-a*HypergeometricU(1+a,1+b,x)");
Expand All @@ -11673,9 +11679,7 @@ public void testHypergeometricU() {
+ "0.76904+I*0.163012,1.05965+I*0.56395,1.44956+I*1.52035,1.97105+I*4.83136,ComplexInfinity," //
+ "2.45436,0.688641,0.312167,0.173724,0.1086,0.0732253,0.0520871,0.0385635}");
check("Table( HypergeometricU(3, 1.0, x), {x,-2.0,2,0.25})", //
"{0.0852414+I*0.212584,0.0312283+I*0.264433,-0.0527303+I*0.306681,-0.171748+I*0.323467,-0.325706+I*0.288932," //
+ "-0.500123+I*0.162311,-0.642219+I*(-0.119092),-0.575265+I*(-0.649898),Indeterminate,0.214115,0.105593,0.0644474," //
+ "0.0436079,0.0314298,0.0236577,0.0183874,0.0146502}");
"{0.0852414+I*0.212584,0.0312283+I*0.264433,-0.0527303+I*0.306681,-0.171748+I*0.323467,-0.325706+I*0.288932,-0.500123+I*0.162311,-0.642219+I*(-0.119092),-0.575265+I*(-0.649898),HypergeometricU(3.0,1.0,0.0),0.214115,0.105593,0.0644474,0.0436079,0.0314298,0.0236577,0.0183874,0.0146502}");
}

@Test
Expand Down
4 changes: 1 addition & 3 deletions symja_android_library/rules/HypergeometricURules.m
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,4 @@
{
HypergeometricU(1,m_,z_) := E^z*z^(1-m)*Gamma(-1+m, z),
HypergeometricU(1/2,1,z_) := (E^(z/2)*BesselK(0, z/2))/Sqrt(Pi),
HypergeometricU(n_Integer, b_, z_) := (-1/Pochhammer(2-b, n-1))*(-Gamma(b - 1, z)*z^(1 - b)*E^z*LaguerreL(n-1, 1-b, -z)+Sum((1/k)*LaguerreL(n-k-1, -b+k+1, -z)*LaguerreL(k-1, b-k-1, z), {k, 1, n - 1}))
/; n>0
HypergeometricU(1/2,1,z_) := (E^(z/2)*BesselK(0, z/2))/Sqrt(Pi)
}

0 comments on commit dc20fb2

Please sign in to comment.