Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updates for Julia 0.6 #62

Merged
merged 3 commits into from
Mar 20, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,13 @@ os:
- osx
- linux
julia:
- nightly
- 0.4
- 0.5
- nightly
notifications:
email: false
script:
- if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
- julia -e 'Pkg.clone(pwd()); Pkg.build("Roots"); Pkg.test("Roots"; coverage=true)';
#script:
# - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
# - julia -e 'Pkg.clone(pwd()); Pkg.build("Roots"); Pkg.test("Roots"; coverage=true)';
after_success:
- julia -e 'cd(Pkg.dir("Example")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';
- julia -e 'cd(Pkg.dir("Roots")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';
2 changes: 1 addition & 1 deletion REQUIRE
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,4 @@ julia 0.4
ForwardDiff 0.2.0
Polynomials 0.0.5
PolynomialFactors 0.0.3
Compat 0.8.0
Compat 0.17.0
1 change: 0 additions & 1 deletion src/Polys/agcd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,6 @@ function precondition(p::Poly,q::Poly)
phis = [2.0^i for i in -5:5]
out = (1,1)
m = ratio(p,q)

for α in alphas
for ϕ in phis
r = ratio(polyval(p, ϕ * variable(p)), α * polyval(q, ϕ * variable(q)))
Expand Down
39 changes: 18 additions & 21 deletions src/Polys/polynomials.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,5 @@
## Extensions to Polynomials.jl

" Return the variable `x` of a polynomial as a polynomial."
variable(p::Poly) = poly(zeros(eltype(p),1), p.var)

"Create a monic polynomial from `p`"
monic(p::Poly) = Poly(p.a/p[degree(p)], p.var)

Expand All @@ -20,25 +17,25 @@ end
Base.convert(::Type{Function}, p::Poly) = x -> Polynomials.polyval(p,x)
## convert a function to a polynomial with error if conversion is not possible
## This needs a hack, as Polynomials.jl defines `/` to return a `div` and not an error for p(x)/q(x)
immutable PolyTest
x
end
import Base: +,-,*,/,^
+(a::PolyTest, b::PolyTest) = PolyTest(a.x + b.x)
+{T<:Number}(a::T, b::PolyTest) = PolyTest(a + b.x)
+{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x + b)
-(a::PolyTest,b::PolyTest) = PolyTest(a.x - b.x)
-{T<:Number}(a::T, b::PolyTest) = PolyTest(a - b.x)
-{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x - b)
-(a::PolyTest) = PolyTest(-a.x)
*(a::PolyTest, b::PolyTest) = PolyTest(a.x * b.x)
*(a::Bool, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(a::T, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(b::PolyTest, a::T) = PolyTest(a * b.x)
/{T<:Number}(b::PolyTest, a::T) = PolyTest(b.x / a)
^{T<:Integer}(b::PolyTest, a::T) = PolyTest(b.x ^ a)
immutable PolyTest
x
end
import Base: +,-,*,/,^
+(a::PolyTest, b::PolyTest) = PolyTest(a.x + b.x)
+{T<:Number}(a::T, b::PolyTest) = PolyTest(a + b.x)
+{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x + b)
-(a::PolyTest,b::PolyTest) = PolyTest(a.x - b.x)
-{T<:Number}(a::T, b::PolyTest) = PolyTest(a - b.x)
-{T<:Number}(a::PolyTest,b::T) = PolyTest(a.x - b)
-(a::PolyTest) = PolyTest(-a.x)
*(a::PolyTest, b::PolyTest) = PolyTest(a.x * b.x)
*(a::Bool, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(a::T, b::PolyTest) = PolyTest(a * b.x)
*{T<:Number}(b::PolyTest, a::T) = PolyTest(a * b.x)
/{T<:Number}(b::PolyTest, a::T) = PolyTest(b.x / a)
^{T<:Integer}(b::PolyTest, a::T) = PolyTest(b.x ^ a)

typealias QQR @compat Union{Int, BigInt, Rational{Int}, Rational{BigInt}, Float64}
const QQR = Union{Int, BigInt, Rational{Int}, Rational{BigInt}, Float64}
function Base.convert{T<:QQR}(::Type{Poly{T}}, f::Function)
try
f(PolyTest(0)) # error if not a poly
Expand Down
28 changes: 14 additions & 14 deletions src/Polys/real_roots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
## variable(p::Poly) = poly(zeros(eltype(p),1), p.var)

## make bounds work for floating point or rational.
_iszero{T <: (@compat Union{Integer, Rational})}(b::T; kwargs...) = b == 0
_iszero{T <: Union{Integer, Rational}}(b::T; kwargs...) = b == 0
_iszero{T<:AbstractFloat}(b::T; xtol=1) = abs(b) <= 2*xtol*eps(T)


Expand Down Expand Up @@ -42,13 +42,13 @@ end


## Our Poly types for which we can find gcd
typealias QQ @compat Union{Int, BigInt, Rational{Int}, Rational{BigInt}}
typealias BB @compat Union{BigInt, Rational{BigInt}}
const QQ = Union{Int, BigInt, Rational{Int}, Rational{BigInt}}
const BB = Union{BigInt, Rational{BigInt}}

## Here computations are exact, as long as we return poly in Q[x]
function Base.divrem{T<:QQ, S<:QQ}(a::Poly{T}, b::Poly{S})
degree(b) == 0 && (b[0] == 0 ? error("b is 0") : return (a/b[0], 0*b))

lc(p::Poly) = p[degree(p)]

a.var == b.var || DomainError() # "symbols must match"
Expand Down Expand Up @@ -94,7 +94,7 @@ end
## titan.princeton.edu/papers/claire/hertz-etal-99.ps
function upperbound(p::Poly)
q, d = monic(p), degree(p)

d == 0 && error("degree 0 is a constant")
d == 1 && abs(q[0])

Expand All @@ -116,7 +116,7 @@ Use Descarte's rule of signs to compute number of *postive* roots of p
¬(f::Function) = x -> !f(x)
function descartes_bound(p::Poly)
bs = filter(¬_iszero, p.a)

## count terms
ctr = 0
@compat bs = sign.(bs)
Expand Down Expand Up @@ -145,7 +145,7 @@ immutable Intervalkc
end

## convenience for two useful functions
DesBound(p,node::Intervalkc) = DesBound(Pkc(p,node))
DesBound(p,node::Intervalkc) = DesBound(Pkc(p,node))
function Pkc(p::Poly, node::Intervalkc)
k,c = node.k, node.c
Pkc(p, k, c)
Expand Down Expand Up @@ -205,7 +205,7 @@ end
## strategy, the push!(st.Internal, i) would be changed in `addSucc`.
function find_in_01(p::Poly)
st = State(p)

initTree(st)
for x in [0,1]
p,k = multiplicity(p,x)
Expand Down Expand Up @@ -255,7 +255,7 @@ function pos_roots(p::Poly)
c = -p[0]/p[1]
if c > 0 return([c]) else return(rts) end
end

## in case 0 is a root, we remove from p. Should be unnecessary, but ...
p,k = multiplicity(p, 0)

Expand Down Expand Up @@ -292,9 +292,9 @@ function real_roots_sqfree(p::Poly)
## Handle simple cases
degree(p) == 0 && error("Roots for degree 0 polynomials are either empty or all of R.")
degree(p) == 1 && return([ -p[0]/p[1] ])

rts = Any[]

p,k = multiplicity(p, 0); k > 0 && push!(rts, 0) # 0
append!(rts, pos_roots(p)) # positive roots
append!(rts, map(-,neg_roots(p))) # negative roots
Expand All @@ -313,11 +313,11 @@ function real_roots{T <: Union{Int, BigInt}}(p::Poly{T})
rts = Real[]
rat_rts = rational_roots(p)
append!(rts, rat_rts)

for rt in rat_rts
f, k = PolynomialFactors.synthetic_division(f, rt)
end

if degree(f) > 0
real_rts = real_roots_sqfree(f)
append!(rts, real_rts)
Expand All @@ -333,7 +333,7 @@ function real_roots(p::Poly)
if degree(p) > 1
u, p, w, residual= agcd(p)
end

real_roots_sqfree(p)
end

Expand Down
38 changes: 19 additions & 19 deletions src/bracketing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,18 +31,18 @@ function _middle(x::Float64, y::Float64)
if !isfinite(x) || !isfinite(y)
return x + y
end

# Always return 0.0 when inputs have opposite sign
if sign(x) != sign(y) && x != 0.0 && y != 0.0
return 0.0
end

negate = x < 0.0 || y < 0.0

xint = reinterpret(UInt64, abs(x))
yint = reinterpret(UInt64, abs(y))
unsigned = reinterpret(Float64, (xint + yint) >> 1)

return negate ? -unsigned : unsigned
end

Expand All @@ -54,11 +54,11 @@ end

####
## find_zero interface. We need to specialize for T<:Float64, and BigSomething
typealias BigSomething @compat Union{BigFloat, BigInt}
const BigSomething = Union{BigFloat, BigInt}

abstract AbstractBisection <: UnivariateZeroMethod
@compat abstract type AbstractBisection <: UnivariateZeroMethod end
type Bisection <: AbstractBisection end
type A42 <: AbstractBisection end
type A42 <: AbstractBisection end

## This is a bit clunky, as we use `a42` for bisection when we don't have `Float64` values.
## As we don't have the `A42` algorithm implemented through `find_zero`, we adjust here.
Expand Down Expand Up @@ -98,7 +98,7 @@ function init_state{T <: Float64}(method::Bisection, fs, x::Vector{T}, bracket)
@compat y0, y2 = fs.f.([x0, x2])

sign(y0) * sign(y2) > 0 && (warn(bracketing_error); throw(ArgumentError))

state = UnivariateZeroState(x2, x0,
y2, y0,
isa(bracket, Nullable) ? bracket : Nullable(convert(Vector{T}, sort(bracket))),
Expand Down Expand Up @@ -139,9 +139,9 @@ function assess_convergence(method::Bisection, fs, state, options)
if x0 > x2
x0, x2 = x2, x0
end

x1 = _middle(x0, x2)

x0 < x1 && x1 < x2 && return false

state.message = ""
Expand All @@ -152,7 +152,7 @@ end

##################################################

"""
"""

Finds the root of a continuous function within a provided
interval [a, b], without requiring derivatives. It is based on algorithm 4.2
Expand All @@ -174,14 +174,14 @@ By John Travers

"""
function a42(f, a, b;
xtol=zero(float(a)),
xtol=zero(float(a)),
maxeval::Int=15,
verbose::Bool=false
)
if a > b
a,b = b,a
end

if a >= b || sign(f(a))*sign(f(b)) >= 0
error("on input a < b and f(a)f(b) < 0 must both hold")
end
Expand All @@ -203,12 +203,12 @@ Solve f(x) = 0 over bracketing interval [a,b] starting at c, with a < c < b

"""
function a42a(f, a, b, c=(0.5)*(a+b);
xtol=zero(float(a)),
xtol=zero(float(a)),
maxeval::Int=15,
verbose::Bool=false
)


try
# re-bracket and check termination
a, b, d = bracket(f, a, b, c, xtol)
Expand Down Expand Up @@ -250,7 +250,7 @@ function a42a(f, a, b, c=(0.5)*(a+b);
end

verbose && println(@sprintf("a=%18.15f, n=%s", float(a), n))

if nextfloat(ch) * prevfloat(ch) <= 0
throw(StateConverged(ch))
end
Expand Down Expand Up @@ -329,7 +329,7 @@ function bracket(f, a, b, c, tol)
faa = f(aa)
fbb = f(bb)
if bb - aa < 2*tole(aa, bb, faa, fbb, tol)
x0 = abs(faa) < abs(fbb) ? aa : bb
x0 = abs(faa) < abs(fbb) ? aa : bb
throw(StateConverged(x0))
end
aa, bb, db
Expand Down Expand Up @@ -372,7 +372,7 @@ end


# approximate zero of f using inverse cubic interpolation
# if the new guess is outside [a, b] we use a quadratic step instead
# if the new guess is outside [a, b] we use a quadratic step instead
# based on algorithm on page 333 of [1]
function ipzero(f, a, b, c, d)
fa = f(a)
Expand Down Expand Up @@ -419,7 +419,7 @@ Searches for zeros of `f` in an interval [a, b].

Basic algorithm used:

* split interval [a,b] into `no_pts` subintervals.
* split interval [a,b] into `no_pts` subintervals.
* For each bracketing interval finds a bracketed zero.
* For other subintervals does a quick search with a derivative-free method.

Expand Down
Loading