Better Float32 support. Populate JWSEXPR.

src/algebra/Makefile.in

@@ -61,8 +61,8 @@ SPAD_SRCS= \
      intaf intalg intaux intclos intden intef integer \
      integrat interval intfact intlocp intpar intpm intrf intspec irexpand \
      irsn ituple jarray32 jarray64 jcfsf jla32 jla64 jnemo jnpoly jobject \
-	 jplot jfsf jfsf2 jf64sf jf64sf2 jnball julia jutils jet jws jwsagg \
-	 jwsexpr jwsnsf jwsutils kl kovacic laplace laurent leadcdet lie \
+	 jplot jfsf jfsf2 jf32sf jf64sf jf32sf2 jf64sf2 jnball julia jutils jet \
+	 jws jwsagg jwsexpr jwsnsf jwsutils kl kovacic laplace laurent leadcdet lie \
 	 limitps lindep lingrob linpen liouv listgcd list lll lmdict lodof \
 	 lodof2 lodo logic mama manip mantepse math_sym mathml mappkg matcat \
 	 matfuns matrix matstor mesh mfinfact mkfunc mkrecord \
@@ -398,11 +398,11 @@ ifdef JULIA_WRAP_SO_TARGET
 JULIALIST = INECF JARBPR JCF32 JCF32VEC JCF32MAT JCF32SMA \
             JCF64 JCF64MAT JCF64MTF JCF64SMA JCF64VEC JCFSF JCFLOAT JCF32LA \
             JCF64LA JCRING JDFRAME JDRAW JF32 JF32AF JF32VEC JF32VEC2 JF32MAT \
-            JF32SMAT  JF64 JF64AF JF64MAT JF64MTF JF64SF JF64SF2 JF64SMAT \
-            JF64VEC JF64VEC2 JFLOAT JFSF JFSF2 JI64 JI64VEC JMATCAT JMATRIX \
-            JMFLOAT JMTYPE JOBJECT JOBAGG JOBBOOL JOBBINT JOBCF64 JOBDICT \
-            JOBDLINK JOBF64 JOBI64 JOBPAIR JOBRING JOBTPLE JOBTYPE JPLOT \
-            JRING JF32LA JF64LA JSTRU JSYM JTYPE JVECCAT JVECTOR JVECTOR2 \
+            JF32SMAT  JF64 JF64AF JF64MAT JF64MTF JF32SF JF64SF JF32SF2 JF64SF2 \
+			JF64SMAT JF64VEC JF64VEC2 JFLOAT JFSF JFSF2 JI64 JI64VEC JMATCAT \
+			JMATRIX JMFLOAT JMTYPE JOBJECT JOBAGG JOBBOOL JOBBINT JOBCF32 JOBCF64 \
+			JOBDICT JOBDLINK JOBF32 JOBF64 JOBI64 JOBPAIR JOBRING JOBTPLE JOBTYPE \
+			JPLOT JRING JF32LA JF64LA JSTRU JSYM JTYPE JVECCAT JVECTOR JVECTOR2 \
             NACF NACFS JUF NAN NACB NARB NCB NCF NCRING NECF NFF NFR NFRAC \
             NINT NMLP NMP NPF NRB NRF NRING NTYPE NULP NUP NZMOD
 JULIACATLIST = JARBPR JCRING JMATCAT JMFLOAT JMTYPE JOBAGG JOBRING JOBTYPE \

src/algebra/exposed.lsp

@@ -163,6 +163,7 @@
   (|JuliaFloat32| . JF32)
   (|JuliaFloat64| . JF64)
   (|JuliaFloat64Matrix| . JF64MAT)
+  (|JuliaFloat32SpecialFunctions| . JF32SF)
   (|JuliaFloat64SpecialFunctions| . JF64SF)
   (|JuliaFloatSpecialFunctions| . JFSF)
   (|JuliaFloat64Vector| . JF64VEC)
@@ -176,7 +177,7 @@
   (|JuliaObjBoolean| . JOBBOOL)
   (|JuliaObjDict| . JOBDICT)
   (|JuliaObjDynamicLinker| . JOBDLINK)
-  (|JuliaObjFloat64| . JOBF64)
+  (|JuliaObjFloat32| . JOBF32) (|JuliaObjFloat64| . JOBF64)
   (|JuliaObjInt64| . JOBI64)
   (|JuliaObjPair| . JOBPAIR)
   (|JuliaObjTuple| . JOBTPLE)
@@ -810,8 +811,10 @@
   (|JuliaFloat32SquareMatrix| . JF32SMAT)
   (|JuliaFloat32Vector| . JF32VEC)
   (|JuliaFloat32VectorFunctions2| . JF32VEC2)
+  (|JuliaFloat32SpecialFunctions2| . JF32SF2)
   (|JuliaFloat64SpecialFunctions2| . JF64SF2)
   (|JuliaFloatSpecialFunctions2| . JFSF2)
+  (|JuliaObjComplexF32| . JOBCF32)
   (|JuliaObjComplexF64| . JOBCF64)
   (|Kernel| . KERNEL)
   (|Kovacic| . KOVACIC)

src/algebra/jf32sf.spad

@@ -0,0 +1,94 @@
+)abbrev package JF32SF JuliaFloat32SpecialFunctions
+++ Special functions computed using Julia's ecosystem.
+++ They are here essentially for "completeness" purpose with
+++ JuliaFloat32. You should use the DoubleFloat's special
+++ functions if available, calling Julia functions is costly.
+JuliaFloat32SpecialFunctions() : Exports == Implementation where
+  JF32 ==> JuliaFloat32
+  JI64 ==> JuliaInt64
+  Exports ==> with
+    ldexp : (JF32, JI64) -> JF32
+    ++ ldexp(x,n) computes x*2^n
+    exp2   : JF32 -> JF32
+    ++ exp2(x) computes the base 2 exponential of x.
+    exp10   : JF32 -> JF32
+    ++ exp10(x) computes the base 10 exponential of x.
+    log1p : JF32 -> JF32
+    ++ log1p(x) computes accurately natural logarithm of 1+x.
+
+    sind   : JF32 -> JF32
+    ++ sind(x) computes sine of x, where x is in degrees.
+    sinpi  : JF32 -> JF32
+    ++ sinpi(x) computes sin(pi*x) more accurately.
+    sinc   : JF32 -> JF32
+    ++ sinc(x) computes sin(pi*x)/(pi*x) if x~=0, and 1 if x=0.
+    cosd  : JF32 -> JF32
+    ++ cosd(x) computes cosine of x, where x is in degrees.
+    cospi : JF32 -> JF32
+    ++ cospi(x) computes cos(pi*x) more accurately.
+    cosc  : JF32 -> JF32
+    ++ cosc(x) computes cos(pi*x)/x−sin(pi*x)/(pi*x^2)
+    ++ if x~=0, and 0 if x=0 i.e. the derivative of sinc(x).
+    tand  : JF32 -> JF32
+    ++ tand(x) computes tangent of x, where x is in degrees.
+    tanpi : JF32 -> JF32
+    ++ tanpi(x) computes tan(pi*x) more accurately.
+    cotd  : JF32 -> JF32
+    ++ cotd(x) computes cotangent of x, where x is in degrees.
+    secd  : JF32 -> JF32
+    ++ secd(x) computes secant of x, where x is in degrees.
+    cscd  : JF32 -> JF32
+    ++ cscd(x) computes cosecant of x, where x is in degrees.
+    hypot : (JF32, JF32) -> JF32
+    ++ hypot(x,y) computes the hypotenuse avoiding overflow and underflow.
+
+    asind : JF32 -> JF32
+    ++ asind(x) computes the inverse sine of x, where output is in degrees.
+    acosd : JF32 -> JF32
+    ++ acosd(x) computes the inverse cosine of x, where output is in degrees.
+    atand : JF32 -> JF32
+    ++ atand(x) computes the inverse tangent of x, where output is in degrees.
+    atand  : (JF32, JF32) -> JF32
+    ++ atand(x, y) computes the inverse tangent of x/y, where output is in degrees.
+    acotd : JF32 -> JF32
+    ++ acotd(x) computes the inverse cotangent of x, where output is in degrees.
+    asecd : JF32 -> JF32
+    ++ asecd(x) computes the inverse secant of x, where output is in degrees.
+    acscd : JF32 -> JF32
+    ++ acscd(x) computes the inverse cosecant of x, where output is in degrees.
+    rad2deg  : JF32 -> JF32
+    ++ rad2deg(x) converts x to degrees, where x is in radians.
+    deg2rad  : JF32 -> JF32
+    ++ deg2rad(x) converts x to radian, where x is in degrees.
+
+  Implementation ==> add
+
+    ldexp(x,n) == jl_dbl_function_dbl_int64("ldexp", x, n)$Lisp
+    exp2(x)    == jlApply("exp2", x)
+    exp10(x)   == jlApply("exp10", x)
+    log1p(x)   == jlApply("log1p", x)
+
+    sinpi(x) == jlApply("sinpi", x)
+    sinc(x)  == jlApply("sinc", x)
+    cospi(x) == jlApply("cospi", x)   
+    cosc(x)  == jlApply("cosc", x)   
+    tanpi(x) == jlApply("tanpi", x)
+    hypot(x,y) == jlApply("hypot", x, y)
+
+    sind(x) == jlApply("sind", x)
+    cosd(x) == jlApply("cosd", x)
+    tand(x) == jlApply("tand", x)
+    cotd(x) == jlApply("cotd", x)
+    secd(x) == jlApply("secd", x)
+    cscd(x) == jlApply("cscd", x)
+
+    asind(x) == jlApply("asind", x)
+    acosd(x) == jlApply("acosd", x)
+    atand(x) == jlApply("atand", x)
+    atand(x,y) == jlApply("atand", x, y)
+    acotd(x) == jlApply("acotd", x)
+    asecd(x) == jlApply("asecd", x)
+    acscd(x) == jlApply("acscd", x)
+
+    rad2deg(x) == jlApply("rad2deg", x)
+    deg2rad(x) == jlApply("deg2rad", x)

src/algebra/jf32sf2.spad

@@ -0,0 +1,245 @@
+)abbrev package JF32SF2 JuliaFloat32SpecialFunctions2
+++ Special functions computed using Julia's ecosystem.
+++ They are here essentially for "completeness" purpose with
+++ JuliaFloat32. You should use the DoubleFloat's special
+++ functions if available, calling Julia functions is costly.
+JuliaFloat32SpecialFunctions2() : Exports == Implementation where
+  JF32 ==> JuliaFloat32
+  JI64 ==> JuliaInt64
+  Exports ==> with
+
+    -- Gamma Function
+
+    Gamma : JF32 -> JF32
+    ++ Gamma(z) computes Gamma function \Gamma(z)
+    logGamma : JF32 -> JF32
+    ++ logGamma(x) computes accurate log(gamma(x)) for large x
+    logabsgamma : JF32 -> JF32
+    ++ logabsgamma(x) computes accurate log(abs(gamma(x))) for large x
+    --logfactorial : JF32 -> JF32
+    --++ logfactorial(x) computes accurate log(factorial(x)) for large x; same as loggamma(x+1) for x > 1, zero otherwise
+    digamma : JF32 -> JF32
+    ++ digamma(x) computes digamma function (i.e. the derivative of loggamma at x)
+    invdigamma : JF32 -> JF32
+    ++ invdigamma(x) computes invdigamma function (i.e. inverse of digamma function at x using fixed-point iteration algorithm)
+    trigamma : JF32 -> JF32
+    ++ trigamma(x) computes trigamma function (i.e the logarithmic second derivative of gamma at x)
+    polygamma : (JI64, JF32) -> JF32
+    ++ polygamma(m,x) computes polygamma function (i.e the (m+1)-th derivative of the loggamma function at x)
+    Gamma : (JF32, JF32) -> JF32
+    ++ Gamma(a,z) computes upper incomplete gamma function \Gamma(a,z)
+    logGamma : (JF32, JF32) -> JF32
+    ++ logGamma(a,z) computes accurate log(gamma(a,x)) for large arguments
+    --gamma_inc : (JF32, JF32, JF32) -> JF32
+    --++ gamma_inc(a,x,IND) computes incomplete gamma function ratio P(a,x)
+    --++ and Q(a,x) (i.e evaluates P(a,x) and Q(a,x) for accuracy specified by IND and returns tuple (p,q))
+    --beta_inc(a,b,x,y) : JF32 -> JF32
+    --++ beta_inc(a,b,x,y) computes incomplete beta function ratio Ix(a,b) and Iy(a,b) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q))
+    gamma_inc_inv : (JF32, JF32, JF32) -> JF32
+    ++ gamma_inc_inv(a,p,q) computes inverse of incomplete gamma function ratio P(a,x) and Q(a,x) (i.e evaluates x given P(a,x)=p and Q(a,x)=q)
+    Beta : (JF32, JF32) -> JF32
+    ++ Beta(x,y) computes beta function at x,y
+    logBeta : (JF32, JF32) -> JF32
+    ++ logBeta(x,y) computes accurate log(beta(x,y)) for large x or y
+    logabsbeta : (JF32, JF32) -> JF32
+    ++ logabsbeta(x,y) computes accurate log(abs(beta(x,y))) for large x or y
+    --logabsbinomial : (JF32, JF32) -> JF32
+    --++ logabsbinomial(x,y) computes accurate log(abs(binomial(n,k))) for large n and k near n/2
+
+    -- Exponential and Trigonometric Integrals
+
+    expint : (JF32, JF32) -> JF32
+    ++ expint(nu, z) computes exponential integral function
+    Ei : JF32 -> JF32
+    ++ Ei(x) computes exponential integral Ei(x)
+    expintx: JF32 -> JF32
+    ++ expintx(x) computes scaled exponential integral function
+    Si : JF32 -> JF32
+    ++ Si(x) computes sine integral Si(x)
+    Ci : JF32 -> JF32
+    ++ Ci(x) computes cosine integral Ci(x)
+
+    -- Error Functions, Dawson’s and Fresnel Integrals
+
+    erf : JF32 -> JF32
+    ++ erf(x) computes error function at x
+    erf : (JF32, JF32) -> JF32
+    ++ erf(x,y) computes accurate version of erf(y) - erf(x)
+    erfc : JF32 -> JF32
+    ++ erfc(x) computes complementary error function, i.e. the accurate version of 1-erf(x) for large x
+    inverseErfc : JF32 -> JF32
+    ++ inverseErfc(x) computes inverse function of erfc.
+    erfcx : JF32 -> JF32
+    ++ erfcx(x) computes scaled complementary error function, i.e. accurate e^(x^2) erfc(x) for large x
+    logerfc : JF32 -> JF32
+    ++ logerfc(x) computes log of the complementary error function, i.e. accurate ln(erfc(x)) for large x
+    logerfcx : JF32 -> JF32
+    ++ logerfcx(x) computes log of the scaled complementary error function, i.e. accurate ln(erfcx(x)) for large negative x
+    erfi : JF32 -> JF32
+    ++ erfi(x) computes imaginary error function defined as -i erf(ix)
+    inverseErf : JF32 -> JF32
+    ++ inverseErf(x) computes inverse function of erf()
+    dawson : JF32 -> JF32
+    ++ dawson(x) computes scaled imaginary error function, a.k.a. Dawson function.
+
+    -- Airy and Related Functions
+
+    airyAi : JF32 -> JF32
+    ++ airyAi(z) computes Airy Ai function at z
+    airyAiPrime : JF32 -> JF32
+    ++ airyAiPrime(z) computes derivative of the Airy Ai function at z
+    airyBi : JF32 -> JF32
+    ++ airyBi(z) computes Airy Bi function at z
+    airyBiPrime : JF32 -> JF32
+    ++ airyBiPrime(z) computes derivative of the Airy Bi function at z
+    airyAix : JF32 -> JF32
+    ++ airyAix(z) computes scaled Airy Ai function and kth derivatives at z
+    airyAiPrimex : JF32 -> JF32
+    ++ airyAiPrimex(z) computes scaled derivative of the Airy Ai function at z
+    airyBix : JF32 -> JF32
+    ++ airyBix(z) computes scaled  Airy Bi function at z
+    airyBiPrimex : JF32 -> JF32
+    ++ airyBiPrimex(z) computes scaled derivative of the Airy Bi function at z
+
+    -- Bessel Functions
+
+    besselJ : (JF32, JF32) -> JF32
+    ++ besselJ(nu,z) computes Bessel function of the first kind of order nu at z
+    besselJ0 : JF32 -> JF32
+    ++ besselJ0(z) computes besselj(0,z)
+    besselJ1 : JF32 -> JF32
+    ++ besselJ1(z) computes besselj(1,z)
+    besselJx : (JF32, JF32) -> JF32
+    ++ besselJx(nu,z) computes scaled Bessel function of the first kind of order nu at z
+    sphericalBesselJ : (JF32, JF32) -> JF32
+    ++ sphericalBesselJ(nu,z) computes Spherical Bessel function of the first kind of order nu at z
+    besselY : (JF32, JF32) -> JF32
+    ++ besselY(nu,z) computes Bessel function of the second kind of order nu at z
+    besselY0 : JF32 -> JF32
+    ++ besselY0(z) computes bessely(0,z)
+    besselY1 : JF32 -> JF32
+    ++ besselY1(z) computes bessely(1,z)
+    besselYx : (JF32, JF32) -> JF32
+    ++ besselYx(nu,z) computes scaled Bessel function of the second kind of order nu at z
+    sphericalBesselY : (JF32, JF32) -> JF32
+    ++ sphericalBesselY(nu,z) computes Spherical Bessel function of the second kind of order nu at z
+    --besselh : (JF32, JF32, JF32) -> JF32
+    --++ besselh(nu,k,z) computes Bessel function of the third kind (a.k.a. Hankel function) of order nu at z; k must be either 1 or 2
+    hankelH1 : (JF32, JF32) -> JF32
+    ++ hankelH1(nu,z) computes besselh(nu, 1, z)
+    hankelH1x : (JF32, JF32) -> JF32
+    ++ hankelH1x(nu,z) computes scaled besselh(nu, 1, z)
+    hankelH2 : (JF32, JF32) -> JF32
+    ++ hankelH2(nu,z) computes besselh(nu, 2, z)
+    hankelH2x : (JF32, JF32) -> JF32
+    ++ hankelH2x(nu,z) computes scaled besselh(nu, 2, z)
+    besselI : (JF32, JF32) -> JF32
+    ++ besselI(nu,z) computes modified Bessel function of the first kind of order nu at z
+    besselIx : (JF32, JF32) -> JF32
+    ++ besselIx(nu,z) computes scaled modified Bessel function of the first kind of order nu at z
+    besselK : (JF32, JF32) -> JF32
+    ++ besselK(nu,z) computes modified Bessel function of the second kind of order nu at z
+    besselKx : (JF32, JF32) -> JF32
+    ++ besselKx(nu,z) computes scaled modified Bessel function of the second kind of order nu at z
+    jinc : JF32 -> JF32
+    ++ jinc(x) computes scaled Bessel function of the first kind divided by x. A.k.a. sombrero or besinc
+
+    -- Elliptic Integrals
+
+    ellipticK : JF32 -> JF32
+    ++ ellipticK(m) computes complete elliptic integral of 1st kind K(m)
+    ellipticE : JF32 -> JF32
+    ++ ellipticE(m) computes complete elliptic integral of 2nd kind E(m)
+
+    -- Zeta and Related Functions
+
+    eta : JF32 -> JF32
+    ++ eta(x) computes Dirichlet eta function at x
+    riemannZeta : JF32 -> JF32
+    ++ riemannZeta(x) computes Riemann zeta function at x
+
+  Implementation ==> add
+    -- Gamma Functions
+
+    Gamma(z) == jlApply("gamma", z)
+    logGamma(x) == jlApply("loggamma", x)
+    logabsgamma(x) == jlApply("logabsgamma", x)
+    --logfactorial(x) == jlApply("logfactorialx)
+    digamma(x) == jlApply("digamma", x)
+    invdigamma(x) == jlApply("invdigamma", x)
+    trigamma(x) == jlApply("trigamma", x)
+    polygamma(m,x) == jl_dbl_function_int64_dbl("polygamma", m, x)$Lisp
+    Gamma(a,z) == jlApply("gamma", a, z)
+    logGamma(a,z) == jlApply("loggamma", a, z)
+    --gamma_inc(a,x,IND) == jlApply("gamma_inc", a, x, IND)
+    --beta_inc(a,b,x,y) == jlApply("beta_inc", x)
+    gamma_inc_inv(a,p,q) == jlApply("gamma_inc_inv", a, p, q)
+    Beta(x,y) == jlApply("beta", x, y)
+    logBeta(x,y) == jlApply("logbeta", x, y)
+    logabsbeta(x,y) == jlApply("logabsbeta", x, y)
+    --logabsbinomial(x,y) == jlApply("logabsbinomial", x, y)
+
+    -- Exponential and Trigonometric Integrals
+
+    expint(nu, z) == jlApply("expint", nu, z)
+    Ei(x) == jlApply("expinti", x)
+    expintx(x) == jlApply("expintx", x)
+    Si(x) == jlApply("sinint", x)
+    Ci(x) == jlApply("cosint", x)
+
+    -- Error Functions, Dawson’s and Fresnel Integrals
+
+    erf(x) == jlApply("erf", x)
+    erf(x,y) == jlApply("erf", x, y)
+    erfc(x) == jlApply("erfc", x)
+    inverseErf(x) == jlApply("erfinv", x)
+    inverseErfc(x) == jlApply("erfcinv", x)
+    erfcx(x) == jlApply("erfcx", x)
+    logerfc(x) == jlApply("logerfc", x)
+    logerfcx(x) == jlApply("logerfcx", x)
+    erfi(x) == jlApply("erfi", x)
+    dawson(x) == jlApply("dawson", x)
+
+    -- Airy and Related Functions
+
+    airyAi(z) == jlApply("airyai", z)
+    airyAiPrime(z) == jlApply("airyaiprime", z)
+    airyBi(z) == jlApply("airybi", z)
+    airyBiPrime(z) == jlApply("airybiprime", z)
+    airyAix(z) == jlApply("airyaix", z)
+    airyAiPrimex(z) == jlApply("airyaiprimex", z)
+    airyBix(z) == jlApply("airybix", z)
+    airyBiPrimex(z) == jlApply("airybiprimex", z)
+
+    -- Bessel Functions
+
+    besselJ(nu,z) == jlApply("besselj", nu, z)
+    besselJ0(z) == jlApply("besselj0", z)
+    besselJ1(z) == jlApply("besselj1", z)
+    besselJx(nu,z) == jlApply("besseljx", nu, z)
+    sphericalBesselJ(nu,z) == jlApply("sphericalbesselj", nu, z)
+    besselY(nu,z) == jlApply("bessely", nu, z)
+    besselY0(z) == jlApply("bessely0", z)
+    besselY1(z) == jlApply("bessely1", z)
+    besselYx(nu,z) == jlApply("besselyx", nu, z)
+    sphericalBesselY(nu,z) == jlApply("sphericalbessely", nu, z)
+    --besselh(nu,k,z) == jlApply("besselh", nu, k, z)
+    hankelH1(nu,z) == jlApply("hankelh1", nu, z)
+    hankelH1x(nu,z) == jlApply("hankelh1x", nu, z)
+    hankelH2(nu,z) == jlApply("hankelh2", nu, z)
+    hankelH2x(nu,z) == jlApply("hankelh2x", nu, z)
+    besselI(nu,z) == jlApply("besseli", nu, z)
+    besselIx(nu,z) == jlApply("besselix", nu, z)
+    besselK(nu,z) == jlApply("besselk", nu, z)
+    besselKx(nu,z) == jlApply("besselkx", nu, z)
+    jinc(x) == jlApply("jinc", x)
+
+    -- Elliptic Integrals
+
+    ellipticK(m) == jlApply("ellipk", m)
+    ellipticE(m) == jlApply("ellipe", m)
+
+    -- Zeta and Related Functions
+
+    eta(x) == jlApply("eta", x)
+    riemannZeta(x) == jlApply("zeta", x)