From 036bfb5fa1c82ed4b8dec59ff08ee18283a78475 Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Sat, 6 Dec 2025 07:14:48 -0500 Subject: [PATCH 1/2] Add RR versions of some flint special functions This way, we can go RR -> RRb -> RR instead of RR -> RRi -> RRb -> RRi -> RR --- M2/Macaulay2/d/actors3.d | 20 +++++--------- M2/Macaulay2/d/ballarith.d | 56 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 13 deletions(-) diff --git a/M2/Macaulay2/d/actors3.d b/M2/Macaulay2/d/actors3.d index 9662babf6e5..c9ac7d42287 100644 --- a/M2/Macaulay2/d/actors3.d +++ b/M2/Macaulay2/d/actors3.d @@ -891,8 +891,7 @@ regularizedGamma(e:Expr):Expr := ( when a.0 is s:RRcell do ( when a.1 - is x:RRcell do toExpr( - midpointRR(regularizedGamma(toRRi(s.v), toRRi(x.v)))) -- # typical value: regularizedGamma, RR, RR, RR + is x:RRcell do toExpr(regularizedGamma(s.v, x.v)) -- # typical value: regularizedGamma, RR, RR, RR is x:RRicell do toExpr(regularizedGamma(toRRi(s.v), x.v)) -- # typical value: regularizedGamma, RR, RRi, RRi is x:CCcell do toExpr(regularizedGamma(toCC(s.v), x.v)) -- # typical value: regularizedGamma, RR, CC, CC else WrongArgRRorCC(2)) @@ -956,7 +955,7 @@ erfc(e:Expr):Expr := ( setupfun("erfc",erfc).Protected=false; inverseErf(e:Expr):Expr := ( when e - is x:RRcell do toExpr(midpointRR(inverseErf(toRRi(x.v)))) -- # typical value: inverseErf, RR, RR + is x:RRcell do toExpr(inverseErf(x.v)) -- # typical value: inverseErf, RR, RR is x:RRicell do toExpr(inverseErf(x.v)) -- # typical value: inverseErf, RRi, RRi else WrongArgRRorRRi()); setupfun("inverseErf",inverseErf).Protected=false; @@ -976,15 +975,13 @@ BesselJ(e:Expr):Expr := ( else WrongArgRRorCC(2)) else ( when s.1 - is x:RRcell do toExpr( - midpointRR(BesselJ(toRRi(n.v), toRRi(x.v)))) + is x:RRcell do toExpr(BesselJ(toRR(n.v), x.v)) is x:RRicell do toExpr(BesselJ(toRRi(n.v,precision(x.v)),x.v)) is x:CCcell do toExpr(BesselJ(toCC(n.v), x.v )) else WrongArgRRorCC(2))) is n:RRcell do ( when s.1 - is x:RRcell do toExpr( - midpointRR(BesselJ(toRRi(n.v), toRRi(x.v)))) + is x:RRcell do toExpr(BesselJ(n.v, x.v)) is x:RRicell do toExpr(BesselJ(toRRi(n.v), x.v)) is x:CCcell do toExpr(BesselJ(toCC(n.v), x.v )) else WrongArgRRorCC(2)) @@ -1017,15 +1014,13 @@ BesselY(e:Expr):Expr := ( else WrongArgRRorCC(2)) else ( when s.1 - is x:RRcell do toExpr( - midpointRR(BesselY(toRRi(n.v), toRRi(x.v)))) + is x:RRcell do toExpr(BesselY(toRR(n.v), x.v)) is x:RRicell do toExpr(BesselY(toRRi(n.v,precision(x.v)),x.v)) is x:CCcell do toExpr(BesselY(toCC(n.v), x.v )) else WrongArgRRorCC(2))) is n:RRcell do ( when s.1 - is x:RRcell do toExpr( - midpointRR(BesselY(toRRi(n.v), toRRi(x.v)))) + is x:RRcell do toExpr(BesselY(n.v, x.v)) is x:RRicell do toExpr(BesselY(toRRi(n.v), x.v)) is x:CCcell do toExpr(BesselY(toCC(n.v), x.v )) else WrongArgRRorCC(2)) @@ -1102,8 +1097,7 @@ regularizedBeta(xx:Expr,yy:Expr,zz:Expr):Expr := ( when yy is y:RRcell do ( when zz - is z:RRcell do toExpr( - midpointRR(regularizedBeta(toRRi(x.v), toRRi(y.v), toRRi(z.v)))) -- # typical value: regularizedBeta, RR, RR, RR, RR + is z:RRcell do toExpr(regularizedBeta(x.v, y.v, z.v)) -- # typical value: regularizedBeta, RR, RR, RR, RR is z:RRicell do toExpr(regularizedBeta(toRRi(x.v), toRRi(y.v), z.v)) -- # typical value: regularizedBeta, RR, RR, RRi, RRi is z:CCcell do toExpr(regularizedBeta(toCC(x.v), toCC(y.v), z.v)) -- # typical value: regularizedBeta, RR, RR, CC, CC else WrongArgRRorCC(3)) diff --git a/M2/Macaulay2/d/ballarith.d b/M2/Macaulay2/d/ballarith.d index dc912047828..1ea927a7862 100644 --- a/M2/Macaulay2/d/ballarith.d +++ b/M2/Macaulay2/d/ballarith.d @@ -52,6 +52,11 @@ toRR(x:RRball, prec:ulong):RR := ( Ccode(int, "arf_get_mpfr(", y, ", arb_midref(", x, "), MPFR_RNDN)"); moveToRRandclear(y)); +moveToRRandclear(x:RRball, prec:ulong):RR := ( + r := toRR(x, prec); + clear(x); + r); + toRRi(x:RRball, prec:ulong):RRi := ( y := newRRimutable(prec); Ccode(void, "arb_get_interval_mpfr((mpfr_ptr)&", y, @@ -89,6 +94,17 @@ export Gamma(z:RRi,w:RRi):RRi := ( clear(y); moveToRRiandclear(r, prec)); +export regularizedGamma(z:RR,w:RR):RR := ( + prec := min(precision(z), precision(w)); + x := toRRball(z); + y := toRRball(w); + r := newRRball(); + Ccode(void, "arb_hypgeom_gamma_upper(", r, ", ", x, ", ", y, ", 1, ", + prec, ")"); + clear(x); + clear(y); + moveToRRandclear(r, prec)); + export regularizedGamma(z:RRi,w:RRi):RRi := ( prec := min(precision(z), precision(w)); x := toRRball(z); @@ -135,6 +151,13 @@ export erfc(x:RRi):RRi := ( clear(y); moveToRRiandclear(r, precision(x))); +export inverseErf(x:RR):RR := ( + y := toRRball(x); + r := newRRball(); + Ccode(void, "arb_hypgeom_erfinv(", r, ", ", y, ", ", precision(x), ")"); + clear(y); + moveToRRandclear(r, precision(x))); + export inverseErf(x:RRi):RRi := ( y := toRRball(x); r := newRRball(); @@ -142,6 +165,16 @@ export inverseErf(x:RRi):RRi := ( clear(y); moveToRRiandclear(r, precision(x))); +export BesselJ(z:RR,w:RR):RR := ( + prec := min(precision(z), precision(w)); + x := toRRball(z); + y := toRRball(w); + r := newRRball(); + Ccode(void, "arb_hypgeom_bessel_j(", r, ", ", x, ", ", y, ", ", prec, ")"); + clear(x); + clear(y); + moveToRRandclear(r, prec)); + export BesselJ(z:RRi,w:RRi):RRi := ( prec := min(precision(z), precision(w)); x := toRRball(z); @@ -152,6 +185,16 @@ export BesselJ(z:RRi,w:RRi):RRi := ( clear(y); moveToRRiandclear(r, prec)); +export BesselY(z:RR,w:RR):RR := ( + prec := min(precision(z), precision(w)); + x := toRRball(z); + y := toRRball(w); + r := newRRball(); + Ccode(void, "arb_hypgeom_bessel_y(", r, ", ", x, ", ", y, ", ", prec, ")"); + clear(x); + clear(y); + moveToRRandclear(r, prec)); + export BesselY(z:RRi,w:RRi):RRi := ( prec := min(precision(z), precision(w)); x := toRRball(z); @@ -174,6 +217,19 @@ export Beta(z:RRi,w:RRi):RRi := ( clear(y); moveToRRiandclear(r, prec)); +export regularizedBeta(u:RR,v:RR,w:RR):RR := ( + prec := min(min(precision(u), precision(v)), precision(w)); + x := toRRball(u); + y := toRRball(v); + z := toRRball(w); + r := newRRball(); + Ccode(void, "arb_hypgeom_beta_lower(", r, ", ", y, ", ", z, ", ", x, + ", 1, ", prec, ")"); + clear(x); + clear(y); + clear(z); + moveToRRandclear(r, prec)); + export regularizedBeta(u:RRi,v:RRi,w:RRi):RRi := ( prec := min(min(precision(u), precision(v)), precision(w)); x := toRRball(u); From 1732cc5dd2f5041f6c9090d76920b402db3e6a7e Mon Sep 17 00:00:00 2001 From: Doug Torrance Date: Mon, 1 Dec 2025 17:12:19 -0500 Subject: [PATCH 2/2] Add polylog function --- M2/Macaulay2/d/actors3.d | 28 ++++++++++++++++++++++++++ M2/Macaulay2/d/ballarith.d | 30 ++++++++++++++++++++++++++++ M2/Macaulay2/m2/exports.m2 | 1 + M2/Macaulay2/m2/typicalvalues.m2 | 2 +- M2/Macaulay2/tests/normal/numbers.m2 | 4 ++++ 5 files changed, 64 insertions(+), 1 deletion(-) diff --git a/M2/Macaulay2/d/actors3.d b/M2/Macaulay2/d/actors3.d index c9ac7d42287..0000252e2a8 100644 --- a/M2/Macaulay2/d/actors3.d +++ b/M2/Macaulay2/d/actors3.d @@ -1145,6 +1145,34 @@ regularizedBeta(e:Expr):Expr := ( else WrongNumArgs(3)); setupfun("regularizedBeta",regularizedBeta).Protected=false; +polylog(e:Expr):Expr := ( + when e + is a:Sequence do ( + if length(a) == 2 then ( + when a.0 + is x:RRcell do ( + when a.1 + is y:RRcell do ( -- # typical value: polylog, RR, RR, RR + if y.v < 1 then toExpr(polylog(x.v, y.v)) + else toExpr(polylog(toCC(x.v), toCC(y.v)))) + is y:RRicell do toExpr(polylog(toRRi(x.v), y.v)) -- # typical value: polylog, RR, RRi, RRi + is y:CCcell do toExpr(polylog(toCC(x.v), y.v)) -- # typical value: polylog, RR, CC, CC + else WrongArgRRorCC(2)) + is x:RRicell do ( + when a.1 + is y:RRcell do toExpr(polylog(x.v, toRRi(y.v))) -- # typical value: polylog, RRi, RR, RRi + is y:RRicell do toExpr(polylog(x.v, y.v)) -- # typical value: polylog, RRi, RRi, RRi + else WrongArgRRorCC(2)) + is x:CCcell do ( + when a.1 + is y:RRcell do toExpr(polylog(x.v, toCC(y.v))) -- # typical value: polylog, CC, RR, CC + is y:CCcell do toExpr(polylog(x.v, y.v)) -- # typical value: polylog, CC, CC, CC + else WrongArgRRorCC(2)) + else WrongArgRRorCC(1)) + else WrongNumArgs(2)) + else WrongNumArgs(2)); +setupfun("polylog", polylog).Protected=false; + cosh(e:Expr):Expr := ( when e is x:CCcell do toExpr(cosh(x.v)) -- # typical value: cosh, CC, CC diff --git a/M2/Macaulay2/d/ballarith.d b/M2/Macaulay2/d/ballarith.d index 1ea927a7862..14bb5dedc5b 100644 --- a/M2/Macaulay2/d/ballarith.d +++ b/M2/Macaulay2/d/ballarith.d @@ -243,6 +243,26 @@ export regularizedBeta(u:RRi,v:RRi,w:RRi):RRi := ( clear(z); moveToRRiandclear(r, prec)); +export polylog(z:RR, w:RR):RR := ( + prec := min(precision(z), precision(w)); + x := toRRball(z); + y := toRRball(w); + r := newRRball(); + Ccode(void, "arb_polylog(", r, ", ", x, ", ", y, ", ", prec, ")"); + clear(x); + clear(y); + moveToRRandclear(r, prec)); + +export polylog(z:RRi, w:RRi):RRi := ( + prec := min(precision(z), precision(w)); + x := toRRball(z); + y := toRRball(w); + r := newRRball(); + Ccode(void, "arb_polylog(", r, ", ", x, ", ", y, ", ", prec, ")"); + clear(x); + clear(y); + moveToRRiandclear(r, prec)); + ------------ -- CCBall -- ------------ @@ -389,3 +409,13 @@ export regularizedBeta(u:CC,v:CC,w:CC):CC := ( clear(y); clear(z); moveToCCandclear(r, prec)); + +export polylog(z:CC, w:CC):CC := ( + prec := min(precision(z), precision(w)); + x := toCCball(z); + y := toCCball(w); + r := newCCball(); + Ccode(void, "acb_polylog(", r, ", ", x, ", ", y, ", ", prec, ")"); + clear(x); + clear(y); + moveToCCandclear(r, prec)); diff --git a/M2/Macaulay2/m2/exports.m2 b/M2/Macaulay2/m2/exports.m2 index 56e7ec34d30..785291ec985 100644 --- a/M2/Macaulay2/m2/exports.m2 +++ b/M2/Macaulay2/m2/exports.m2 @@ -1009,6 +1009,7 @@ export { "poincare", "poincareN", "polarize", + "polylog", "position", "positions", "power", diff --git a/M2/Macaulay2/m2/typicalvalues.m2 b/M2/Macaulay2/m2/typicalvalues.m2 index 76256290af4..576033c16b8 100644 --- a/M2/Macaulay2/m2/typicalvalues.m2 +++ b/M2/Macaulay2/m2/typicalvalues.m2 @@ -130,7 +130,7 @@ then generateTypicalValues(currentFileDirectory | "../d/") ----------------------------------------------------------------------------- -- numerical functions that will be wrapped -redefs := hashTable apply({acos, agm, asin, atan, atan2, Beta, cos, cosh, cot, coth, csc, csch, Digamma, eint, erf, erfc, exp, expm1, Gamma, inverseErf, inverseRegularizedBeta, inverseRegularizedGamma, log, log1p, regularizedBeta, regularizedGamma, sec, sech, sin, sinh, sqrt, tan, tanh, zeta}, +redefs := hashTable apply({acos, agm, asin, atan, atan2, Beta, cos, cosh, cot, coth, csc, csch, Digamma, eint, erf, erfc, exp, expm1, Gamma, inverseErf, inverseRegularizedBeta, inverseRegularizedGamma, log, log1p, polylog, regularizedBeta, regularizedGamma, sec, sech, sin, sinh, sqrt, tan, tanh, zeta}, f -> f => method()); variants := new MutableHashTable; diff --git a/M2/Macaulay2/tests/normal/numbers.m2 b/M2/Macaulay2/tests/normal/numbers.m2 index 940b511f767..6779ee9a0b4 100644 --- a/M2/Macaulay2/tests/normal/numbers.m2 +++ b/M2/Macaulay2/tests/normal/numbers.m2 @@ -1013,6 +1013,10 @@ assert small abs(BesselY(1, ii) + 0.565159103992484 - 0.383186043874565*ii) assert (abs(BesselY(ii, 1) + 0.476556612479964 + 1.50506915911039*ii) < 1e-14) assert small abs(BesselY(ii, ii) + 0.665181892391867 - 0.395137431337008*ii) +assert small(polylog(1, 1/2) - log 2) +assert small(polylog(2, 1/2) - 1/12*pi^2 + 1/2*(log 2)^2) +assert(polylog(3, 1) == zeta 3) + assert Equation(truncate 1, 1) assert Equation(truncate 1.9, 1) assert Equation(truncate(-1.9), -1)