diff --git a/.travis.yml b/.travis.yml index 4320596..70fe265 100644 --- a/.travis.yml +++ b/.travis.yml @@ -14,8 +14,6 @@ julia: notifications: email: false -coveralls: true - before_install: - sudo apt-add-repository -y ppa:leo.robol/mpsolve - sudo apt-get update -qq diff --git a/src/MPSolve.jl b/src/MPSolve.jl index 5b016f4..6aac3a2 100644 --- a/src/MPSolve.jl +++ b/src/MPSolve.jl @@ -161,7 +161,16 @@ function get_degree(context::MContext) ccall((:mps_context_get_degree, :libmps), Cint, (Ref{Cvoid},), context.cntxt_ptr) end -function get_roots_bigfloat(context::MContext) +function get_roots(context::MContext, output_precision::Int) + if output_precision <= 53 + roots,radii = get_roots_f64(context) + else + roots,radii = get_roots_big(context) + end + return roots, radii +end + +function get_roots_big(context::MContext) degree = get_degree(context) roots_m = Array{Mpsc}(undef, degree) for n=1:degree @@ -184,7 +193,7 @@ function get_roots_bigfloat(context::MContext) return (roots, radii) end -function get_roots_float64(context::MContext) +function get_roots_f64(context::MContext) degree = get_degree(context) roots_c = Array{Cplx}(undef, degree) radii = Array{Cdouble}(undef, degree) @@ -204,9 +213,9 @@ end """ - (approximations, radii) = mps_roots(coefficients, output_precision) + (approximations, radii) = mps_roots(coefficients, output_precision=53) -Approximate the roots of the polynomial specified by the array +Approximate the roots of a polynomial specified by the array of its coefficients. Output precision is specified in bits. # Example @@ -221,14 +230,33 @@ julia> all(map(x->abs((x^N - 1)/(N*x^(N - 1))), app) < rad) true ``` """ -function mps_roots(coeffs, output_precision::Integer=53) - context = MContext(coeffs) +function mps_roots(coefficients::AbstractVector, output_precision::Integer=53) + return mps_roots(coefficients; output_precision=output_precision) +end + +""" + (approximations, radii) = mps_roots(A, B, output_precision=53) + +Approximate the roots of a secular equation. Output precision is specified in bits. + +# Example +```jldoctest +julia> N = 64; + +julia> + +julia> + +``` +""" +function mps_roots(A::AbstractVector, B::AbstractVector, output_precision::Integer=53) + return mps_roots(A, B; output_precision=output_precision) +end + +function mps_roots(coeffs...; output_precision::Integer=53) + context = MContext(coeffs...) solve_poly(context, output_precision) - if output_precision <= 53 - (roots, radii) = get_roots_float64(context) - else - (roots, radii) = get_roots_bigfloat(context) - end + (roots, radii) = get_roots(context, output_precision) free_context(context) (roots, radii) end diff --git a/src/MpsNumberTypes.jl b/src/MpsNumberTypes.jl index 3af8bf6..91d1fe9 100644 --- a/src/MpsNumberTypes.jl +++ b/src/MpsNumberTypes.jl @@ -65,38 +65,54 @@ struct Mpf _mp_d::Ptr{Limb} function Mpf(b::BigFloat) + f = Ref{Mpf}() + ccall( + (:__gmpf_init_set, :libgmp), + Cvoid, + (Ref{Mpf}, Ref{BigFloat}), + f, + b + ) + end + + function Mpf(b::BigInt) f = Ref{Mpf}() ccall((:__gmpf_init, :libgmp), Cvoid, (Ref{Mpf},), f) ccall( - (:mpfr_get_f, :libmpfr), - Cint, (Ref{Mpf},Ref{BigFloat},MPFRRoundingMode), - f, b, MPFRRoundNearest + (:__gmpf_set_z, :libgmp), + Cvoid, (Ref{Mpf}, Ref{BigInt}), + f, + b ) return f[] end - + function Mpf(d::Float64) f = Ref{Mpf}() ccall((:__gmpf_init_set_d, :libgmp), Cvoid, (Ref{Mpf},Cdouble), f, d) return f[] end - function Mpf(i::Int64) + function Mpf(i::Int32) f = Ref{Mpf}() ccall((:__gmpf_init_set_si, :libgmp), Cvoid, (Ref{Mpf}, Clong), f, i) return f[] end end -Base.convert(::Type{Mpf}, x::T) where T<:Union{Int32,Int16,Int8} = Mpf(Int64(x)) +Base.convert(::Type{Mpf}, x::Signed) = Mpf(BigInt(x)) +Base.convert(::Type{Mpf}, x::T) where T<:Union{Int32,Int16,Int8} = Mpf(Int32(x)) Base.convert(::Type{Mpf}, x::AbstractFloat) = Mpf(x) function BigFloat(f::Mpf) b = BigFloat() ccall( (:mpfr_set_f, :libmpfr), - Cint, (Ref{BigFloat},Ref{Mpf},MPFRRoundingMode), - b, f, MPFRRoundNearest + Cint, + (Ref{BigFloat},Ref{Mpf},MPFRRoundingMode), + b, + f, + MPFRRoundNearest ) return b end diff --git a/test/runtests.jl b/test/runtests.jl index 30b005f..6eb96fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,8 +25,8 @@ function roots2secular_coeffs(roots) return L.*F,D end -function solve_test(rts, args; output_prec=53) - (app, rad) = mps_roots(args, output_prec) +function solve_test(rts, args...; output_prec=53) + (app, rad) = mps_roots(args..., output_prec) T = real(eltype(app)) for i = 1:length(rts) (err, ind) = findmin(map(abs, app .- rts[i])) @@ -65,7 +65,7 @@ end function test_secular_roots_unity(n) E = unity_roots(n) A,B = roots2secular_coeffs(E) - solve_test(E, (A, B), output_prec=256) + solve_test(E, A, B, output_prec=256) end """ Test solving the Wilkinson polynomial