Merge pull request #285 from jverzani/N_speedup

Author: jverzani

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

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

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)