Skip to content

Commit

Permalink
Fixed the conversion of integers to Mpf()
Browse files Browse the repository at this point in the history
  • Loading branch information
kagalenko-m-b committed Nov 18, 2020
1 parent a5148a0 commit e3ec988
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 24 deletions.
2 changes: 0 additions & 2 deletions .travis.yml
Expand Up @@ -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
Expand Down
50 changes: 39 additions & 11 deletions src/MPSolve.jl
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down
32 changes: 24 additions & 8 deletions src/MpsNumberTypes.jl
Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions test/runtests.jl
Expand Up @@ -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]))
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit e3ec988

Please sign in to comment.