diff --git a/README.md b/README.md index 175cb6c..74bb173 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ -[![Build Status](https://travis-ci.org/kagalenko-m-b/MPSolve.jl.svg?branch=master)](https://travis-ci.org/kagalenko-m-b/MPSolve.jl) -[![Coverage Status](https://coveralls.io/repos/github/kagalenko-m-b/MPSolve.jl/badge.svg)](https://coveralls.io/github/kagalenko-m-b/MPSolve.jl) +[![Build Status](https://travis-ci.com/kagalenko-m-b/MPSolve.jl.svg?branch=master)](https://travis-ci.com/kagalenko-m-b/MPSolve.jl) +![Coverage Status](https://coveralls.io/repos/github/kagalenko-m-b/MPSolve.jl/badge.svg?branch=master)(https://coveralls.io/github/kagalenko-m-b/MPSolve.jl?branch=master) # MPSolve.jl This repository contains the Julia interface to the MPSolve diff --git a/src/MPSolve.jl b/src/MPSolve.jl index 6aac3a2..de332bf 100644 --- a/src/MPSolve.jl +++ b/src/MPSolve.jl @@ -4,7 +4,7 @@ include("MpsNumberTypes.jl") using .MpsNumberTypes -export mps_roots,lagrange_barycentric +export mps_roots,mps_barycentric_coeffs abstract type PolyType end struct Monomial <: PolyType end @@ -107,7 +107,7 @@ function set_coefficients!( end function set_input_poly(context::MContext) - ccall((:mps_context_set_input_poly, :libmps), Cvoid, + ccall((:mps_context_set_input_poly, :libmps), Cvoid, (Ref{Cvoid}, Ref{Cvoid}), context.cntxt_ptr, context.poly) end @@ -261,6 +261,55 @@ function mps_roots(coeffs...; output_precision::Integer=53) (roots, radii) end + +raw""" + D,L,F = barycentric_coeffs_mps(f::Function; M::Int=0, D::AbstractVector=[]) + +Coefficients of barycentric form of Lagrange interpolation of the form used by MPSolve + + f(z) -> + + / L[0] L[M] \ + -> | ----------*F[0] + ... + ----------*F[M] - 1 | + \ z - D[0] z - D[M] / + + / / L[0] L[M] \ + / | ---------- + ... + ---------- | + / \ z - D[0] z - D[M] / + + If D not specified, set them to the roots of unity: + + / im*pi*k \ + D[k] = exp | --------- | , F[k] = f(D[k]) + \ M / + + For polynomial f() of the degree M, f(z) === f(C,D,L,F,z), where + + f(C,D,L,F,z)= any(z.==D) ? C*F[findfirst(z.==D)] : C*(eta(D,L.*F,z) - 1)/eta(D,L,z) + + See Berrut, Trefethen, "Barycentric Lagrange Interpolation," SIAM Review, + Vol. 46, No. 3, pp. 501-517, 2004 +""" +function mps_barycentric_coeffs(f::Function; M::Int=0, D::AbstractVector=[]) + F0 = f(0) + T = promote_type(typeof(F0), Irrational) + if isempty(D) + D = exp.(2im*T(pi)*(1:M)/M) + L = D/M + elseif M == 0 + M == length(D) + L = [1/(prod(D[k] .- D[1:k - 1])*prod(D[k] .- D[k + 1:end])) for k in 1:M] + else + throw(ArgumentError("specify either M or D")) + end + F = f.(D) + C = F0*sum(L./D) - sum((L.*F)./D) #;println("C=$(C)") + C != 0 || error("function appears to be a polynomial of the degree less than $(M)") + F /= C + return D,L,F +end + + function free_poly(context::MContext{Monomial}) ccall( (:mps_monomial_poly_free, :libmps), Cvoid, diff --git a/src/MpsNumberTypes.jl b/src/MpsNumberTypes.jl index 91d1fe9..541bbcb 100644 --- a/src/MpsNumberTypes.jl +++ b/src/MpsNumberTypes.jl @@ -2,8 +2,18 @@ module MpsNumberTypes using Base.GMP: Limb using Base.MPFR: MPFRRoundingMode,MPFRRoundNearest - -export Mpq,Mpf,Mpsc,Rdpe,Cplx,mps_clear! +import Base: convert, big +export Mpz,Mpq,Mpf,Mpsc,Rdpe,Cplx,mpsf_precision,mpsf_setprecision,mps_clear! + +function __init__() + try + # set GMP floating types precision to current precsion of julia's BigFloat + mpsf_setprecision(precision(BigFloat)) + catch ex + Base.showerror_nostdio(ex, "WARNING: Error during initialization of MpsNumberTypes") + end + nothing +end struct Mpz alloc::Cint @@ -26,7 +36,9 @@ function BigInt(m::Mpz) ccall((:__gmpz_set, :libgmp), Cvoid, (Ref{BigInt},Ref{Mpz}), b, m) return b end - +big(::Type{Mpz}) = BigInt +big(v::Mpz) = BigInt(v) + struct Mpq num::Mpz den::Mpz @@ -49,15 +61,19 @@ struct Mpq return q[] end end +mps_clear!(m::Mpz) = ccall((:__gmpz_clear, :libgmp), Cvoid, (Ref{Mpz},), m) +mps_clear!(m::Mpq) = ccall((:__gmpq_clear, :libgmp), Cvoid, (Ref{Mpq},), m) + Mpq(f::AbstractFloat) = Mpq(big(f)) Base.convert(::Type{Mpq}, x::T) where T<:Union{Signed,Rational} = Mpq(x) Base.Rational(q::Mpq) = Rational(BigInt(q.num), BigInt(q.den)) (::Type{T})(q::Mpq) where T<: AbstractFloat = T(Rational(q)) -mps_clear!(m::Mpz) = ccall((:__gmpz_clear, :libgmp), Cvoid, (Ref{Mpz},), m) -mps_clear!(m::Mpq) = ccall((:__gmpq_clear, :libgmp), Cvoid, (Ref{Mpq},), m) +big(::Type{Mpq}) = BigFloat +big(v::Mpq) = BigFloat(v) # Arbitrary precision floating point type from gmp.h used by MPSolve # is different from Julia's BigFloat + struct Mpf _mp_prec::Cint _mp_size::Cint @@ -66,13 +82,17 @@ struct Mpf function Mpf(b::BigFloat) f = Ref{Mpf}() + ccall((:__gmpf_init, :libgmp), Cvoid, (Ref{Mpf},), f) ccall( - (:__gmpf_init_set, :libgmp), - Cvoid, - (Ref{Mpf}, Ref{BigFloat}), + (:mpfr_get_f, :libmpfr), + Cint, + (Ref{Mpf},Ref{BigFloat},MPFRRoundingMode), f, - b + b, + MPFRRoundNearest ) + + return f[] end function Mpf(b::BigInt) @@ -84,6 +104,7 @@ struct Mpf f, b ) + return f[] end @@ -92,10 +113,11 @@ struct Mpf ccall((:__gmpf_init_set_d, :libgmp), Cvoid, (Ref{Mpf},Cdouble), f, d) return f[] end - + function Mpf(i::Int32) f = Ref{Mpf}() ccall((:__gmpf_init_set_si, :libgmp), Cvoid, (Ref{Mpf}, Clong), f, i) + return f[] end end @@ -117,6 +139,12 @@ function BigFloat(f::Mpf) return b end +big(::Type{Mpf}) = BigFloat +big(v::Mpf) = BigFloat(v) + +mpsf_precision() = ccall((:__gmpf_get_default_prec, :libgmp), Culong, ()) +mpsf_setprecision(p) = ccall((:__gmpf_set_default_prec, :libgmp), Cvoid, (Culong,), p) + struct Mpsc r::Mpf i::Mpf diff --git a/test/runtests.jl b/test/runtests.jl index 6eb96fa..e1489bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ -# Testing for the MPSolve.jl package. +# Testing for the MPSolve.jl package. -using MPSolve +using MPSolve: Mpz,Mpq,Mpf,mps_clear!,mps_roots,mps_barycentric_coeffs using Test unity_roots(n) = [exp(j*2im*BigFloat(pi)/n) for j = 1:n] @@ -16,21 +16,21 @@ function roots2coeffs(roots) coeffs end -function roots2secular_coeffs(roots) - M = length(roots) - T = eltype(roots) - D = exp.(im*T(pi)*range(1, stop=2M-1, length=M)/T(M)) - L = -D/M - F = [-prod(d .- roots) for d in D] - return L.*F,D -end +# function roots2secular_coeffs(roots) +# M = length(roots) +# T = eltype(roots) +# D = exp.(im*T(pi)*range(1, stop=2M-1, length=M)/T(M)) +# L = -D/M +# F = [-prod(d .- roots) for d in D] +# return L.*F,D +# end -function solve_test(rts, args...; output_prec=53) +function solve_test(roots, 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])) - @test err <= max(rad[ind], sqrt(eps(T))) + for rt in roots + (err, ind) = findmin(map(abs, app .- rt)) + @test err <= max(rad[ind], 1e4*eps(T)) end end @@ -63,9 +63,10 @@ function test_roots_of_unity_bigfloat(n) end function test_secular_roots_unity(n) - E = unity_roots(n) - A,B = roots2secular_coeffs(E) - solve_test(E, A, B, output_prec=256) + roots = unity_roots(n) + f(x) = prod(x .- roots) + D,L,F = mps_barycentric_coeffs(f, M=n) + solve_test(roots, L.*F, D, output_prec=256) end """ Test solving the Wilkinson polynomial @@ -77,9 +78,10 @@ function test_wilkinson(n) end function test_secular_wilkinson(n) - roots = big.(range(1, stop=n)) - A,B = roots2secular_coeffs(roots) - solve_test(roots, A, B) + roots = range(-n, stop=n*big(1.0), length=n)/n + f(x) = prod(x .- roots) + D,L,F = mps_barycentric_coeffs(f, M=n) + solve_test(roots, L.*F, D) end """ Test if solving a polynomial with complex integer coefficients @@ -102,13 +104,31 @@ if VERSION < v"1.3" error("This module only works with Julia version 1.3 or greater") end -for N = 99:100 - test_roots_of_unity(N) - test_roots_of_unity_fp(N) - test_roots_of_unity_bigint(N) - test_roots_of_unity_bigfloat(N) - # test_wilkinson(N) - test_secular_roots_unity(N) +@testset "Types for interfacing with MPSolve" begin + v = [Mpz(), Mpz(Int32(1234567890)), convert(Mpz, big"123456789012345678901234567890")] + @test eltype(v) == Mpz + @test eltype(big.(v)) == BigInt + @test big.(v) == [0, 1234567890, 123456789012345678901234567890] + mps_clear!(v) + # + v = Mpq(big(pi)) + @test typeof(v) == Mpq && typeof(big(v)) == BigFloat && big(v) == big(pi) + mps_clear!(v) + # + v = Mpf(big(pi)) + @test typeof(v) == Mpf && typeof(big(v)) == BigFloat && big(v) == big(pi) +end + +@testset "MPSolve routines" begin + for n in 99:100 + test_roots_of_unity(n) + test_roots_of_unity_fp(n) + test_roots_of_unity_bigint(n) + test_roots_of_unity_bigfloat(n) + test_wilkinson(n) + test_secular_roots_unity(n) + test_secular_wilkinson(n) + end + test_complex_int() + test_complex_bigint() end -test_complex_int() -test_complex_bigint()