From dc20fb20f8321d5596f864add03410a536dd0e67 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 14 Jan 2024 13:54:10 +0100 Subject: [PATCH] use apfloat hypergeometricU for evaluations --- .../core/builtin/HypergeometricFunctions.java | 88 +++++++++++++------ .../core/expression/ComplexNum.java | 20 +++-- .../org/matheclipse/core/expression/Num.java | 24 +++-- .../system/rules/HypergeometricURules.java | 5 +- .../core/system/LowercaseTestCase.java | 14 +-- .../rules/HypergeometricURules.m | 4 +- 6 files changed, 103 insertions(+), 52 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java index 87092b4568..7be6be6653 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java @@ -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; } @@ -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); diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java index 81d60683cb..25e90a78a6 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java @@ -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} */ diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java index 892d5ff206..0d3005f963 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java @@ -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} */ diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/rules/HypergeometricURules.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/rules/HypergeometricURules.java index 82ca58a3bf..0cb05342d1 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/rules/HypergeometricURules.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/rules/HypergeometricURules.java @@ -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)))) ); } diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java index ab942e2058..7c9893013f 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java @@ -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)"); @@ -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 diff --git a/symja_android_library/rules/HypergeometricURules.m b/symja_android_library/rules/HypergeometricURules.m index ba18a33737..b9079d54ff 100644 --- a/symja_android_library/rules/HypergeometricURules.m +++ b/symja_android_library/rules/HypergeometricURules.m @@ -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) } \ No newline at end of file