Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use FLINT ball arithmetic for special functions #3286

Merged
merged 8 commits into from
Aug 6, 2024
Merged

Use FLINT ball arithmetic for special functions #3286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
7cbba23
Add basic conversion to/from flint's ball arithmetic types
d-torrance Jun 5, 2024
6f64dce
Add support for special functions w/ complex args
d-torrance Jun 5, 2024
a1ae809
Add support for special functions w/ real interval args
d-torrance Jun 5, 2024
226fd40
Move some special functions from boost math -> flint arb
d-torrance Jun 6, 2024
5a79a78
Fix tests that are failing after boostmath -> flint arb switch
d-torrance Jun 6, 2024
d42ffbc
Add missing documentation keys for numerical functions
d-torrance Jul 17, 2024
2d9cef3
Fix bug getting precision for 3-argument numeric methods.
d-torrance Jul 19, 2024
78045b0
Add unit tests for special functions with added CC/RRi support
d-torrance Jul 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions M2/Macaulay2/d/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ set(DLIST
pthread0.d
stdiop0.d
gmp.d
ballarith.d
engine.dd
xml.d # removed unless WITH_XML
stdio0.d
Expand Down
1 change: 1 addition & 0 deletions M2/Macaulay2/d/Makefile.files.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ M2_DFILES += interrupts.d
M2_DFILES += pthread0.d
M2_DFILES += stdiop0.d
M2_DFILES += gmp.d
M2_DFILES += ballarith.d
M2_DFILES += engine.dd
ifeq (@XML@,yes)
M2_DFILES += xml.d
Expand Down
235 changes: 210 additions & 25 deletions M2/Macaulay2/d/actors3.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

use evaluate;
use actors;
use ballarith;

isOption(e:Expr):bool := (
when e
Expand Down Expand Up @@ -844,25 +845,72 @@ setupfun("expm1",expm1).Protected=false;
eint(e:Expr):Expr := (
when e
is x:RRcell do toExpr(eint(x.v)) -- # typical value: eint, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(eint(x.v)) -- # typical value: eint, RRi, RRi
is x:CCcell do toExpr(eint(x.v)) -- # typical value: eint, CC, CC
else WrongArgRRorRRiorCC()
);
setupfun("eint",eint).Protected=false;
Gamma(e:Expr):Expr := (
when e
is x:RRcell do toExpr(Gamma(x.v)) -- # typical value: Gamma, RR, RR
is x:RRicell do toExpr(Gamma(x.v)) -- # typical value: Gamma, RRi, RRi
is x:CCcell do toExpr(Gamma(x.v)) -- # typical value: Gamma, CC, CC
is a:Sequence do
if length(a) != 2 then WrongNumArgs(2)
else when a.0 is s:RRcell do
when a.1 is x:RRcell do toExpr(Gamma(s.v, x.v)) -- # typical value: Gamma, RR, RR, RR
else WrongArgRR(2)
else WrongArgRR(1)
else WrongArgRR()
if length(a) != 2 then WrongNumArgs(1,2)
else (
when a.0
is s:RRcell do (
when a.1
is x:RRcell do toExpr(Gamma( s.v , x.v)) -- # typical value: Gamma, RR, RR, RR
is x:RRicell do toExpr(Gamma(toRRi(s.v), x.v)) -- # typical value: Gamma, RR, RRi, RRi
is x:CCcell do toExpr(Gamma( toCC(s.v), x.v)) -- # typical value: Gamma, RR, CC, CC
else WrongArgRRorRRiorCC(2))
is s:RRicell do (
when a.1
is x:RRcell do toExpr(Gamma(s.v, toRRi(x.v))) -- # typical value: Gamma, RRi, RR, RRi
is x:RRicell do toExpr(Gamma(s.v, x.v )) -- # typical value: Gamma, RRi, RRi, RRi
else WrongArgRRorRRi(2))
is s:CCcell do (
when a.1
is x:RRcell do toExpr(Gamma(s.v, toCC(x.v))) -- # typical value: Gamma, CC, RR, CC
is x:CCcell do toExpr(Gamma(s.v, x.v)) -- # typical value: Gamma, CC, CC, CC
else WrongArgRRorCC(2))
else WrongArgRRorCC(2))
else WrongArgRRorRRiorCC()
);
setupfun("Gamma",Gamma).Protected=false;
regularizedGamma(e:Expr):Expr := (
when e
is a:Sequence do (
if length(a) != 2 then WrongNumArgs(2)
else (
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: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 WrongArgRRorRRiorCC(2))
is s:RRicell do (
when a.1
is x:RRcell do toExpr(regularizedGamma(s.v, toRRi(x.v))) -- # typical value: regularizedGamma, RRi, RR, RRi
is x:RRicell do toExpr(regularizedGamma(s.v, x.v)) -- # typical value: regularizedGamma, RRi, RRi, RRi
else WrongArgRRorRRi(2))
is s:CCcell do (
when a.1
is x:RRcell do toExpr(regularizedGamma(s.v, toCC(x.v))) -- # typical value: regularizedGamma, CC, RR, CC
is x:CCcell do toExpr(regularizedGamma(s.v, x.v)) -- # typical value: regularizedGamma, CC, CC, CC
else WrongArgRRorCC(2))
else WrongArgRRorCC(2)))
else WrongNumArgs(2));
setupfun("regularizedGamma",regularizedGamma).Protected=false;
Digamma(e:Expr):Expr := (
when e
is x:RRcell do toExpr(Digamma(x.v)) -- # typical value: Digamma, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(Digamma(x.v)) -- # typical value: Digamma, RRi, RRi
is x:CCcell do toExpr(Digamma(x.v)) -- # typical value: Digamma, CC, CC
else WrongArgRRorRRiorCC()
);
setupfun("Digamma",Digamma).Protected=false;
export lgamma(x:RR):Expr := (
Expand All @@ -873,38 +921,80 @@ export lgamma(x:RR):Expr := (
lgamma(e:Expr):Expr := (
when e
is x:RRcell do lgamma(x.v) -- # typical value: lgamma, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(lgamma(x.v)) -- # typical value: lgamma, RRi, RRi
is x:CCcell do toExpr(lgamma(x.v)) -- # typical value: lgamma, CC, CC
else WrongArgRRorRRiorCC()
);
setupfun("lgamma",lgamma);
zeta(e:Expr):Expr := (
when e
is x:RRcell do toExpr(zeta(x.v)) -- # typical value: zeta, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(zeta(x.v)) -- # typical value: zeta, RRi, RRi
is x:CCcell do toExpr(zeta(x.v)) -- # typical value: zeta, CC, CC
else WrongArgRRorRRiorCC()
);
setupfun("zeta",zeta).Protected=false;
erf(e:Expr):Expr := (
when e
is x:RRcell do toExpr(erf(x.v)) -- # typical value: erf, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(erf(x.v)) -- # typical value: erf, RRi, RRi
is x:CCcell do toExpr(erf(x.v)) -- # typical value: erf, CC, CC
else WrongArgRRorRRiorCC()
);
setupfun("erf",erf).Protected=false;
erfc(e:Expr):Expr := (
when e
is x:RRcell do toExpr(erfc(x.v)) -- # typical value: erfc, RR, RR
else WrongArgRR()
is x:RRicell do toExpr(erfc(x.v)) -- # typical value: erfc, RRi, RRi
is x:CCcell do toExpr(erfc(x.v)) -- # typical value: erfc, CC, CC
else WrongArgRRorRRiorCC()
);
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:RRicell do toExpr(inverseErf(x.v)) -- # typical value: inverseErf, RRi, RRi
else WrongArgRRorRRi());
setupfun("inverseErf",inverseErf).Protected=false;
BesselJ(n:long,x:RR):RR := (
if n == long(0) then j0(x)
else if n == long(1) then j1(x)
else jn(n,x));
BesselJ(e:Expr):Expr := (
when e is s:Sequence do (
when s.0 is n:ZZcell do if !isLong(n.v) then WrongArg(1,"a small integer") else (
when s.1
is x:RRcell do toExpr(BesselJ(toLong(n.v),x.v))
else WrongArgRR(2))
else WrongArgZZ(1))
when s.0
is n:ZZcell do (
if isLong(n.v) then (
when s.1
is x:RRcell do toExpr(BesselJ(toLong(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 WrongArgRRorRRiorCC(2))
else (
when s.1
is x:RRcell do toExpr(
midpointRR(BesselJ(toRRi(n.v), toRRi(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 WrongArgRRorRRiorCC(2)))
is n:RRcell do (
when s.1
is x:RRcell do toExpr(
midpointRR(BesselJ(toRRi(n.v), toRRi(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 WrongArgRRorRRiorCC(2))
is n:RRicell do (
when s.1
is x:RRcell do toExpr(BesselJ(n.v, toRRi(x.v)))
is x:RRicell do toExpr(BesselJ(n.v, x.v))
else WrongArgRRorRRi(2))
is n:CCcell do (
when s.1
is x:RRcell do toExpr(BesselJ(n.v, toCC(x.v)))
is x:CCcell do toExpr(BesselJ(n.v, x.v ))
else WrongArgRRorCC(2))
else WrongArg(1, "an integer, real number or interval, or complex number"))
else WrongNumArgs(2));
setupfun("BesselJ",BesselJ).Protected=false;
BesselY(n:long,x:RR):RR := (
Expand All @@ -913,11 +1003,39 @@ BesselY(n:long,x:RR):RR := (
else yn(n,x));
BesselY(e:Expr):Expr := (
when e is s:Sequence do (
when s.0 is n:ZZcell do if !isLong(n.v) then WrongArg(1,"a small integer") else (
when s.1
is x:RRcell do toExpr(BesselY(toLong(n.v),x.v))
else WrongArgRR(2))
else WrongArgZZ(1))
when s.0
is n:ZZcell do (
if isLong(n.v) then (
when s.1
is x:RRcell do toExpr(BesselY(toLong(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 WrongArgRRorRRiorCC(2))
else (
when s.1
is x:RRcell do toExpr(
midpointRR(BesselY(toRRi(n.v), toRRi(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 WrongArgRRorRRiorCC(2)))
is n:RRcell do (
when s.1
is x:RRcell do toExpr(
midpointRR(BesselY(toRRi(n.v), toRRi(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 WrongArgRRorRRiorCC(2))
is n:RRicell do (
when s.1
is x:RRcell do toExpr(BesselY(n.v, toRRi(x.v)))
is x:RRicell do toExpr(BesselY(n.v, x.v))
else WrongArgRRorRRi(2))
is n:CCcell do (
when s.1
is x:RRcell do toExpr(BesselY(n.v, toCC(x.v)))
is x:CCcell do toExpr(BesselY(n.v, x.v ))
else WrongArgRRorCC(2))
else WrongArg(1, "an integer, real number or interval, or complex number"))
else WrongNumArgs(2));
setupfun("BesselY",BesselY).Protected=false;
atan2(yy:Expr,xx:Expr):Expr := (
Expand Down Expand Up @@ -954,16 +1072,83 @@ Beta(yy:Expr,xx:Expr):Expr := (
when yy
is y:RRcell do (
when xx
is x:RRcell do toExpr(Beta(y.v,x.v)) -- # typical value: Beta, RR, RR, RR
else WrongArgRR(1))
else WrongArgRR(2)
is x:RRcell do toExpr(Beta(y.v, x.v)) -- # typical value: Beta, RR, RR, RR
is x:RRicell do toExpr(Beta(toRRi(y.v), x.v)) -- # typical value: Beta, RR, RRi, RRi
is x:CCcell do toExpr(Beta(toCC(y.v), x.v)) -- # typical value: Beta, RR, CC, CC
else WrongArgRRorRRiorCC(2))
is y:RRicell do (
when xx
is x:RRcell do toExpr(Beta(y.v, toRRi(x.v))) -- # typical value: Beta, RRi, RR, RRi
is x:RRicell do toExpr(Beta(y.v, x.v)) -- # typical value: Beta, RRi, RRi, RRi
else WrongArgRRorRRi(2))
is y:CCcell do (
when xx
is x:RRcell do toExpr(Beta(y.v, toCC(x.v))) -- # typical value: Beta, CC, RR, CC
is x:CCcell do toExpr(Beta(y.v, x.v)) -- # typical value: Beta, CC, CC, CC
else WrongArgRRorCC(2))
else WrongArgRRorRRiorCC(1)
);
Beta(e:Expr):Expr := (
when e is s:Sequence do if length(s) == 2 then Beta(s.0,s.1)
else WrongNumArgs(2)
else WrongNumArgs(2));
setupfun("Beta",Beta).Protected=false;

regularizedBeta(xx:Expr,yy:Expr,zz:Expr):Expr := (
when xx
is x:RRcell do (
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: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 WrongArgRRorRRiorCC(3))
is y:RRicell do (
when zz
is z:RRcell do toExpr(regularizedBeta(toRRi(x.v), y.v, toRRi(z.v))) -- # typical value: regularizedBeta, RR, RRi, RR, RRi
is z:RRicell do toExpr(regularizedBeta(toRRi(x.v), y.v, z.v)) -- # typical value: regularizedBeta, RR, RRi, RRi, RRi
else WrongArgRRorRRi(3))
is y:CCcell do (
when zz
is z:RRcell do toExpr(regularizedBeta(toCC(x.v), y.v, toCC(z.v))) -- # typical value: regularizedBeta, RR, CC, RR, CC
is z:CCcell do toExpr(regularizedBeta(toCC(x.v), y.v, z.v)) -- # typical value: regularizedBeta, RR, CC, CC, CC
else WrongArgRRorCC(3))
else WrongArgRRorRRiorCC(2))
is x:RRicell do (
when yy
is y:RRcell do (
when zz
is z:RRcell do toExpr(regularizedBeta(x.v, toRRi(y.v), toRRi(z.v))) -- # typical value: regularizedBeta, RRi, RR, RR, RRi
is z:RRicell do toExpr(regularizedBeta(x.v, toRRi(y.v), z.v)) -- # typical value: regularizedBeta, RRi, RR, RRi, RRi
else WrongArgRRorRRi(3))
is y:RRicell do (
when zz
is z:RRcell do toExpr(regularizedBeta((x.v), y.v, toRRi(z.v))) -- # typical value: regularizedBeta, RRi, RRi, RR, RRi
is z:RRicell do toExpr(regularizedBeta(x.v, y.v, z.v)) -- # typical value: regularizedBeta, RRi, RRi, RRi, RRi
else WrongArgRRorRRi(3))
else WrongArgRRorRRi(2))
is x:CCcell do (
when yy
is y:RRcell do (
when zz
is z:RRcell do toExpr(regularizedBeta(x.v, toCC(y.v), toCC(z.v))) -- # typical value: regularizedBeta, CC, RR, RR, CC
is z:CCcell do toExpr(regularizedBeta(x.v, toCC(y.v), z.v)) -- # typical value: regularizedBeta, CC, RR, CC, CC
else WrongArgRRorRRiorCC(3))
is y:CCcell do (
when zz
is z:RRcell do toExpr(regularizedBeta(x.v, y.v, toCC(z.v))) -- # typical value: regularizedBeta, CC, CC, RR, CC
is z:CCcell do toExpr(regularizedBeta(x.v, y.v, z.v)) -- # typical value: regularizedBeta, CC, CC, CC, CC
else WrongArgRRorCC(3))
else WrongArgRRorCC(2))
else WrongArgRRorRRiorCC(1));
regularizedBeta(e:Expr):Expr := (
when e is s:Sequence do if length(s) == 3 then regularizedBeta(s.0,s.1,s.2)
else WrongNumArgs(3)
else WrongNumArgs(3));
setupfun("regularizedBeta",regularizedBeta).Protected=false;

cosh(e:Expr):Expr := (
when e
is x:CCcell do toExpr(cosh(x.v)) -- # typical value: cosh, CC, CC
Expand Down
Loading
Loading