diff --git a/src/mathops.jl b/src/mathops.jl index 783335a..da9815c 100644 --- a/src/mathops.jl +++ b/src/mathops.jl @@ -83,7 +83,7 @@ Base.one(::Type{T}) where {T<:BasicType} = BasicType(Basic(1)) ## Math constants ## no oo! -for op in [:IM, :PI, :E, :EulerGamma, :Catalan, :oo, :zoo, :NAN] +for op in [:IM, :PI, :E, :EulerGamma, :Catalan, :GoldenRatio, :oo, :zoo, :NAN] @eval begin const $op = Basic(C_NULL) end @@ -108,6 +108,7 @@ function init_constants() @init_constant E E @init_constant EulerGamma EulerGamma @init_constant Catalan Catalan + @init_constant GoldenRatio GoldenRatio @init_constant oo infinity @init_constant zoo complex_infinity @init_constant NAN nan diff --git a/src/numerics.jl b/src/numerics.jl index 017ade1..52734fa 100644 --- a/src/numerics.jl +++ b/src/numerics.jl @@ -4,7 +4,7 @@ import Base: trunc, ceil, floor, round function evalf(b::Basic, bits::Integer=53, real::Bool=false) - !isfinite(b) && return b + (b == oo || b == zoo || b == NAN) && return b c = Basic() bits > 53 && real && (have_mpfr || throw(ArgumentError("libsymengine has to be compiled with MPFR for this feature"))) bits > 53 && !real && (have_mpc || throw(ArgumentError("libsymengine has to be compiled with MPC for this feature"))) @@ -16,89 +16,42 @@ function evalf(b::Basic, bits::Integer=53, real::Bool=false) end end -## Conversions from SymEngine -> Julia at the ccall level -function convert(::Type{BigInt}, b::BasicType{Val{:Integer}}) +## Conversions from SymEngine.Basic -> Julia at the ccall level +function _convert(::Type{BigInt}, b::Basic) a = BigInt() - c = Basic(b) - ccall((:integer_get_mpz, libsymengine), Nothing, (Ref{BigInt}, Ref{Basic}), a, c) + _convert_bigint!(a, b) return a end +function _convert_bigint!(a::BigInt, b::Basic) # non-allocating (sometimes) + is_a_Integer(b) || throw(ArgumentError("Not an integer")) + ccall((:integer_get_mpz, libsymengine), Nothing, (Ref{BigInt}, Ref{Basic}), a, b) + a +end -function convert(::Type{BigFloat}, b::BasicType{Val{:RealMPFR}}) - c = Basic(b) +function _convert(::Type{Int}, b::Basic) + is_a_Integer(b) || throw(ArgumentError("Not an integer")) + ccall((:integer_get_si, libsymengine), Int, (Ref{Basic},), b) +end + +function _convert(::Type{BigFloat}, b::Basic) a = BigFloat() - ccall((:real_mpfr_get, libsymengine), Nothing, (Ref{BigFloat}, Ref{Basic}), a, c) + _convert_bigfloat!(a, b) return a end -function convert(::Type{Cdouble}, b::BasicType{Val{:RealDouble}}) - c = Basic(b) - return ccall((:real_double_get_d, libsymengine), Cdouble, (Ref{Basic},), c) +function _convert_bigfloat!(a::BigFloat, b::Basic) # non-allocating + is_a_RealMPFR(b) || throw("Not a big value") + ccall((:real_mpfr_get, libsymengine), Nothing, (Ref{BigFloat}, Ref{Basic}), a, b) + a end -if SymEngine.libversion >= VersionNumber("0.4.0") - - function real(b::BasicComplexNumber) - c = Basic(b) - a = Basic() - ccall((:complex_base_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function imag(b::BasicComplexNumber) - c = Basic(b) - a = Basic() - ccall((:complex_base_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - -else - - function real(b::BasicType{Val{:ComplexDouble}}) - c = Basic(b) - a = Basic() - ccall((:complex_double_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function imag(b::BasicType{Val{:ComplexDouble}}) - c = Basic(b) - a = Basic() - ccall((:complex_double_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function real(b::BasicType{Val{:Complex}}) - c = Basic(b) - a = Basic() - ccall((:complex_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function imag(b::BasicType{Val{:Complex}}) - c = Basic(b) - a = Basic() - ccall((:complex_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function real(b::BasicType{Val{:ComplexMPC}}) - c = Basic(b) - a = Basic() - ccall((:complex_mpc_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - - function imag(b::BasicType{Val{:ComplexMPC}}) - c = Basic(b) - a = Basic() - ccall((:complex_mpc_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c) - return a - end - +function _convert(::Type{Cdouble}, b::Basic) + is_a_RealDouble(b) || throw(ArgumentError("Not a real double")) + return ccall((:real_double_get_d, libsymengine), Cdouble, (Ref{Basic},), b) end + ################################################## # N """ @@ -109,10 +62,11 @@ Convert a SymEngine numeric value into a Julian number N(a::Integer) = a N(a::Rational) = a N(a::Complex) = a -N(b::Basic) = N(BasicType(b)) -function N(b::BasicType{Val{:Integer}}) - a = convert(BigInt, b) +N(b::Basic) = N(get_symengine_class_val(b), b) + +function N(::Val{:Integer}, b::Basic) + a = _convert(BigInt, b) if (a.size > 1 || a.size < -1) return a elseif (a.size == 0) @@ -130,35 +84,92 @@ function N(b::BasicType{Val{:Integer}}) end end -N(b::BasicType{Val{:Rational}}) = Rational(N(numerator(b)), N(denominator(b))) # TODO: conditionally wrap rational_get_mpq from cwrapper.h -N(b::BasicType{Val{:RealDouble}}) = convert(Cdouble, b) -N(b::BasicType{Val{:RealMPFR}}) = convert(BigFloat, b) -N(b::BasicType{Val{:NaN}}) = NaN -function N(b::BasicType{Val{:Infty}}) - b == oo && return Inf - b == -oo && return -Inf - b == zoo && return Complex(Inf, Inf) +# TODO: conditionally wrap rational_get_mpq from cwrapper.h +N(::Val{:Rational}, b::Basic) = Rational(N(numerator(b)), N(denominator(b))) +N(::Val{:RealDouble}, b::Basic) = _convert(Cdouble, b) +N(::Val{:RealMPFR}, b::Basic) = _convert(BigFloat, b) +N(::Val{:Complex}, b::Basic) = complex(N(real(b)), N(imag(b))) +N(::Val{:ComplexMPC}, b::Basic) = complex(N(real(b)), N(imag(b))) +N(::Val{:ComplexDouble}, b::Basic) = complex(N(real(b)), N(imag(b))) + +N(::Val{:NaN}, b::Basic) = NaN +function N(::Val{:Infty}, b::Basic) + if b == oo + return Inf + elseif b == zoo + return Complex(Inf,Inf) + elseif b == -oo + return -Inf + else + throw(ArgumentError("Unknown infinity symbol")) + end end -## Mapping of SymEngine Constants into julia values -constant_map = Dict("pi" => π, "eulergamma" => γ, "exp(1)" => e, "catalan" => catalan, - "goldenratio" => φ) - -N(b::BasicType{Val{:Constant}}) = constant_map[toString(b)] +function N(::Val{:Constant}, b::Basic) + if b == PI + return π + elseif b == EulerGamma + return γ + elseif b == E + return ℯ + elseif b == Catalan + return catalan + elseif b == GoldenRatio + return φ + else + throw(ArgumentError("Unknown constant")) + end +end -N(b::BasicComplexNumber) = complex(N(real(b)), N(imag(b))) -function N(b::BasicType) - b = convert(Basic, b) - fs = free_symbols(b) - if length(fs) > 0 +function N(::Val{<:Any}, b::Basic) + is_constant(b) || throw(ArgumentError("Object can have no free symbols")) - end out = evalf(b) - imag(out) == Basic(0.0) ? real(out) : out + imag(out) == Basic(0.0) ? N(real(out)) : N(out) end +## deprecate N(::BasicType) +N(b::BasicType{T}) where {T} = N(convert(Basic, b), T) -## Conversions SymEngine -> Julia +## define convert(T, x) methods leveraging N() when needed +function convert(::Type{Float64}, x::Basic) + is_a_RealDouble(x) && return _convert(Cdouble, x) + convert(Float64, N(evalf(x, 53, true))) +end + +function convert(::Type{BigFloat}, x::Basic) + is_a_RealMPFR(x) && return _convert(BigFloat, x) + convert(BigFloat, N(evalf(x, precision(BigFloat), true))) +end + +function convert(::Type{Complex{Float64}}, x::Basic) + z = is_a_ComplexDouble(x) ? x : evalf(x, 53, false) + a,b = _real(z), _imag(z) + u,v = _convert(Cdouble, a), _convert(Cdouble, b) + return complex(u,v) +end + +function convert(::Type{Complex{BigFloat}}, x::Basic) + z = is_a_ComplexMPC(x) ? x : evalf(x, precision(BigFloat), false) + a,b = _real(z), _imag(z) + u,v = _convert(BigFloat, a), _convert(BigFloat, b) + return complex(u,v) +end + +convert(::Type{Number}, x::Basic) = x +convert(::Type{T}, x::Basic) where {T <: Real} = convert(T, N(x)) +convert(::Type{Complex{T}}, x::Basic) where {T <: Real} = convert(Complex{T}, N(x)) + +# Constructors no longer fall back to `convert` methods +Base.Int64(x::Basic) = convert(Int64, x) +Base.Int32(x::Basic) = convert(Int32, x) +Base.Float32(x::Basic) = convert(Float32, x) +Base.Float64(x::Basic) = convert(Float64, x) +Base.BigInt(x::Basic) = convert(BigInt, x) +Base.Real(x::Basic) = convert(Real, x) + + +## Rational -- p/q parts function as_numer_denom(x::Basic) a, b = Basic(), Basic() ccall((:basic_as_numer_denom, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}, Ref{Basic}), a, b, x) @@ -170,6 +181,28 @@ denominator(x::SymbolicType) = as_numer_denom(x)[2] numerator(x::SymbolicType) = as_numer_denom(x)[1] ## Complex +# b::Basic -> a::Basic +function _real(b::Basic) + if is_a_RealDouble(b) || is_a_RealMPFR(b) || is_a_Integer(b) || is_a_Rational(b) + return b + end + if !(is_a_Complex(b) || is_a_ComplexDouble(b) || is_a_ComplexMPC(b)) + throw(ArgumentError("Not a complex number")) + end + a = Basic() + ccall((:complex_base_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, b) + return a +end + +function _imag(b::Basic) + if !(is_a_Complex(b) || is_a_ComplexDouble(b) || is_a_ComplexMPC(b)) + throw(ArgumentError("Not a complex number")) + end + a = Basic() + ccall((:complex_base_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, b) + return a +end + real(x::Basic) = Basic(real(SymEngine.BasicType(x))) real(x::SymEngine.BasicType) = x @@ -186,22 +219,6 @@ conj(x::Basic) = Basic(conj(SymEngine.BasicType(x))) # To allow future extension, we define the fallback on `BasicType``. conj(x::BasicType) = 2 * real(x.x) - x.x -## define convert(T, x) methods leveraging N() -convert(::Type{Float64}, x::Basic) = convert(Float64, N(evalf(x, 53, true))) -convert(::Type{BigFloat}, x::Basic) = convert(BigFloat, N(evalf(x, precision(BigFloat), true))) -convert(::Type{Complex{Float64}}, x::Basic) = convert(Complex{Float64}, N(evalf(x, 53, false))) -convert(::Type{Complex{BigFloat}}, x::Basic) = convert(Complex{BigFloat}, N(evalf(x, precision(BigFloat), false))) -convert(::Type{Number}, x::Basic) = x -convert(::Type{T}, x::Basic) where {T <: Real} = convert(T, N(x)) -convert(::Type{Complex{T}}, x::Basic) where {T <: Real} = convert(Complex{T}, N(x)) - -# Constructors no longer fall back to `convert` methods -Base.Int64(x::Basic) = convert(Int64, x) -Base.Int32(x::Basic) = convert(Int32, x) -Base.Float32(x::Basic) = convert(Float32, x) -Base.Float64(x::Basic) = convert(Float64, x) -Base.BigInt(x::Basic) = convert(BigInt, x) -Base.Real(x::Basic) = convert(Real, x) ## For generic programming in Julia float(x::Basic) = float(N(x)) @@ -265,7 +282,9 @@ trunc(::Type{T},x::Basic, args...) where {T <: Integer} = convert(T, trunc(x,arg round(x::Basic; kwargs...) = Basic(round(N(x); kwargs...)) round(::Type{T},x::Basic; kwargs...) where {T <: Integer} = convert(T, round(x; kwargs...)) +prec(x::Basic) = prec(BasicType(x)) prec(x::BasicType{Val{:RealMPFR}}) = ccall((:real_mpfr_get_prec, libsymengine), Clong, (Ref{Basic},), x) +prec(::BasicType) = throw(ArgumentError("Method not applicable")) # eps eps(x::Basic) = eps(BasicType(x)) @@ -276,3 +295,26 @@ eps(::Type{BasicType{Val{:RealDouble}}}) = 2^-52 eps(::Type{BasicType{Val{:ComplexDouble}}}) = 2^-52 eps(x::BasicType{Val{:RealMPFR}}) = evalf(Basic(2), prec(x), true) ^ (-prec(x)+1) eps(x::BasicType{Val{:ComplexMPFR}}) = eps(real(x)) + +## convert from BasicType +function convert(::Type{BigInt}, b::BasicType{Val{:Integer}}) + _convert(BigInt, Basic(b)) +end + +function convert(::Type{BigFloat}, b::BasicType{Val{:RealMPFR}}) + _convert(BigInt, Basic(b)) +end + +function convert(::Type{Cdouble}, b::BasicType{Val{:RealDouble}}) + _convert(Cdouble, Basic(b)) +end + +## real/imag for BasicType +function real(b::BasicComplexNumber) + _real(Basic(b)) +end + +function imag(b::BasicComplexNumber) + _imag(Basic(b)) +end +## end deprecate diff --git a/src/types.jl b/src/types.jl index 64b9f1d..0d34bb3 100644 --- a/src/types.jl +++ b/src/types.jl @@ -122,9 +122,13 @@ function _get_symengine_classes() end const symengine_classes = _get_symengine_classes() +const symengine_classes_val = [Val(c) for c in SymEngine.symengine_classes] +const symengine_classes_val_type = [Val{c} for c in SymEngine.symengine_classes] "Get SymEngine class of an object (e.g. 1=>:Integer, 1//2 =:Rational, sin(x) => :Sin, ..." get_symengine_class(s::Basic) = symengine_classes[get_type(s) + 1] +get_symengine_class_val(s::Basic) = symengine_classes_val[get_type(s) + 1] +get_symengine_class_val_type(s::Basic) = symengine_classes_val_type[get_type(s) + 1] ## Construct symbolic objects @@ -221,8 +225,9 @@ SymbolicType = Union{Basic, BasicType} convert(::Type{Basic}, x::BasicType) = x.x Basic(x::BasicType) = x.x -BasicType(val::Basic) = BasicType{Val{get_symengine_class(val)}}(val) -convert(::Type{BasicType{T}}, val::Basic) where {T} = BasicType{Val{get_symengine_class(val)}}(val) +BasicType(val::Basic) = BasicType{get_symengine_class_val_type(val)}(val) +convert(::Type{BasicType{T}}, val::Basic) where {T} = + BasicType{get_symengine_class_val_type(val)}(val) # Needed for julia v0.4.7 convert(::Type{T}, x::Basic) where {T<:BasicType} = BasicType(x) @@ -263,6 +268,14 @@ BasicTrigFunction = Union{[SymEngine.BasicType{Val{i}} for i in trig_types]...} ### + +"Is expression constant" +function is_constant(ex::Basic) + syms = CSetBasic() + ccall((:basic_free_symbols, libsymengine), Nothing, (Ref{Basic}, Ptr{Cvoid}), ex, syms.ptr) + Base.length(syms) == 0 +end + "Is expression a symbol" function is_symbol(x::SymbolicType) res = ccall((:is_a_Symbol, libsymengine), Cuint, (Ref{Basic},), x) diff --git a/test/runtests.jl b/test/runtests.jl index a3ba7b8..aebb04e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -182,7 +182,7 @@ A = [x 2] @test lambdify(A)(1) == [1 2] @test isa(convert.(Expr, [0 x x+1]), Array{Expr}) -## N +## N, convert, _convert for val in samples @test N(Basic(val)) == val end @@ -191,6 +191,58 @@ for val in [π, γ, e, φ, catalan] @test N(Basic(val)) == val end +for a in Basic.((1, big(1))) + @test SymEngine.isinteger(a) + @test N(a) isa Int + @test (@allocated convert(Int, a)) > 0 + @test (@allocated SymEngine._convert(Int, a)) == 0 +end + +for a in Basic.((big(10)^100,)) + @test SymEngine.isinteger(a) + @test N(a) isa BigInt + @test (@allocated convert(BigInt, a)) > 0 + @test (@allocated SymEngine._convert(BigInt, a)) > 0 + + a′ = Basic(big(10)^10) + _b = BigInt() + @test (@allocated SymEngine._convert_bigint!(_b, a′)) == 0 + +end + + + +for a in Basic.((1.0, 1.23, Inf)) + @test N(a) isa Float64 + @test float(a) == N(a) + @test (@allocated convert(Float64, a)) >= 0 + @test (@allocated SymEngine._convert(Float64, a)) == 0 +end + +for a in Basic.((big(1.2)^100,)) + @test N(a) isa BigFloat + @test (@allocated convert(BigFloat, a)) > 0 + @test (@allocated SymEngine._convert(BigFloat, a)) > 0 + _b = BigFloat() + @test (@allocated SymEngine._convert_bigfloat!(_b, a)) == 0 +end + +for a in Basic.((1//2,)) + @test SymEngine.is_a_Rational(a) + @test N(a) isa Rational + @test (@allocated convert(Rational{Int}, a)) > 0 +end + +for (a,b) in (PI=>π, E=>ℯ, + GoldenRatio => Base.MathConstants.golden, + Catalan => Base.MathConstants.catalan, + oo=>Inf, zoo => complex(Inf, Inf), + ) + @test N(a) == b +end + +@test isnan(N(NAN)) + @test !isnan(x) @test isnan(Basic(0)/0)