diff -Nru julia-0.3.0/base/array.jl julia-0.3.0/base/array.jl --- julia-0.3.0/base/array.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/array.jl 2014-06-03 19:33:18.000000000 +0000 @@ -1145,7 +1145,7 @@ # similar to Matlab's ismember # returns a vector containing the highest index in b for each value in a that is a member of b -function indexin{T}(a::AbstractArray{T}, b::AbstractArray{T}) +function indexin(a::AbstractArray, b::AbstractArray) bdict = Dict(b, 1:length(b)) [get(bdict, i, 0) for i in a] end diff -Nru julia-0.3.0/base/combinatorics.jl julia-0.3.0/base/combinatorics.jl --- julia-0.3.0/base/combinatorics.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/combinatorics.jl 2014-06-03 19:33:18.000000000 +0000 @@ -41,6 +41,14 @@ factorial(n::Union(Int8,Uint8,Int16,Uint16,Int32,Uint32)) = factorial(int64(n)) end +function gamma(n::Union(Int8,Uint8,Int16,Uint16,Int32,Uint32,Int64,Uint64)) + n < 0 && throw(DomainError()) + n == 0 && return Inf + n <= 2 && return 1.0 + n > 20 && return gamma(float64(n)) + @inbounds return float64(_fact_table64[n-1]) +end + function factorial(n::Integer) n < 0 && throw(DomainError()) local f::typeof(n*n), i::typeof(n*n) @@ -378,10 +386,14 @@ length(f::FixedPartitions) = npartitions(f.n,f.m) -partitions(n::Integer, m::Integer) = (@assert 2 <= m <= n; FixedPartitions(n,m)) +partitions(n::Integer, m::Integer) = n >= 1 && m >= 1 ? FixedPartitions(n,m) : throw(DomainError()) start(f::FixedPartitions) = Int[] -done(f::FixedPartitions, s::Vector{Int}) = !isempty(s) && s[1]-1 <= s[end] +function done(f::FixedPartitions, s::Vector{Int}) + f.m <= f.n || return true + isempty(s) && return false + return f.m == 1 || s[1]-1 <= s[end] +end next(f::FixedPartitions, s::Vector{Int}) = (xs = nextfixedpartition(f.n,f.m,s); (xs,xs)) function nextfixedpartition(n, m, bs) @@ -439,7 +451,7 @@ partitions(s::AbstractVector) = SetPartitions(s) start(p::SetPartitions) = (n = length(p.s); (zeros(Int32, n), ones(Int32, n-1), n, 1)) -done(p::SetPartitions, s) = !isempty(s) && s[1][1] > 0 +done(p::SetPartitions, s) = s[1][1] > 0 next(p::SetPartitions, s) = nextsetpartition(p.s, s...) function nextsetpartition(s::AbstractVector, a, b, n, m) @@ -503,13 +515,19 @@ length(p::FixedSetPartitions) = nfixedsetpartitions(length(p.s),p.m) -partitions(s::AbstractVector,m::Int) = (@assert 2 <= m <= length(s); FixedSetPartitions(s,m)) +partitions(s::AbstractVector,m::Int) = length(s) >= 1 && m >= 1 ? FixedSetPartitions(s,m) : throw(DomainError()) -start(p::FixedSetPartitions) = (n = length(p.s);m=p.m; (vcat(ones(Int, n-m),1:m), vcat(1,n-m+2:n), n)) -# state consists of vector a of length n describing to which partition every element of s belongs and -# a vector b of length m describing the first index b[i] that belongs to partition i +function start(p::FixedSetPartitions) + n = length(p.s) + m = p.m + m <= n ? (vcat(ones(Int, n-m),1:m), vcat(1,n-m+2:n), n) : (Int[], Int[], n) +end +# state consists of: +# vector a of length n describing to which partition every element of s belongs +# vector b of length n describing the first index b[i] that belongs to partition i +# integer n -done(p::FixedSetPartitions, s) = !isempty(s) && s[1][1] > 1 +done(p::FixedSetPartitions, s) = length(s[1]) == 0 || s[1][1] > 1 next(p::FixedSetPartitions, s) = nextfixedsetpartition(p.s,p.m, s...) function nextfixedsetpartition(s::AbstractVector, m, a, b, n) @@ -523,6 +541,11 @@ part = makeparts(s,a) + if m == 1 + a[1] = 2 + return (part, (a, b, n)) + end + if a[end] != m a[end] += 1 else @@ -565,7 +588,7 @@ end -# For a list of integers i1, i2, i3, find the smallest +# For a list of integers i1, i2, i3, find the smallest # i1^n1 * i2^n2 * i3^n3 >= x # for integer n1, n2, n3 function nextprod(a::Vector{Int}, x) @@ -579,7 +602,7 @@ p::widen(Int) = mx[1] # initial value of product in this case best = p icarry = 1 - + while v[end] < mx[end] if p >= x best = p < best ? p : best # keep the best found yet @@ -621,7 +644,7 @@ p::widen(Int) = first best = p icarry = 1 - + while v[end] < mx[end] while p <= x best = p > best ? p : best diff -Nru julia-0.3.0/base/datafmt.jl julia-0.3.0/base/datafmt.jl --- julia-0.3.0/base/datafmt.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/datafmt.jl 2014-06-03 19:33:18.000000000 +0000 @@ -245,13 +245,18 @@ comments = get(optsd, :comments, true) comment_char = get(optsd, :comment_char, '#') dims = get(optsd, :dims, nothing) - has_header = get(optsd, :has_header, false) + + has_header = get(optsd, :header, get(optsd, :has_header, false)) + haskey(optsd, :has_header) && (optsd[:has_header] != has_header) && error("conflicting values for header and has_header") + + skipstart = get(optsd, :skipstart, 0) + (skipstart >= 0) || error("invalid value for skipstart") offset_handler = (dims == nothing) ? DLMOffsets(sbuff) : DLMStore(T, dims, has_header, sbuff, auto, eol) for retry in 1:2 try - dims = dlm_parse(sbuff, eol, dlm, '"', comment_char, ign_empty, quotes, comments, offset_handler) + dims = dlm_parse(sbuff, eol, dlm, '"', comment_char, ign_empty, quotes, comments, skipstart, offset_handler) break catch ex if isa(ex, TypeError) && (ex.func == :store_cell) @@ -274,15 +279,17 @@ return readdlm_string(sbuff, dlm, T, eol, auto, optsd) end -const valid_opts = [:has_header, :ignore_invalid_chars, :use_mmap, :quotes, :comments, :dims, :comment_char] -const valid_opt_types = [Bool, Bool, Bool, Bool, Bool, NTuple{2,Integer}, Char] +const valid_opts = [:header, :has_header, :ignore_invalid_chars, :use_mmap, :quotes, :comments, :dims, :comment_char, :skipstart] +const valid_opt_types = [Bool, Bool, Bool, Bool, Bool, Bool, NTuple{2,Integer}, Char, Integer] +const deprecated_opts = [ :has_header => :header ] function val_opts(opts) - d = Dict{Symbol,Union(Bool,NTuple{2,Integer},Char)}() - for opt in opts - !in(opt[1], valid_opts) && error("unknown option $(opt[1])") - opt_typ = valid_opt_types[findfirst(valid_opts, opt[1])] - !isa(opt[2], opt_typ) && error("$(opt[1]) should be of type $opt_typ") - d[opt[1]] = opt[2] + d = Dict{Symbol,Union(Bool,NTuple{2,Integer},Char,Integer)}() + for (opt_name, opt_val) in opts + !in(opt_name, valid_opts) && error("unknown option $opt_name") + opt_typ = valid_opt_types[findfirst(valid_opts, opt_name)] + !isa(opt_val, opt_typ) && error("$opt_name should be of type $opt_typ") + d[opt_name] = opt_val + haskey(deprecated_opts, opt_name) && warn("$opt_name is deprecated, use $(deprecated_opts[opt_name]) instead") end d end @@ -319,14 +326,17 @@ colval{S<:String}(sval::S, cells::Array, row::Int, col::Int, tmp64::Array{Float64,1}) = true -dlm_parse(s::ASCIIString, eol::Char, dlm::Char, qchar::Char, cchar::Char, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, dh::DLMHandler) = dlm_parse(s.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, dh) -function dlm_parse{T,D}(dbuff::T, eol::D, dlm::D, qchar::D, cchar::D, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, dh::DLMHandler) +dlm_parse(s::ASCIIString, eol::Char, dlm::Char, qchar::Char, cchar::Char, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, skipstart::Int, dh::DLMHandler) = + dlm_parse(s.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, skipstart, dh) + +function dlm_parse{T,D}(dbuff::T, eol::D, dlm::D, qchar::D, cchar::D, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, skipstart::Int, dh::DLMHandler) all_ascii = (D <: Uint8) || (isascii(eol) && isascii(dlm) && (!allow_quote || isascii(qchar)) && (!allow_comments || isascii(cchar))) - (T <: UTF8String) && all_ascii && (return dlm_parse(dbuff.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, dh)) + (T <: UTF8String) && all_ascii && (return dlm_parse(dbuff.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, skipstart, dh)) ncols = nrows = col = 0 is_default_dlm = (dlm == convert(D, invalid_dlm)) error_str = "" - state = 0 # 0: begin field, 1: quoted field, 2: unquoted field, 3: second quote (could either be end of field or escape character), 4: comment + # 0: begin field, 1: quoted field, 2: unquoted field, 3: second quote (could either be end of field or escape character), 4: comment, 5: skipstart + state = (skipstart > 0) ? 5 : 0 is_eol = is_dlm = is_cr = is_quote = is_comment = expct_col = false idx = 1 try @@ -442,6 +452,12 @@ error_str = escape_string("unexpected character '$(char(val))' after quoted field at row $(nrows+1) column $(col+1)") break end + elseif 5 == state # skip start + if is_eol + col_start_idx = idx + skipstart -= 1 + (0 == skipstart) && (state = 0) + end end was_cr = is_cr end diff -Nru julia-0.3.0/base/dict.jl julia-0.3.0/base/dict.jl --- julia-0.3.0/base/dict.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/dict.jl 2014-06-03 19:33:18.000000000 +0000 @@ -11,28 +11,96 @@ !is(v, secret_table_token) && (v == p[2]) end -function show{K,V}(io::IO, t::Associative{K,V}) - if isempty(t) - print(io, typeof(t),"()") - else - if K === Any && V === Any - delims = ['{','}'] - else - delims = ['[',']'] - end +function summary(t::Associative) + n = length(t) + string(typeof(t), " with ", n, (n==1 ? " entry" : " entries")) +end + +function showcompact{K,V}(io::IO, t::Associative{K,V}) + print(io, summary(t)) + if !isempty(t) + print(io, ": ") + delims = (K == V == Any) ? ('{', '}') : ('[', ']') print(io, delims[1]) first = true - for (k, v) = t - first || print(io, ',') - first = false - show(io, k) + for (k, v) in t + first || (print(io, ','); first = false) + showcompact(io, k) print(io, "=>") - show(io, v) + showcompact(io, v) end print(io, delims[2]) end end +function _truncate_at_width_or_chars(str, width, chars="", truncmark="…") + truncwidth = strwidth(truncmark) + (width <= 0 || width < truncwidth) && return "" + + wid = truncidx = lastidx = 0 + idx = start(str) + while !done(str, idx) + lastidx = idx + c, idx = next(str, idx) + wid += charwidth(c) + wid >= width - truncwidth && truncidx == 0 && (truncidx = lastidx) + (wid >= width || c in chars) && break + end + + str[lastidx] in chars && (lastidx = prevind(str, lastidx)) + truncidx == 0 && (truncidx = lastidx) + if lastidx < sizeof(str) + return bytestring(SubString(str, 1, truncidx) * truncmark) + else + return bytestring(str) + end +end + +function showdict{K,V}(io::IO, t::Associative{K,V}, limit_output::Bool = false, + rows = tty_rows()-3, cols = tty_cols()) + print(io, summary(t)) + isempty(t) && return + print(io, ":") + + if limit_output + rows < 2 && (print(io, " …"); return) + cols < 12 && (cols = 12) # Minimum widths of 2 for key, 4 for value + cols -= 6 # Subtract the widths of prefix " " separator " => " + rows -= 2 # Subtract the summary and final ⋮ continuation lines + + # determine max key width to align the output, caching the strings + ks = Array(String, min(rows, length(t))) + keylen = 0 + for (i, k) in enumerate(keys(t)) + i > rows && break + ks[i] = sprint(show, k) + keylen = clamp(length(ks[i]), keylen, div(cols, 3)) + end + end + + for (i, (k, v)) in enumerate(t) + print(io, "\n ") + limit_output && i > rows && (print(io, rpad("⋮", keylen), " => ⋮"); break) + + if limit_output + key = rpad(_truncate_at_width_or_chars(ks[i], keylen, "\r\n"), keylen) + else + key = sprint(show, k) + end + print(io, key) + print(io, " => ") + + val = sprint(show, v) + if limit_output + val = _truncate_at_width_or_chars(val, cols - keylen, "\r\n") + end + print(io, val) + end +end + +show{K,V}(io::IO, t::Associative{K,V}) = showdict(io, t, false) +showlimited{K,V}(io::IO, t::Associative{K,V}) = showdict(io, t, true) + immutable KeyIterator{T<:Associative} dict::T end @@ -40,6 +108,34 @@ dict::T end +summary{T<:Union(KeyIterator,ValueIterator)}(iter::T) = + string(T.name, " for a ", summary(iter.dict)) + +show(io::IO, iter::Union(KeyIterator,ValueIterator)) = showkv(io, iter, false) +showlimited(io::IO, iter::Union(KeyIterator,ValueIterator)) = showkv(io, iter, true) + +function showkv{T<:Union(KeyIterator,ValueIterator)}(io::IO, iter::T, limit_output::Bool = false, + rows = tty_rows()-3, cols = tty_cols()) + print(io, summary(iter)) + isempty(iter) && return + print(io, ". ", T<:KeyIterator ? "Keys" : "Values", ":") + if limit_output + rows < 2 && (print(io, " …"); return) + cols < 4 && (cols = 4) + cols -= 2 # For prefix " " + rows -= 2 # For summary and final ⋮ continuation lines + end + + for (i, v) in enumerate(iter) + print(io, "\n ") + limit_output && i >= rows && (print(io, "⋮"); break) + + str = sprint(show, v) + limit_output && (str = _truncate_at_width_or_chars(str, cols, "\r\n")) + print(io, str) + end +end + length(v::Union(KeyIterator,ValueIterator)) = length(v.dict) isempty(v::Union(KeyIterator,ValueIterator)) = isempty(v.dict) eltype(v::KeyIterator) = eltype(v.dict)[1] diff -Nru julia-0.3.0/base/float.jl julia-0.3.0/base/float.jl --- julia-0.3.0/base/float.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/float.jl 2014-06-03 19:33:18.000000000 +0000 @@ -72,7 +72,7 @@ if WORD_SIZE == 64 iround(x::Float32) = iround(float64(x)) - itrunc(x::Float32) = itrunc(float64(x)) + itrunc(x::Float32) = box(Int64,fptosi(Int64,unbox(Float32,x))) iround(x::Float64) = box(Int64,fpsiround(unbox(Float64,x))) itrunc(x::Float64) = box(Int64,fptosi(unbox(Float64,x))) else @@ -175,10 +175,10 @@ ==(x::Int64 , y::Float64) = eqfsi64(unbox(Float64,y),unbox(Int64,x)) ==(x::Uint64 , y::Float64) = eqfui64(unbox(Float64,y),unbox(Uint64,x)) -==(x::Float32, y::Int64 ) = eqfsi64(unbox(Float64,float64(x)),unbox(Int64,y)) -==(x::Float32, y::Uint64 ) = eqfui64(unbox(Float64,float64(x)),unbox(Uint64,y)) -==(x::Int64 , y::Float32) = eqfsi64(unbox(Float64,float64(y)),unbox(Int64,x)) -==(x::Uint64 , y::Float32) = eqfui64(unbox(Float64,float64(y)),unbox(Uint64,x)) +==(x::Float32, y::Int64 ) = eqfsi64(unbox(Float32,x),unbox(Int64,y)) +==(x::Float32, y::Uint64 ) = eqfui64(unbox(Float32,x),unbox(Uint64,y)) +==(x::Int64 , y::Float32) = eqfsi64(unbox(Float32,y),unbox(Int64,x)) +==(x::Uint64 , y::Float32) = eqfui64(unbox(Float32,y),unbox(Uint64,x)) < (x::Float64, y::Int64 ) = ltfsi64(unbox(Float64,x),unbox(Int64,y)) < (x::Float64, y::Uint64 ) = ltfui64(unbox(Float64,x),unbox(Uint64,y)) diff -Nru julia-0.3.0/base/inference.jl julia-0.3.0/base/inference.jl --- julia-0.3.0/base/inference.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/inference.jl 2014-06-03 19:33:18.000000000 +0000 @@ -934,6 +934,10 @@ t0 = abstract_eval(e.args[1], vtypes, sv) # intersect with Any to remove Undef t = typeintersect(t0, Any) + if isa(t,DataType) && typeseq(t,t.name.primary) + # remove unnecessary typevars + t = t.name.primary + end if is(t,None) && Undef<:t0 # the first time we see this statement the variable will probably # be Undef; return None so this doesn't contribute to the type diff -Nru julia-0.3.0/base/linalg/dense.jl julia-0.3.0/base/linalg/dense.jl --- julia-0.3.0/base/linalg/dense.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/linalg/dense.jl 2014-06-03 19:33:18.000000000 +0000 @@ -30,13 +30,12 @@ stride1(x::Array) = 1 stride1(x::StridedVector) = stride(x, 1)::Int -Base.sum_seq{T<:BlasFloat}(::Base.AbsFun, a::Union(Array{T},StridedVector{T}), ifirst::Int, ilast::Int) = - BLAS.asum(ilast-ifirst+1, pointer(a, ifirst), stride1(a)) +import Base: mapreduce_seq_impl, AbsFun, Abs2Fun, AddFun -# This appears to show a benefit from a larger block size -Base.sum_pairwise_blocksize(::Base.Abs2Fun) = 4096 +mapreduce_seq_impl{T<:BlasFloat}(::AbsFun, ::AddFun, a::Union(Array{T},StridedVector{T}), ifirst::Int, ilast::Int) = + BLAS.asum(ilast-ifirst+1, pointer(a, ifirst), stride1(a)) -function Base.sum_seq{T<:BlasFloat}(::Base.Abs2Fun, a::Union(Array{T},StridedVector{T}), ifirst::Int, ilast::Int) +function mapreduce_seq_impl{T<:BlasFloat}(::Abs2Fun, ::AddFun, a::Union(Array{T},StridedVector{T}), ifirst::Int, ilast::Int) n = ilast-ifirst+1 px = pointer(a, ifirst) incx = stride1(a) diff -Nru julia-0.3.0/base/math.jl julia-0.3.0/base/math.jl --- julia-0.3.0/base/math.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/math.jl 2014-06-03 19:33:18.000000000 +0000 @@ -48,102 +48,6 @@ ex end -function sinpi(x::Real) - if isinf(x) - return throw(DomainError()) - elseif isnan(x) - return nan(x) - end - - rx = copysign(float(rem(x,2)),x) - arx = abs(rx) - - if arx < oftype(rx,0.25) - return sin(pi*rx) - elseif arx <= oftype(rx,0.75) - arx = oftype(rx,0.5) - arx - return copysign(cos(pi*arx),rx) - elseif arx < oftype(rx,1.25) - rx = (one(rx) - arx)*sign(rx) - return sin(pi*rx) - elseif arx <= oftype(rx,1.75) - arx = oftype(rx,1.5) - arx - return -copysign(cos(pi*arx),rx) - else - rx = rx - copysign(oftype(rx,2.0),rx) - return sin(pi*rx) - end -end - -function cospi(x::Real) - if isinf(x) - return throw(DomainError()) - elseif isnan(x) - return nan(x) - end - - rx = abs(float(rem(x,2))) - - if rx <= oftype(rx,0.25) - return cos(pi*rx) - elseif rx < oftype(rx,0.75) - rx = oftype(rx,0.5) - rx - return sin(pi*rx) - elseif rx <= oftype(rx,1.25) - rx = one(rx) - rx - return -cos(pi*rx) - elseif rx < oftype(rx,1.75) - rx = rx - oftype(rx,1.5) - return sin(pi*rx) - else - rx = oftype(rx,2.0) - rx - return cos(pi*rx) - end -end - -sinpi(x::Integer) = zero(x) -cospi(x::Integer) = isodd(x) ? -one(x) : one(x) - -function sinpi(z::Complex) - zr, zi = reim(z) - if !isfinite(zi) && zr == 0 return complex(zr, zi) end - if isnan(zr) && !isfinite(zi) return complex(zr, zi) end - if !isfinite(zr) && zi == 0 return complex(oftype(zr, NaN), zi) end - if !isfinite(zr) && isfinite(zi) return complex(oftype(zr, NaN), oftype(zi, NaN)) end - if !isfinite(zr) && !isfinite(zi) return complex(zr, oftype(zi, NaN)) end - pizi = pi*zi - complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi)) -end - -function cospi(z::Complex) - zr, zi = reim(z) - if !isfinite(zi) && zr == 0 - return complex(isnan(zi) ? zi : oftype(zi, Inf), - isnan(zi) ? zr : zr*-sign(zi)) - end - if !isfinite(zr) && isinf(zi) - return complex(oftype(zr, Inf), oftype(zi, NaN)) - end - if isinf(zr) - return complex(oftype(zr, NaN), zi==0 ? -copysign(zi, zr) : oftype(zi, NaN)) - end - if isnan(zr) && zi==0 return complex(zr, abs(zi)) end - pizi = pi*zi - complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi)) -end -@vectorize_1arg Number sinpi -@vectorize_1arg Number cospi - - -sinc(x::Number) = x==0 ? one(x) : oftype(x,sinpi(x)/(pi*x)) -sinc(x::Integer) = x==0 ? one(x) : zero(x) -sinc{T<:Integer}(x::Complex{T}) = sinc(float(x)) -@vectorize_1arg Number sinc -cosc(x::Number) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x) -cosc(x::Integer) = cosc(float(x)) -cosc{T<:Integer}(x::Complex{T}) = cosc(float(x)) -@vectorize_1arg Number cosc - rad2deg(z::Real) = oftype(z, 57.29577951308232*z) deg2rad(z::Real) = oftype(z, 0.017453292519943295*z) rad2deg(z::Integer) = rad2deg(float(z)) @@ -151,95 +55,6 @@ @vectorize_1arg Real rad2deg @vectorize_1arg Real deg2rad -for (finv, f) in ((:sec, :cos), (:csc, :sin), (:cot, :tan), - (:sech, :cosh), (:csch, :sinh), (:coth, :tanh), - (:secd, :cosd), (:cscd, :sind), (:cotd, :tand)) - @eval begin - ($finv){T<:Number}(z::T) = one(T) / (($f)(z)) - ($finv){T<:Number}(z::AbstractArray{T}) = one(T) ./ (($f)(z)) - end -end - -for (fa, fainv) in ((:asec, :acos), (:acsc, :asin), (:acot, :atan), - (:asech, :acosh), (:acsch, :asinh), (:acoth, :atanh)) - @eval begin - ($fa){T<:Number}(y::T) = ($fainv)(one(T) / y) - ($fa){T<:Number}(y::AbstractArray{T}) = ($fainv)(one(T) ./ y) - end -end - -function sind(x::Real) - if isinf(x) - return throw(DomainError()) - elseif isnan(x) - return nan(x) - end - - rx = copysign(float(rem(x,360)),x) - arx = abs(rx) - - if arx < oftype(rx,45.0) - return sin(deg2rad(rx)) - elseif arx <= oftype(rx,135.0) - arx = oftype(rx,90.0) - arx - return copysign(cos(deg2rad(arx)),rx) - elseif arx < oftype(rx,225.0) - rx = (oftype(rx,180.0) - arx)*sign(rx) - return sin(deg2rad(rx)) - elseif arx <= 315.0 - arx = oftype(rx,270.0) - arx - return -copysign(cos(deg2rad(arx)),rx) - else - rx = rx - copysign(oftype(rx,360.0),rx) - return sin(deg2rad(rx)) - end -end -@vectorize_1arg Real sind - -function cosd(x::Real) - if isinf(x) - return throw(DomainError()) - elseif isnan(x) - return nan(x) - end - - rx = abs(float(rem(x,360))) - - if rx <= oftype(rx,45.0) - return cos(deg2rad(rx)) - elseif rx < oftype(rx,135.0) - rx = oftype(rx,90.0) - rx - return sin(deg2rad(rx)) - elseif rx <= oftype(rx,225.0) - rx = oftype(rx,180.0) - rx - return -cos(deg2rad(rx)) - elseif rx < oftype(rx,315.0) - rx = rx - oftype(rx,270.0) - return sin(deg2rad(rx)) - else - rx = oftype(rx,360.0) - rx - return cos(deg2rad(rx)) - end -end -@vectorize_1arg Real cosd - -tand(x::Real) = sind(x) / cosd(x) -@vectorize_1arg Real tand - -for (fd, f) in ((:sind, :sin), (:cosd, :cos), (:tand, :tan)) - @eval begin - ($fd)(z) = ($f)(deg2rad(z)) - end -end - -for (fd, f) in ((:asind, :asin), (:acosd, :acos), (:atand, :atan), - (:asecd, :asec), (:acscd, :acsc), (:acotd, :acot)) - @eval begin - ($fd)(y) = rad2deg(($f)(y)) - @vectorize_1arg Real $fd - end -end - log(b,x) = log(x)./log(b) # type specific math functions @@ -331,26 +146,6 @@ end end -gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm), Float64, (Float64,), x), x) -gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm), Float32, (Float32,), x), x) -gamma(x::Real) = gamma(float(x)) -@vectorize_1arg Number gamma - -function lgamma_r(x::Float64) - signp = Array(Int32, 1) - y = ccall((:lgamma_r,libm), Float64, (Float64, Ptr{Int32}), x, signp) - return y, signp[1] -end -function lgamma_r(x::Float32) - signp = Array(Int32, 1) - y = ccall((:lgammaf_r,libm), Float32, (Float32, Ptr{Int32}), x, signp) - return y, signp[1] -end -lgamma_r(x::Real) = lgamma_r(float(x)) - -lfact(x::Real) = (x<=1 ? zero(float(x)) : lgamma(x+one(x))) -@vectorize_1arg Number lfact - max{T<:FloatingPoint}(x::T, y::T) = ifelse((y > x) | (x != x), y, x) @vectorize_2arg Real max @@ -412,6 +207,18 @@ modf(x) = rem(x,one(x)), trunc(x) +const _modff_temp = Float32[0] +function modf(x::Float32) + f = ccall((:modff,libm), Float32, (Float32,Ptr{Float32}), x, _modff_temp) + f, _modff_temp[1] +end + +const _modf_temp = Float64[0] +function modf(x::Float64) + f = ccall((:modf,libm), Float64, (Float64,Ptr{Float64}), x, _modf_temp) + f, _modf_temp[1] +end + ^(x::Float64, y::Float64) = nan_dom_err(ccall((:pow,libm), Float64, (Float64,Float64), x, y), x+y) ^(x::Float32, y::Float32) = nan_dom_err(ccall((:powf,libm), Float32, (Float32,Float32), x, y), x+y) @@ -420,290 +227,6 @@ ^(x::Float32, y::Integer) = box(Float32, powi_llvm(unbox(Float32,x), unbox(Int32,int32(y)))) -# special functions - -for jy in ("j","y"), nu in (0,1) - jynu = Expr(:quote, symbol(string(jy,nu))) - jynuf = Expr(:quote, symbol(string(jy,nu,"f"))) - bjynu = symbol(string("bessel",jy,nu)) - if jy == "y" - @eval begin - $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm), Float64, (Float64,), x), x) - $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x) - end - else - @eval begin - $bjynu(x::Float64) = ccall(($jynu,libm), Float64, (Float64,), x) - $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x) - end - end - @eval begin - $bjynu(x::Real) = $bjynu(float(x)) - $bjynu(x::Complex) = $(symbol(string("bessel",jy)))($nu,x) - @vectorize_1arg Number $bjynu - end -end - - -type AmosException <: Exception - info::Int32 -end - -let - const ai::Array{Float64,1} = Array(Float64,2) - const ae::Array{Int32,1} = Array(Int32,2) -global airy -function airy(k::Int, z::Complex128) - id = int32(k==1 || k==3) - if k == 0 || k == 1 - ccall((:zairy_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), - &id, &1, - pointer(ai,1), pointer(ai,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(ai[1],ai[2]) - else - throw(AmosException(ae[2])) - end - elseif k == 2 || k == 3 - ccall((:zbiry_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), - &real(z), &imag(z), - &id, &1, - pointer(ai,1), pointer(ai,2), - pointer(ae,1)) - if ae[1] == 0 || ae[1] == 3 # ignore underflow - return complex(ai[1],ai[2]) - else - throw(AmosException(ae[2])) - end - else - error("invalid argument") - end -end -end - -airy(z) = airy(0,z) -@vectorize_1arg Number airy -airyprime(z) = airy(1,z) -@vectorize_1arg Number airyprime -airyai(z) = airy(0,z) -@vectorize_1arg Number airyai -airyaiprime(z) = airy(1,z) -@vectorize_1arg Number airyaiprime -airybi(z) = airy(2,z) -@vectorize_1arg Number airybi -airybiprime(z) = airy(3,z) -@vectorize_1arg Number airybiprime - -airy(k::Number, x::FloatingPoint) = oftype(x, real(airy(k, complex(x)))) -airy(k::Number, x::Real) = airy(k, float(x)) -airy(k::Number, z::Complex64) = complex64(airy(k, complex128(z))) -airy(k::Number, z::Complex) = airy(convert(Int,k), complex128(z)) -@vectorize_2arg Number airy - -const cy = Array(Float64,2) -const ae = Array(Int32,2) -const wrk = Array(Float64,2) - -function _besselh(nu::Float64, k::Integer, z::Complex128) - ccall((:zbesh_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &1, &k, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besseli(nu::Float64, z::Complex128) - ccall((:zbesi_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &1, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besselj(nu::Float64, z::Complex128) - ccall((:zbesj_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &1, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besselk(nu::Float64, z::Complex128) - ccall((:zbesk_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &1, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _bessely(nu::Float64, z::Complex128) - ccall((:zbesy_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, - Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), - &real(z), &imag(z), &nu, &1, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(wrk,1), - pointer(wrk,2), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function besselh(nu::Float64, k::Integer, z::Complex128) - if nu < 0 - s = (k == 1) ? 1 : -1 - return _besselh(-nu, k, z) * complex(cospi(nu),-s*sinpi(nu)) - end - return _besselh(nu, k, z) -end - -function besseli(nu::Float64, z::Complex128) - if nu < 0 - return _besseli(-nu,z) - 2_besselk(-nu,z)*sinpi(nu)/pi - else - return _besseli(nu, z) - end -end - -function besselj(nu::Float64, z::Complex128) - if nu < 0 - return _besselj(-nu,z)cos(pi*nu) + _bessely(-nu,z)*sinpi(nu) - else - return _besselj(nu, z) - end -end - -besselj(nu::Integer, x::FloatingPoint) = typemin(Int32) <= nu <= typemax(Int32) ? - oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x)) : - besselj(float64(nu), x) - -besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ? - ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) : - besselj(float64(nu), x) - -besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z) - -function bessely(nu::Float64, z::Complex128) - if nu < 0 - return _bessely(-nu,z)*cospi(nu) - _besselj(-nu,z)*sinpi(nu) - else - return _bessely(nu, z) - end -end - -besselh(nu, z) = besselh(nu, 1, z) -besselh(nu::Real, k::Integer, z::Complex64) = complex64(besselh(float64(nu), k, complex128(z))) -besselh(nu::Real, k::Integer, z::Complex) = besselh(float64(nu), k, complex128(z)) -besselh(nu::Real, k::Integer, x::Real) = besselh(float64(nu), k, complex128(x)) -@vectorize_2arg Number besselh - -besseli(nu::Real, z::Complex64) = complex64(besseli(float64(nu), complex128(z))) -besseli(nu::Real, z::Complex) = besseli(float64(nu), complex128(z)) -besseli(nu::Real, x::Integer) = besseli(nu, float64(x)) -function besseli(nu::Real, x::FloatingPoint) - if x < 0 && !isinteger(nu) - throw(DomainError()) - end - oftype(x, real(besseli(float64(nu), complex128(x)))) -end -@vectorize_2arg Number besseli - -function besselj(nu::FloatingPoint, x::FloatingPoint) - if isinteger(nu) - if typemin(Int32) <= nu <= typemax(Int32) - return besselj(int(nu), x) - end - elseif x < 0 - throw(DomainError()) - end - oftype(x, real(besselj(float64(nu), complex128(x)))) -end - -besselj(nu::Real, z::Complex64) = complex64(besselj(float64(nu), complex128(z))) -besselj(nu::Real, z::Complex) = besselj(float64(nu), complex128(z)) -besselj(nu::Real, x::Integer) = besselj(nu, float(x)) -@vectorize_2arg Number besselj - -besselk(nu::Real, z::Complex64) = complex64(besselk(float64(nu), complex128(z))) -besselk(nu::Real, z::Complex) = besselk(float64(nu), complex128(z)) -besselk(nu::Real, x::Integer) = besselk(nu, float64(x)) -function besselk(nu::Real, x::FloatingPoint) - if x < 0 - throw(DomainError()) - end - if x == 0 - return oftype(x, Inf) - end - oftype(x, real(besselk(float64(nu), complex128(x)))) -end -@vectorize_2arg Number besselk - -bessely(nu::Real, z::Complex64) = complex64(bessely(float64(nu), complex128(z))) -bessely(nu::Real, z::Complex) = bessely(float64(nu), complex128(z)) -bessely(nu::Real, x::Integer) = bessely(nu, float64(x)) -function bessely(nu::Real, x::FloatingPoint) - if x < 0 - throw(DomainError()) - end - if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32) - return bessely(int(nu), x) - end - oftype(x, real(bessely(float64(nu), complex128(x)))) -end -function bessely(nu::Integer, x::FloatingPoint) - if x < 0 - throw(DomainError()) - end - return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x)) -end -function bessely(nu::Integer, x::Float32) - if x < 0 - throw(DomainError()) - end - return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x) -end -@vectorize_2arg Number bessely - -hankelh1(nu, z) = besselh(nu, 1, z) -@vectorize_2arg Number hankelh1 - -hankelh2(nu, z) = besselh(nu, 2, z) -@vectorize_2arg Number hankelh2 - - function angle_restrict_symm(theta) const P1 = 4 * 7.8539812564849853515625e-01 const P2 = 4 * 3.7748947079307981766760e-08 @@ -717,742 +240,6 @@ return r end -const clg_coeff = [76.18009172947146, - -86.50532032941677, - 24.01409824083091, - -1.231739572450155, - 0.1208650973866179e-2, - -0.5395239384953e-5] - -function clgamma_lanczos(z) - const sqrt2pi = 2.5066282746310005 - - y = x = z - temp = x + 5.5 - zz = log(temp) - zz = zz * (x+0.5) - temp -= zz - ser = complex(1.000000000190015, 0) - for j=1:6 - y += 1.0 - zz = clg_coeff[j]/y - ser += zz - end - zz = sqrt2pi*ser / x - return log(zz) - temp -end - -function lgamma(z::Complex) - if real(z) <= 0.5 - a = clgamma_lanczos(1-z) - b = log(sinpi(z)) - const logpi = 1.14472988584940017 - z = logpi - b - a - else - z = clgamma_lanczos(z) - end - complex(real(z), angle_restrict_symm(imag(z))) -end - -gamma(z::Complex) = exp(lgamma(z)) - -# Derivatives of the digamma function -function psifn(x::Float64, n::Int, kode::Int, m::Int) -# Translated from http://www.netlib.org/slatec/src/dpsifn.f -# Note: Underflow handling at 380 in original is skipped - const nmax = 100 - ans = Array(Float64, m) -#----------------------------------------------------------------------- -# bernoulli numbers -#----------------------------------------------------------------------- - const b = [1.00000000000000000e+00, - -5.00000000000000000e-01,1.66666666666666667e-01, - -3.33333333333333333e-02,2.38095238095238095e-02, - -3.33333333333333333e-02,7.57575757575757576e-02, - -2.53113553113553114e-01,1.16666666666666667e+00, - -7.09215686274509804e+00,5.49711779448621554e+01, - -5.29124242424242424e+02,6.19212318840579710e+03, - -8.65802531135531136e+04,1.42551716666666667e+06, - -2.72982310678160920e+07,6.01580873900642368e+08, - -1.51163157670921569e+10,4.29614643061166667e+11, - -1.37116552050883328e+13,4.88332318973593167e+14, - -1.92965793419400681e+16] - trm = Array(Float64, 22) - trmr = Array(Float64, 100) -#***first executable statement dpsifn - if x <= 0.0 throw(DomainError()) end - if n < 0 error("n must be non-negative") end - if kode < 1 | kode > 2 error("kode must be one or two") end - if m < 1 error("m must be larger than one") end - mm = m - const nx = min(-exponent(realmin(Float64)) + 1, exponent(realmax(Float64))) - const r1m5 = log10(2) - const r1m4 = Base.eps(Float64) * 0.5 - const wdtol = max(r1m4, 0.5e-18) -#----------------------------------------------------------------------- -# elim = approximate exponential over and underflow limit -#----------------------------------------------------------------------- - const elim = 2.302*(nx*r1m5 - 3.0) - xln = log(x) - nn = n + mm - 1 - fn = nn - t = (fn + 1)*xln -#----------------------------------------------------------------------- -# overflow and underflow test for small and large x -#----------------------------------------------------------------------- - if abs(t) > elim - if t <= 0.0 error("n too large") end - error("overflow; x too small or n+m-1 too large or both") - end - if x < wdtol - ans[1] = x^(-n - 1) - if mm != 1 - k = 1 - for i = 2:mm - ans[k + 1] = ans[k]/x - k += 1 - end - end - if n != 0 return ans end - if kode == 2 ans[1] = ans[1] + xln end - return ans - end -#----------------------------------------------------------------------- -# compute xmin and the number of terms of the series, fln+1 -#----------------------------------------------------------------------- - rln = r1m5 * precision(x) - rln = min(rln, 18.06) - fln = max(rln, 3.0) - 3.0 - yint = 3.50 + 0.40*fln - slope = 0.21 + fln*(0.0006038*fln + 0.008677) - xm = yint + slope*fn - mx = itrunc(xm) + 1 - xmin = mx - if n != 0 - xm = -2.302*rln - min(0.0,xln) - arg = xm/n - arg = min(0.0,arg) - eps = exp(arg) - xm = 1.0 - eps - if abs(arg) < 1.0e-3 xm = -arg end - fln = x*xm/eps - xm = xmin - x - if (xm > 7.0) & (fln < 15.0) - nn = itrunc(fln) + 1 - np = n + 1 - t1 = (n + 1)*xln - t = exp(-t1) - s = t - den = x - for i = 1:nn - den += 1.0 - trm[i] = den^(-np) - s += trm[i] - end - ans[1] = s - if n == 0 - if kode == 2 ans[1] = s + xln end - end - if mm == 1 return ans end -#----------------------------------------------------------------------- -# generate higher derivatives, j.gt.n -#----------------------------------------------------------------------- - tol = wdtol/5.0 - for j = 2:mm - t = t/x - s = t - tols = t*tol - den = x - for i = 1:nn - den += 1.0 - trm[i] = trm[i]/den - s += trm[i] - if trm[i] < tols break end - end - ans[j] = s - end - return ans - end - end - - xdmy = x - xdmln = xln - xinc = 0.0 - if x < xmin - nx = itrunc(x) - xinc = xmin - nx - xdmy = x + xinc - xdmln = log(xdmy) - end -#----------------------------------------------------------------------- -# generate w(n+mm-1,x) by the asymptotic expansion -#----------------------------------------------------------------------- - t = fn*xdmln - t1 = xdmln + xdmln - t2 = t + xdmln - tk = max(abs(t), abs(t1), abs(t2)) - if tk > elim error("underflow") end - tss = exp(-t) - tt = 0.5/xdmy - t1 = tt - tst = wdtol*tt - if nn != 0 t1 = tt + 1.0/fn end - rxsq = 1.0/(xdmy*xdmy) - ta = 0.5*rxsq - t = (fn + 1)*ta - s = t*b[3] - if abs(s) >= tst - tk = 2.0 - for k = 4:22 - t = t*((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0))*rxsq - trm[k] = t*b[k] - if abs(trm[k]) < tst break end - s += trm[k] - tk += 2.0 - end - end - s = (s + t1)*tss - while true - if xinc != 0.0 -#----------------------------------------------------------------------- -# backward recur from xdmy to x -#----------------------------------------------------------------------- - nx = itrunc(xinc) - np = nn + 1 - if nx > nmax error("n too large") end - if nn == 0 break end - xm = xinc - 1.0 - fx = x + xm -#----------------------------------------------------------------------- -# this loop should not be changed. fx is accurate when x is small -#----------------------------------------------------------------------- - for i = 1:nx - trmr[i] = fx^(-np) - s += trmr[i] - xm -= 1.0 - fx = x + xm - end - end - ans[mm] = s - if fn == 0 - if kode != 2 - ans[1] = s - xdmln - return ans - end - if xdmy == x return ans end - xq = xdmy/x - ans[1] = s - log(xq) - return ans - end -#----------------------------------------------------------------------- -# generate lower derivatives, j.lt.n+mm-1 -#----------------------------------------------------------------------- - if mm == 1 return ans end - for j = 2:mm - fn -= 1 - tss *= xdmy - t1 = tt - if fn != 0 t1 = tt + 1.0/fn end - t = (fn + 1)*ta - s = t*b[3] - if abs(s) >= tst - tk = 4 + fn - for k = 4:22 #110 - trm[k] = trm[k]*(fn + 1)/tk - if abs(trm[k]) < tst break end - s += trm[k] - tk += 2.0 - end - end - s = (s + t1)*tss - if xinc != 0.0 - if fn == 0 break end - xm = xinc - 1.0 - fx = x + xm - for i = 1:nx - trmr[i] = trmr[i]*fx - s += trmr[i] - xm -= 1.0 - fx = x + xm - end - end - mx = mm - j + 1 - ans[mx] = s - if fn == 0 - if kode != 2 - ans[1] = s - xdmln - return ans - end - if xdmy == x return ans end - xq = xdmy/x - ans[1] = s - log(xq) - return ans - end - end - if fn == 0 break end - return ans - end -#----------------------------------------------------------------------- -# recursion for n = 0 -#----------------------------------------------------------------------- - for i = 1:nx - s += 1.0/(x + nx - i) - end - if kode != 2 - ans[1] = s - xdmln - return ans - end - if xdmy == x return ans end - xq = xdmy/x - ans[1] = s - log(xq) - return ans -end -polygamma(k::Int, x::Float64) = (2rem(k,2) - 1)*psifn(x, k, 1, 1)[1]/gamma(k + 1) -polygamma(k::Int, x::Float32) = float32(polygamma(k, float64(x))) -polygamma(k::Int, x::Real) = polygamma(k, float64(x)) - -# Translation of psi.c from cephes -function digamma(x::Float64) - negative = false - nz = 0.0 - - if x <= 0.0 - negative = true - q = x - p = floor(q) - if p == q - return NaN - end - - nz = q - p - if nz != 0.5 - if nz > 0.5 - p += 1.0 - nz = q - p - end - nz = pi / tan(pi * nz) - else - nz = 0.0 - end - x = 1.0 - x - end - - if x <= 10.0 && x == floor(x) - y = 0.0 - for i = 1:x-1 - y += 1.0 / i - end - y -= γ # γ == -digamma(1) == 0.5772156649015328606065121; - - if negative - y -= nz - end - return y - end - - w = 0.0 - while x < 10.0 - w += 1.0 / x - x += 1.0 - end - - if x < 1.0e17 - z = 1.0 / (x*x) - y = @horner(z, 8.33333333333333333333e-2, -8.33333333333333333333e-3, 3.96825396825396825397e-3, - -4.16666666666666666667e-3, 7.57575757575757575758e-3,-2.10927960927960927961e-2, - 8.33333333333333333333e-2) - y *= z - else - y = 0.0 - end - - y = log(x) - 0.5/x - y - w - - if negative - y -= nz - end - - return y -end -digamma(x::Float32) = float32(digamma(float64(x))) -digamma(x::Real) = digamma(float64(x)) -@vectorize_1arg Real digamma - -trigamma(x::Real) = polygamma(1, x) -@vectorize_1arg Real trigamma - -# Inverse digamma function -# -# Implementation of fixed point algorithm described in -# "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000 -function invdigamma(y::Float64) - # Closed form initial estimates - if y >= -2.22 - x_old = exp(y) + 0.5 - x_new = x_old - else - x_old = -1.0 / (y - digamma(1.0)) - x_new = x_old - end - - # Fixed point algorithm - delta = Inf - iteration = 0 - while delta > 1e-12 && iteration < 25 - iteration += 1 - x_new = x_old - (digamma(x_old) - y) / trigamma(x_old) - delta = abs(x_new - x_old) - x_old = x_new - end - - return x_new -end -invdigamma(x::Float32) = float32(invdigamma(float64(x))) -invdigamma(x::Real) = invdigamma(float64(x)) -@vectorize_1arg Real invdigamma - -function beta(x::Number, w::Number) - yx, sx = lgamma_r(x) - yw, sw = lgamma_r(w) - yxw, sxw = lgamma_r(x+w) - return copysign(exp(yx + yw - yxw), sx*sw*sxw) -end -lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w) -@vectorize_2arg Number beta -@vectorize_2arg Number lbeta - -const eta_coeffs = - [.99999999999999999997, - -.99999999999999999821, - .99999999999999994183, - -.99999999999999875788, - .99999999999998040668, - -.99999999999975652196, - .99999999999751767484, - -.99999999997864739190, - .99999999984183784058, - -.99999999897537734890, - .99999999412319859549, - -.99999996986230482845, - .99999986068828287678, - -.99999941559419338151, - .99999776238757525623, - -.99999214148507363026, - .99997457616475604912, - -.99992394671207596228, - .99978893483826239739, - -.99945495809777621055, - .99868681159465798081, - -.99704078337369034566, - .99374872693175507536, - -.98759401271422391785, - .97682326283354439220, - -.95915923302922997013, - .93198380256105393618, - -.89273040299591077603, - .83945793215750220154, - -.77148960729470505477, - .68992761745934847866, - -.59784149990330073143, - .50000000000000000000, - -.40215850009669926857, - .31007238254065152134, - -.22851039270529494523, - .16054206784249779846, - -.10726959700408922397, - .68016197438946063823e-1, - -.40840766970770029873e-1, - .23176737166455607805e-1, - -.12405987285776082154e-1, - .62512730682449246388e-2, - -.29592166263096543401e-2, - .13131884053420191908e-2, - -.54504190222378945440e-3, - .21106516173760261250e-3, - -.76053287924037718971e-4, - .25423835243950883896e-4, - -.78585149263697370338e-5, - .22376124247437700378e-5, - -.58440580661848562719e-6, - .13931171712321674741e-6, - -.30137695171547022183e-7, - .58768014045093054654e-8, - -.10246226511017621219e-8, - .15816215942184366772e-9, - -.21352608103961806529e-10, - .24823251635643084345e-11, - -.24347803504257137241e-12, - .19593322190397666205e-13, - -.12421162189080181548e-14, - .58167446553847312884e-16, - -.17889335846010823161e-17, - .27105054312137610850e-19] - -function eta(z::Union(Float64,Complex128)) - if z == 0 - return oftype(z, 0.5) - end - re, im = reim(z) - if im==0 && re < 0 && re==round(re/2)*2 - return zero(z) - end - reflect = false - if re < 0.5 - z = 1-z - reflect = true - end - s = zero(z) - for n = length(eta_coeffs):-1:1 - c = eta_coeffs[n] - p = n^-z - s += c * p - end - if reflect - z2 = 2.0^z - b = 2.0 - (2.0*z2) - f = z2 - 2 - piz = pi^z - - b = b/f/piz - - return s * gamma(z) * b * cospi(z/2) - end - return s -end - -eta(x::Integer) = eta(float64(x)) -eta(x::Real) = oftype(x,eta(float64(x))) -eta(z::Complex) = oftype(z,eta(complex128(z))) -@vectorize_1arg Number eta - -function zeta(z::Number) - zz = 2^z - eta(z) * zz/(zz-2) -end -@vectorize_1arg Number zeta - -@unix_only if WORD_SIZE == 64 -# TODO: complex return only on 64-bit for now -for f in (:erf, :erfc, :erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f - @eval begin - ($fname)(z::Complex128) = complex128(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) - ($fname)(z::Complex64) = complex64(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), complex128(z), float64(eps(Float32)))) - ($fname)(z::Complex) = ($fname)(complex128(z)) - end -end -end -for f in (:erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f - @eval begin - ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), x) - ($fname)(x::Float32) = float32(ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), float64(x))) - ($fname)(x::Integer) = ($fname)(float(x)) - @vectorize_1arg Number $fname - end -end - -# Compute the inverse of the error function: erf(erfinv(x)) == x, -# using the rational approximants tabulated in: -# J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev -# approximations for the inverse of the error function," Math. Comp. 30, -# pp. 827--830 (1976). -# http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7 -# http://www.jstor.org/stable/2005402 -function erfinv(x::Float64) - a = abs(x) - if a >= 1.0 - if x == 1.0 - return inf(Float64) - elseif x == -1.0 - return -inf(Float64) - end - throw(DomainError()) - elseif a <= 0.75 # Table 17 in Blair et al. - t = x*x - 0.5625 - return x * @horner(t, 0.16030_49558_44066_229311e2, - -0.90784_95926_29603_26650e2, - 0.18644_91486_16209_87391e3, - -0.16900_14273_46423_82420e3, - 0.65454_66284_79448_7048e2, - -0.86421_30115_87247_794e1, - 0.17605_87821_39059_0) / - @horner(t, 0.14780_64707_15138_316110e2, - -0.91374_16702_42603_13936e2, - 0.21015_79048_62053_17714e3, - -0.22210_25412_18551_32366e3, - 0.10760_45391_60551_23830e3, - -0.20601_07303_28265_443e2, - 0.1e1) - elseif a <= 0.9375 # Table 37 in Blair et al. - t = x*x - 0.87890625 - return x * @horner(t, -0.15238_92634_40726_128e-1, - 0.34445_56924_13612_5216, - -0.29344_39867_25424_78687e1, - 0.11763_50570_52178_27302e2, - -0.22655_29282_31011_04193e2, - 0.19121_33439_65803_30163e2, - -0.54789_27619_59831_8769e1, - 0.23751_66890_24448) / - @horner(t, -0.10846_51696_02059_954e-1, - 0.26106_28885_84307_8511, - -0.24068_31810_43937_57995e1, - 0.10695_12997_33870_14469e2, - -0.23716_71552_15965_81025e2, - 0.24640_15894_39172_84883e2, - -0.10014_37634_97830_70835e2, - 0.1e1) - else # Table 57 in Blair et al. - t = 1.0 / sqrt(-log(1.0 - a)) - return @horner(t, 0.10501_31152_37334_38116e-3, - 0.10532_61131_42333_38164_25e-1, - 0.26987_80273_62432_83544_516, - 0.23268_69578_89196_90806_414e1, - 0.71678_54794_91079_96810_001e1, - 0.85475_61182_21678_27825_185e1, - 0.68738_08807_35438_39802_913e1, - 0.36270_02483_09587_08930_02e1, - 0.88606_27392_96515_46814_9) / - (copysign(t, x) * - @horner(t, 0.10501_26668_70303_37690e-3, - 0.10532_86230_09333_27531_11e-1, - 0.27019_86237_37515_54845_553, - 0.23501_43639_79702_53259_123e1, - 0.76078_02878_58012_77064_351e1, - 0.11181_58610_40569_07827_3451e2, - 0.11948_78791_84353_96667_8438e2, - 0.81922_40974_72699_07893_913e1, - 0.40993_87907_63680_15361_45e1, - 0.1e1)) - end -end - -function erfinv(x::Float32) - a = abs(x) - if a >= 1.0f0 - if x == 1.0f0 - return inf(Float32) - elseif x == -1.0f0 - return -inf(Float32) - end - throw(DomainError()) - elseif a <= 0.75f0 # Table 10 in Blair et al. - t = x*x - 0.5625f0 - return x * @horner(t, -0.13095_99674_22f2, - 0.26785_22576_0f2, - -0.92890_57365f1) / - @horner(t, -0.12074_94262_97f2, - 0.30960_61452_9f2, - -0.17149_97799_1f2, - 0.1f1) - elseif a <= 0.9375f0 # Table 29 in Blair et al. - t = x*x - 0.87890625f0 - return x * @horner(t, -0.12402_56522_1f0, - 0.10688_05957_4f1, - -0.19594_55607_8f1, - 0.42305_81357f0) / - @horner(t, -0.88276_97997f-1, - 0.89007_43359f0, - -0.21757_03119_6f1, - 0.1f1) - else # Table 50 in Blair et al. - t = 1.0f0 / sqrt(-log(1.0f0 - a)) - return @horner(t, 0.15504_70003_116f0, - 0.13827_19649_631f1, - 0.69096_93488_87f0, - -0.11280_81391_617f1, - 0.68054_42468_25f0, - -0.16444_15679_1f0) / - (copysign(t, x) * - @horner(t, 0.15502_48498_22f0, - 0.13852_28141_995f1, - 0.1f1)) - end -end - -erfinv(x::Integer) = erfinv(float(x)) -@vectorize_1arg Real erfinv - -# Inverse complementary error function: use Blair tables for y = 1-x, -# exploiting the greater accuracy of y (vs. x) when y is small. -function erfcinv(y::Float64) - if y > 0.0625 - return erfinv(1.0 - y) - elseif y <= 0.0 - if y == 0.0 - return inf(Float64) - end - throw(DomainError()) - elseif y >= 1e-100 # Table 57 in Blair et al. - t = 1.0 / sqrt(-log(y)) - return @horner(t, 0.10501_31152_37334_38116e-3, - 0.10532_61131_42333_38164_25e-1, - 0.26987_80273_62432_83544_516, - 0.23268_69578_89196_90806_414e1, - 0.71678_54794_91079_96810_001e1, - 0.85475_61182_21678_27825_185e1, - 0.68738_08807_35438_39802_913e1, - 0.36270_02483_09587_08930_02e1, - 0.88606_27392_96515_46814_9) / - (t * - @horner(t, 0.10501_26668_70303_37690e-3, - 0.10532_86230_09333_27531_11e-1, - 0.27019_86237_37515_54845_553, - 0.23501_43639_79702_53259_123e1, - 0.76078_02878_58012_77064_351e1, - 0.11181_58610_40569_07827_3451e2, - 0.11948_78791_84353_96667_8438e2, - 0.81922_40974_72699_07893_913e1, - 0.40993_87907_63680_15361_45e1, - 0.1e1)) - else # Table 80 in Blair et al. - t = 1.0 / sqrt(-log(y)) - return @horner(t, 0.34654_29858_80863_50177e-9, - 0.25084_67920_24075_70520_55e-6, - 0.47378_13196_37286_02986_534e-4, - 0.31312_60375_97786_96408_3388e-2, - 0.77948_76454_41435_36994_854e-1, - 0.70045_68123_35816_43868_271e0, - 0.18710_42034_21679_31668_683e1, - 0.71452_54774_31351_45428_3e0) / - (t * @horner(t, 0.34654_29567_31595_11156e-9, - 0.25084_69079_75880_27114_87e-6, - 0.47379_53129_59749_13536_339e-4, - 0.31320_63536_46177_68848_0813e-2, - 0.78073_48906_27648_97214_733e-1, - 0.70715_04479_95337_58619_993e0, - 0.19998_51543_49112_15105_214e1, - 0.15072_90269_27316_80008_56e1, - 0.1e1)) - end -end - -function erfcinv(y::Float32) - if y > 0.0625f0 - return erfinv(1.0f0 - y) - elseif y <= 0.0f0 - if y == 0.0f0 - return inf(Float32) - end - throw(DomainError()) - else # Table 50 in Blair et al. - t = 1.0f0 / sqrt(-log(y)) - return @horner(t, 0.15504_70003_116f0, - 0.13827_19649_631f1, - 0.69096_93488_87f0, - -0.11280_81391_617f1, - 0.68054_42468_25f0, - -0.16444_15679_1f0) / - (t * - @horner(t, 0.15502_48498_22f0, - 0.13852_28141_995f1, - 0.1f1)) - end -end - -erfcinv(x::Integer) = erfcinv(float(x)) -@vectorize_1arg Real erfcinv - ## mod2pi-related calculations ## function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64) @@ -1539,4 +326,11 @@ mod2pi(fx) end +# More special functions + +include("special/trig.jl") +include("special/bessel.jl") +include("special/erf.jl") +include("special/gamma.jl") + end # module diff -Nru julia-0.3.0/base/pcre.jl julia-0.3.0/base/pcre.jl --- julia-0.3.0/base/pcre.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/pcre.jl 2014-06-03 19:33:18.000000000 +0000 @@ -99,6 +99,7 @@ end extra end + study(re::Ptr{Void}) = study(re, int32(0)) free_study(extra::Ptr{Void}) = @@ -106,17 +107,24 @@ free(regex::Ptr{Void}) = ccall(unsafe_load(cglobal((:pcre_free, :libpcre),Ptr{Void})), Void, (Ptr{Void},), regex) -exec(regex::Ptr{Void}, extra::Ptr{Void}, str::SubString, offset::Integer, options::Integer, cap::Bool) = - exec(regex, extra, str.string, str.offset, offset, sizeof(str), options, cap) -exec(regex::Ptr{Void}, extra::Ptr{Void}, str::ByteString, offset::Integer, options::Integer, cap::Bool) = - exec(regex, extra, str, 0, offset, sizeof(str), options, cap) +function exec(regex::Ptr{Void}, extra::Ptr{Void}, str::SubString, offset::Integer, + options::Integer, ovec::Vector{Int32}) + return exec(regex, extra, str.string, str.offset, offset, sizeof(str), + options, ovec) +end + +function exec(regex::Ptr{Void}, extra::Ptr{Void}, str::ByteString, offset::Integer, + options::Integer, ovec::Vector{Int32}) + return exec(regex, extra, str, 0, offset, sizeof(str), options, ovec) +end + function exec(regex::Ptr{Void}, extra::Ptr{Void}, - str::ByteString, shift::Integer, offset::Integer, len::Integer, options::Integer, cap::Bool) + str::ByteString, shift::Integer, offset::Integer, + len::Integer, options::Integer, + ovec::Vector{Int32}) if offset < 0 || len < offset || len+shift > sizeof(str) error(BoundsError) end - ncap = info(regex, extra, INFO_CAPTURECOUNT, Int32) - ovec = Array(Int32, 3(ncap+1)) n = ccall((:pcre_exec, :libpcre), Int32, (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32, Int32, Int32, Ptr{Int32}, Int32), @@ -125,7 +133,7 @@ if n < -1 error("error $n") end - cap ? ((n > -1 ? ovec[1:2(ncap+1)] : Array(Int32,0)), ncap) : n > -1 + return n > -1 end end # module diff -Nru julia-0.3.0/base/printf.jl julia-0.3.0/base/printf.jl --- julia-0.3.0/base/printf.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/printf.jl 2014-06-03 19:33:18.000000000 +0000 @@ -794,7 +794,7 @@ write(fmt, '0') elseif precision > 0 write(fmt, '.') - print(fmt, precision+1) + print(fmt, precision) end write(fmt, 'R') write(fmt, c) diff -Nru julia-0.3.0/base/range.jl julia-0.3.0/base/range.jl --- julia-0.3.0/base/range.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/range.jl 2014-06-03 19:33:18.000000000 +0000 @@ -44,7 +44,7 @@ remain = oftype(T, unsigned(diff) % step) end else - remain = diff % step + remain = steprem(start,stop,step) end last = stop - remain end @@ -54,6 +54,8 @@ end end +steprem(start,stop,step) = (stop-start) % step + StepRange{T,S}(start::T, step::S, stop::T) = StepRange{T,S}(start, step, stop) immutable UnitRange{T<:Real} <: OrdinalRange{T,Int} @@ -491,15 +493,15 @@ ## sorting ## issorted(r::UnitRange) = true -issorted(r::Range) = step(r) >= 0 +issorted(r::Range) = step(r) >= zero(step(r)) sort(r::UnitRange) = r sort!(r::UnitRange) = r -sort{T<:Real}(r::Range{T}) = issorted(r) ? r : reverse(r) +sort(r::Range) = issorted(r) ? r : reverse(r) sortperm(r::UnitRange) = 1:length(r) -sortperm{T<:Real}(r::Range{T}) = issorted(r) ? (1:1:length(r)) : (length(r):-1:1) +sortperm(r::Range) = issorted(r) ? (1:1:length(r)) : (length(r):-1:1) function sum{T<:Real}(r::Range{T}) l = length(r) @@ -515,7 +517,7 @@ end function in(x, r::Range) - n = step(r) == 0 ? 1 : iround((x-first(r))/step(r))+1 + n = step(r) == zero(step(r)) ? 1 : iround((x-first(r))/step(r))+1 n >= 1 && n <= length(r) && r[n] == x end diff -Nru julia-0.3.0/base/reduce.jl julia-0.3.0/base/reduce.jl --- julia-0.3.0/base/reduce.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/reduce.jl 2014-06-03 19:33:18.000000000 +0000 @@ -18,6 +18,8 @@ type MulFun <: Func{2} end type AndFun <: Func{2} end type OrFun <: Func{2} end +type MaxFun <: Func{2} end +type MinFun <: Func{2} end evaluate(::IdFun, x) = x evaluate(::AbsFun, x) = abs(x) @@ -30,239 +32,161 @@ evaluate(::MulFun, x, y) = x * y evaluate(::AndFun, x, y) = x & y evaluate(::OrFun, x, y) = x | y +evaluate(::MaxFun, x, y) = scalarmax(x, y) +evaluate(::MinFun, x, y) = scalarmin(x, y) evaluate(f::Callable, x, y) = f(x, y) -###### Generic reduction functions ###### -# Note that getting type-stable results from reduction functions, -# or at least having type-stable loops, is nontrivial (#6069). +###### Generic (map)reduce functions ###### -## foldl +# r_promote: promote x to the type of reduce(op, [x]) +r_promote(op, x) = x +r_promote(::AddFun, x) = x + zero(x) +r_promote(::MulFun, x) = x * one(x) -function _foldl(op, v0, itr, i) + +## foldl && mapfoldl + +function mapfoldl_impl(f, op, v0, itr, i) if done(itr, i) return v0 else (x, i) = next(itr, i) - v = evaluate(op, v0, x) + v = evaluate(op, v0, evaluate(f, x)) while !done(itr, i) (x, i) = next(itr, i) - v = evaluate(op, v, x) + v = evaluate(op, v, evaluate(f, x)) end return v end end -function foldl(op::Callable, v0, itr, i) - is(op, +) && return _foldl(AddFun(), v0, itr, i) - is(op, *) && return _foldl(MulFun(), v0, itr, i) - is(op, &) && return _foldl(AndFun(), v0, itr, i) - is(op, |) && return _foldl(OrFun(), v0, itr, i) - return _foldl(op, v0, itr, i) -end +mapfoldl(f, op, v0, itr) = mapfoldl_impl(f, op, v0, itr, start(itr)) -foldl(op::Callable, v0, itr) = foldl(op, v0, itr, start(itr)) +function mapfoldl(f, op::Function, v0, itr) + is(op, +) ? mapfoldl(f, AddFun(), v0, itr) : + is(op, *) ? mapfoldl(f, MulFun(), v0, itr) : + is(op, &) ? mapfoldl(f, AndFun(), v0, itr) : + is(op, |) ? mapfoldl(f, OrFun(), v0, itr) : + mapfoldl_impl(f, op, v0, itr, start(itr)) +end -function foldl(op::Callable, itr) +function mapfoldl(f, op, itr) i = start(itr) done(itr, i) && error("Argument is empty.") - (v0, i) = next(itr, i) - return foldl(op, v0, itr, i) + (x, i) = next(itr, i) + v0 = evaluate(f, x) + mapfoldl_impl(f, op, v0, itr, i) end -## foldr +foldl(op, v0, itr) = mapfoldl(IdFun(), op, v0, itr) +foldl(op, itr) = mapfoldl(IdFun(), op, itr) -function foldr(op::Callable, v0, itr, i=endof(itr)) - # use type stable procedure +## foldr & mapfoldr + +function mapfoldr_impl(f, op, v0, itr, i::Integer) if i == 0 return v0 else - v = op(itr[i], v0) + x = itr[i] + v = evaluate(op, evaluate(f, x), v0) while i > 1 x = itr[i -= 1] - v = op(x, v) + v = evaluate(op, evaluate(f, x), v) end return v end end -foldr(op::Callable, itr) = (i = endof(itr); foldr(op, itr[i], itr, i-1)) +mapfoldr(f, op, v0, itr) = mapfoldr_impl(f, op, v0, itr, endof(itr)) +mapfoldr(f, op, itr) = (i = endof(itr); mapfoldr_impl(f, op, evaluate(f, itr[i]), itr, i-1)) -## reduce +foldr(op, v0, itr) = mapfoldr(IdFun(), op, v0, itr) +foldr(op, itr) = mapfoldr(IdFun(), op, itr) -reduce(op::Callable, v, itr) = foldl(op, v, itr) +## reduce & mapreduce -function reduce(op::Callable, itr) # this is a left fold - if is(op, +) - return sum(itr) - elseif is(op, *) - return prod(itr) - elseif is(op, |) - return any(itr) - elseif is(op, &) - return all(itr) - end - return foldl(op, itr) -end - -# pairwise reduction, requires n > 1 (to allow type-stable loop) -function r_pairwise(op::Callable, A::AbstractArray, i1,n) - if n < 128 - @inbounds v = op(A[i1], A[i1+1]) - for i = i1+2:i1+n-1 - @inbounds v = op(v,A[i]) - end - return v - else - n2 = div(n,2) - return op(r_pairwise(op,A, i1,n2), r_pairwise(op,A, i1+n2,n-n2)) +# mapreduce_***_impl require ifirst < ilast +function mapreduce_seq_impl(f, op, A::AbstractArray, ifirst::Int, ilast::Int) + @inbounds fx1 = evaluate(f, A[ifirst]) + @inbounds fx2 = evaluate(f, A[ifirst+=1]) + @inbounds v = evaluate(op, fx1, fx2) + while ifirst < ilast + @inbounds fx = evaluate(f, A[ifirst+=1]) + v = evaluate(op, v, fx) end + return v end -function reduce(op::Callable, A::AbstractArray) - n = length(A) - n == 0 ? error("argument is empty") : n == 1 ? A[1] : r_pairwise(op,A, 1,n) -end - -function reduce(op::Callable, v0, A::AbstractArray) - n = length(A) - n == 0 ? v0 : n == 1 ? op(v0, A[1]) : op(v0, r_pairwise(op,A, 1,n)) -end - - -###### Specific reduction functions ###### - -## in & contains - -function in(x, itr) - for y in itr - if y == x - return true - end +function mapreduce_pairwise_impl(f, op, A::AbstractArray, ifirst::Int, ilast::Int, blksize::Int) + if ifirst + blksize > ilast + return mapreduce_seq_impl(f, op, A, ifirst, ilast) + else + imid = (ifirst + ilast) >>> 1 + v1 = mapreduce_seq_impl(f, op, A, ifirst, imid) + v2 = mapreduce_seq_impl(f, op, A, imid+1, ilast) + return evaluate(op, v1, v2) end - return false end -const ∈ = in -∉(x, itr)=!∈(x, itr) -∋(itr, x)= ∈(x, itr) -∌(itr, x)=!∋(itr, x) -function contains(itr, x) - depwarn("contains(collection, item) is deprecated, use in(item, collection) instead", :contains) - in(x, itr) -end +mapreduce(f, op, itr) = mapfoldl(f, op, itr) +mapreduce(f, op, v0, itr) = mapfoldl(f, op, v0, itr) +mapreduce_impl(f, op, A::AbstractArray, ifirst::Int, ilast::Int) = + mapreduce_seq_impl(f, op, A, ifirst, ilast) + +# handling empty arrays +mr_empty(f, op, T) = error("Reducing over an empty array is not allow.") +mr_empty(::IdFun, op::AddFun, T) = r_promote(op, zero(T)) +mr_empty(::AbsFun, op::AddFun, T) = r_promote(op, abs(zero(T))) +mr_empty(::Abs2Fun, op::AddFun, T) = r_promote(op, abs2(zero(T))) +mr_empty(::IdFun, op::MulFun, T) = r_promote(op, one(T)) +mr_empty(::AbsFun, op::MaxFun, T) = abs(zero(T)) +mr_empty(::Abs2Fun, op::MaxFun, T) = abs2(zero(T)) +mr_empty(f, op::AndFun, T) = true +mr_empty(f, op::OrFun, T) = false -function contains(eq::Function, itr, x) - for y in itr - if eq(y, x) - return true +function _mapreduce{T}(f, op, A::AbstractArray{T}) + n = length(A) + if n == 0 + return mr_empty(f, op, T) + elseif n == 1 + return r_promote(op, evaluate(f, A[1])) + elseif n < 16 + @inbounds fx1 = evaluate(f, A[1]) + @inbounds fx2 = evaluate(f, A[2]) + s = evaluate(op, fx1, fx2) + i = 2 + while i < n + @inbounds fx = evaluate(f, A[i+=1]) + s = evaluate(op, s, fx) end + return s + else + return mapreduce_impl(f, op, A, 1, n) end - return false end +mapreduce(f, op, A::AbstractArray) = _mapreduce(f, op, A) +mapreduce(f, op, a::Number) = evaluate(f, a) -## countnz & count - -function countnz(itr) - n = 0 - for x in itr - if x != 0 - n += 1 - end - end - return n +function mapreduce(f, op::Function, A::AbstractArray) + is(op, +) ? _mapreduce(f, AddFun(), A) : + is(op, *) ? _mapreduce(f, MulFun(), A) : + is(op, &) ? _mapreduce(f, AndFun(), A) : + is(op, |) ? _mapreduce(f, OrFun(), A) : + _mapreduce(f, op, A) end -function countnz(a::AbstractArray) - n = 0 - for i = 1:length(a) - @inbounds x = a[i] - if x != 0 - n += 1 - end - end - return n -end +reduce(op, v0, itr) = mapreduce(IdFun(), op, v0, itr) +reduce(op, itr) = mapreduce(IdFun(), op, itr) +reduce(op, a::Number) = a -function count(pred::Function, itr) - n = 0 - for x in itr - if pred(x) - n += 1 - end - end - return n -end +###### Specific reduction functions ###### ## sum -# result type inference for sum - -sumtype{T}(::Type{T}) = typeof(zero(T) + zero(T)) -sumzero{T}(::Type{T}) = zero(T) + zero(T) -addzero(x) = x + zero(x) - -typealias SumResultNumber Union(Uint,Uint64,Uint128,Int,Int64,Int128,Float32,Float64,Complex64,Complex128) - -sumtype{T<:SumResultNumber}(::Type{T}) = T -sumzero{T<:SumResultNumber}(::Type{T}) = zero(T) -addzero(x::SumResultNumber) = x - -sumzero{T<:AbstractArray}(::Type{T}) = error("Summing over an empty collection of arrays is not allowed.") -addzero(a::AbstractArray) = a - -# general sum over iterables - -function _sum(f, itr, s) # deal with non-empty cases - # pre-condition: s = start(itr) && !done(itr, s) - (v, s) = next(itr, s) - done(itr, s) && return addzero(evaluate(f, v)) # adding zero for type stability - # specialize for length > 1 to have type-stable loop - (x, s) = next(itr, s) - result = evaluate(f, v) + evaluate(f, x) - while !done(itr, s) - (x, s) = next(itr, s) - result += evaluate(f, x) - end - return result -end - -function sum(itr) - s = start(itr) - if done(itr, s) - if applicable(eltype, itr) - return sumzero(eltype(itr)) - else - throw(ArgumentError("sum(itr) is undefined for empty collections; instead, do isempty(itr) ? z : sum(itr), where z is the correct type of zero for your sum")) - end - end - _sum(IdFun(), itr, s) -end - -function sum(f::Union(Function,Func{1}), itr) - s = start(itr) - done(itr, s) && error("Argument is empty.") - _sum(f, itr, s) -end - -sum(x::Number) = x -sum(A::AbstractArray{Bool}) = countnz(A) - -sumabs(itr) = sum(AbsFun(), itr) -sumabs2(itr) = sum(Abs2Fun(), itr) - -sumabs(x::Number) = abs(x) -sumabs2(x::Number) = abs2(x) - -# Note: sum_seq uses four accumulators, so each accumulator gets at most 256 numbers -sum_pairwise_blocksize(f) = 1024 - -# a fast implementation of sum in sequential order (from left to right). -# to allow type-stable loops, requires length > 1 -function sum_seq(f, a::AbstractArray, ifirst::Int, ilast::Int) - +function mapreduce_seq_impl(f, op::AddFun, a::AbstractArray, ifirst::Int, ilast::Int) @inbounds if ifirst + 6 >= ilast # length(a) < 8 i = ifirst s = evaluate(f, a[i]) + evaluate(f, a[i+1]) @@ -272,18 +196,11 @@ end return s - else # length(a) >= 8 - - # more effective utilization of the instruction - # pipeline through manually unrolling the sum - # into four-way accumulation. Benchmark shows - # that this results in about 2x speed-up. - + else # length(a) >= 8, manual unrolling s1 = evaluate(f, a[ifirst]) + evaluate(f, a[ifirst + 4]) s2 = evaluate(f, a[ifirst + 1]) + evaluate(f, a[ifirst + 5]) s3 = evaluate(f, a[ifirst + 2]) + evaluate(f, a[ifirst + 6]) s4 = evaluate(f, a[ifirst + 3]) + evaluate(f, a[ifirst + 7]) - i = ifirst + 8 il = ilast - 3 while i <= il @@ -293,65 +210,28 @@ s4 += evaluate(f, a[i+3]) i += 4 end - while i <= ilast s1 += evaluate(f, a[i]) i += 1 end - return s1 + s2 + s3 + s4 - end + end end -# Pairwise (cascade) summation of A[i1:i1+n-1], which has O(log n) error growth -# [vs O(n) for a simple loop] with negligible performance cost if -# the base case is large enough. See, e.g.: -# http://en.wikipedia.org/wiki/Pairwise_summation -# Higham, Nicholas J. (1993), "The accuracy of floating point -# summation", SIAM Journal on Scientific Computing 14 (4): 783–799. -# In fact, the root-mean-square error growth, assuming random roundoff -# errors, is only O(sqrt(log n)), which is nearly indistinguishable from O(1) -# in practice. See: -# Manfred Tasche and Hansmartin Zeuner, Handbook of -# Analytic-Computational Methods in Applied Mathematics (2000). -# -# sum_impl requires length(a) > 1 -# -function sum_impl(f, a::AbstractArray, ifirst::Int, ilast::Int) - if ifirst + sum_pairwise_blocksize(f) >= ilast - sum_seq(f, a, ifirst, ilast) - else - imid = (ifirst + ilast) >>> 1 - sum_impl(f, a, ifirst, imid) + sum_impl(f, a, imid+1, ilast) - end -end -sum_impl{T<:Integer}(f::Union(IdFun,AbsFun,Abs2Fun), a::AbstractArray{T}, ifirst::Int, ilast::Int) = - sum_seq(f, a, ifirst, ilast) +# Note: sum_seq uses four accumulators, so each accumulator gets at most 256 numbers +sum_pairwise_blocksize(f) = 1024 -function sum(f::Union(Function,Func{1}), a::AbstractArray) - n = length(a) - n == 0 && error("Argument is empty.") - n == 1 && return addzero(evaluate(f, a[1])) - sum_impl(f, a, 1, n) -end - -for (fname, func, cutoff) in ((:sum, :IdFun, 16), (:sumabs, :AbsFun, 32), (:sumabs2, :Abs2Fun, 32)) - @eval function $fname{T}(a::AbstractArray{T}) - n = length(a) - n == 0 && return addzero(evaluate($func(), zero(T))) - n == 1 && return addzero(evaluate($func(), a[1])) - if n < $cutoff - # It is important that this is inlined to provide good - # performance for small inputs - @inbounds s = evaluate($func(), a[1]) + evaluate($func(), a[2]) - for i = 3:length(a) - @inbounds s += evaluate($func(), a[i]) - end - return s - end - sum_impl($func(), a, 1, n) - end -end +# This appears to show a benefit from a larger block size +sum_pairwise_blocksize(::Abs2Fun) = 4096 + +mapreduce_impl(f, op::AddFun, A::AbstractArray, ifirst::Int, ilast::Int) = + mapreduce_pairwise_impl(f, op, A, ifirst, ilast, sum_pairwise_blocksize(f)) + +sum(f::Union(Function,Func{1}), a) = mapreduce(f, AddFun(), a) +sum(a) = mapreduce(IdFun(), AddFun(), a) +sum(a::AbstractArray{Bool}) = countnz(a) +sumabs(a) = mapreduce(AbsFun(), AddFun(), a) +sumabs2(a) = mapreduce(Abs2Fun(), AddFun(), a) # Kahan (compensated) summation: O(1) error growth, at the expense # of a considerable increase in computational expense. @@ -360,8 +240,8 @@ if n == 0 return sumzero(T) end - s = addzero(A[1]) - c = sumzero(T) + c = zero(T) + s = A[1] + c for i in 2:n @inbounds Ai = A[i] t = s + Ai @@ -378,108 +258,15 @@ ## prod -prodtype{T}(::Type{T}) = typeof(zero(T) * zero(T)) -prodone{T}(::Type{T}) = one(T) * one(T) -multone(x) = x * one(x) - -function _prod(f, itr, s) - (x, s) = next(itr, s) - v = evaluate(f, x) - done(itr, s) && return multone(v) # multiplying by one for type stability - # specialize for length > 1 to have type-stable loop - (x, s) = next(itr, s) - result = v * evaluate(f, x) - while !done(itr, s) - (x, s) = next(itr, s) - result *= evaluate(f, x) - end - return result -end - -function prod(itr) - s = start(itr) - if done(itr, s) - if applicable(eltype, itr) - T = eltype(itr) - return prodone(T) - else - throw(ArgumentError("prod(itr) is undefined for empty collections; instead, do isempty(itr) ? o : prod(itr), where o is the correct type of identity for your product")) - end - end - _prod(IdFun(), itr, s) -end - -function prod(f::Function, itr) - s = start(itr) - done(itr, s) && error("Argument is empty.") - _prod(f, itr, s) -end - -prod(x::Number) = x +prod(f::Union(Function,Func{1}), a) = mapreduce(f, MulFun(), a) +prod(a) = mapreduce(IdFun(), MulFun(), a) prod(A::AbstractArray{Bool}) = error("use all() instead of prod() for boolean arrays") -function prod_impl{T}(f, A::AbstractArray{T}, first::Int, last::Int) - # pre-condition: last > first - i = first - @inbounds v = evaluate(f, A[i]) - @inbounds result = v * evaluate(f, A[i+=1]) - while i < last - @inbounds result *= evaluate(f, A[i+=1]) - end - return result -end - -function prod{T}(A::AbstractArray{T}) - n = length(A) - n == 0 && return prodone(T) - n == 1 && return multone(A[1]) - prod_impl(IdFun(), A, 1, n) -end - -function prod(f::Function, A::AbstractArray) - n = length(A) - n == 0 && error("Argument is empty.") - n == 1 && return multone(evaluate(f, A[1])) - prod_impl(f, A, 1, n) -end - - -## maximum & minimum - -function maximum(f::Union(Function,Func{1}), itr) - s = start(itr) - if done(itr, s) - error("argument is empty") - end - (x, s) = next(itr, s) - v = evaluate(f, x) - while !done(itr, s) - (x, s) = next(itr, s) - v = scalarmax(v, evaluate(f, x)) - end - return v -end - -function minimum(f::Union(Function,Func{1}), itr) - s = start(itr) - if done(itr, s) - error("argument is empty") - end - (x, s) = next(itr, s) - v = evaluate(f, x) - while !done(itr, s) - (x, s) = next(itr, s) - v = scalarmin(v, evaluate(f, x)) - end - return v -end +## maximum & minimum -maximum(itr) = maximum(IdFun(), itr) -minimum(itr) = minimum(IdFun(), itr) - -function maximum_impl(f, A::AbstractArray, first::Int, last::Int) +function mapreduce_impl(f, op::MaxFun, A::AbstractArray, first::Int, last::Int) # locate the first non NaN number v = evaluate(f, A[first]) i = first + 1 @@ -487,7 +274,6 @@ @inbounds v = evaluate(f, A[i]) i += 1 end - while i <= last @inbounds x = evaluate(f, A[i]) if x > v @@ -498,7 +284,7 @@ v end -function minimum_impl(f, A::AbstractArray, first::Int, last::Int) +function mapreduce_impl(f, op::MinFun, A::AbstractArray, first::Int, last::Int) # locate the first non NaN number v = evaluate(f, A[first]) i = first + 1 @@ -506,7 +292,6 @@ @inbounds v = evaluate(f, A[i]) i += 1 end - while i <= last @inbounds x = evaluate(f, A[i]) if x < v @@ -517,181 +302,144 @@ v end -function maximum(f::Union(Function,Func{1}), A::AbstractArray) - n = length(A) - n == 0 && error("Argument is empty.") - n == 1 && return evaluate(f, A[1]) - maximum_impl(f, A, 1, n) -end - -function minimum(f::Union(Function,Func{1}), A::AbstractArray) - n = length(A) - n == 0 && error("Argument is empty.") - n == 1 && return evaluate(f, A[1]) - minimum_impl(f, A, 1, n) -end +maximum(f::Union(Function,Func{1}), a) = mapreduce(f, MaxFun(), a) +minimum(f::Union(Function,Func{1}), a) = mapreduce(f, MinFun(), a) -maximum(A::AbstractArray) = maximum(IdFun(), A) -minimum(A::AbstractArray) = minimum(IdFun(), A) - -minabs(A::AbstractArray) = minimum(AbsFun(), A) - -# maxabs accepts empty array -function maxabs(A::AbstractArray) - n = length(A) - n == 0 && return abs(zero(T)) - n == 1 && return abs(A[1]) - maximum_impl(AbsFun(), A, 1, n) -end - - -maximum(x::Real) = x -minimum(x::Real) = x -maxabs(x::Number) = abs(x) -minabs(x::Number) = abs(x) +maximum(a) = mapreduce(IdFun(), MaxFun(), a) +minimum(a) = mapreduce(IdFun(), MinFun(), a) +maxabs(a) = mapreduce(AbsFun(), MaxFun(), a) +minabs(a) = mapreduce(AbsFun(), MinFun(), a) ## extrema extrema(r::Range) = (minimum(r), maximum(r)) +extrema(x::Real) = (x, x) function extrema(itr) s = start(itr) - if done(itr, s) - error("argument is empty") - end + done(itr, s) && error("argument is empty") (v, s) = next(itr, s) - vmin = v - vmax = v - while !done(itr, s) + while v != v && !done(itr, s) (x, s) = next(itr, s) - if x == x - if x > vmax - vmax = x - elseif x < vmin - vmin = x - end - end - end - return (vmin, vmax) -end - -function extrema{T<:Real}(A::AbstractArray{T}) - if isempty(A); error("argument must not be empty"); end - n = length(A) - - # locate the first non NaN number - v = A[1] - i = 2 - while v != v && i <= n - @inbounds v = A[i] - i += 1 + v = x end - vmin = v vmax = v - while i <= n - @inbounds v = A[i] - if v > vmax - vmax = v - elseif v < vmin - vmin = v + while !done(itr, s) + (x, s) = next(itr, s) + if x > vmax + vmax = x + elseif x < vmin + vmin = x end - i += 1 end - return (vmin, vmax) end - ## all & any -function all(itr) +function mapfoldl(f, ::AndFun, itr) for x in itr - if !x + if !evaluate(f, x) return false end end return true end -function any(itr) +function mapfoldl(f, ::OrFun, itr) for x in itr - if x + if evaluate(f, x) return true end end return false end -function any(pred::Union(Function,Func{1}), itr) - for x in itr - if pred(x) - return true +function mapreduce_impl(f, op::AndFun, A::AbstractArray, ifirst::Int, ilast::Int) + while ifirst <= ilast + @inbounds x = A[ifirst] + if !evaluate(f, x) + return false end + ifirst += 1 end - return false + return true end -function all(pred::Union(Function,Func{1}), itr) - for x in itr - if !pred(x) - return false +function mapreduce_impl(f, op::OrFun, A::AbstractArray, ifirst::Int, ilast::Int) + while ifirst <= ilast + @inbounds x = A[ifirst] + if evaluate(f, x) + return true end + ifirst += 1 end - return true + return false end +all(a) = mapreduce(IdFun(), AndFun(), a) +any(a) = mapreduce(IdFun(), OrFun(), a) -###### mapreduce ###### +all(pred::Union(Function,Func{1}), a) = mapreduce(pred, AndFun(), a) +any(pred::Union(Function,Func{1}), a) = mapreduce(pred, OrFun(), a) -function mapreduce(f::Callable, op::Callable, itr) - s = start(itr) - if done(itr, s) - error("argument is empty") - end - (x, s) = next(itr, s) - v = f(x) - if done(itr, s) - return v - else # specialize for length > 1 to have a hopefully type-stable loop - (x, s) = next(itr, s) - result = op(v, f(x)) - while !done(itr, s) - (x, s) = next(itr, s) - result = op(result, f(x)) + +## in & contains + +immutable EqX{T} <: Func{1} + x::T +end +EqX{T}(x::T) = EqX{T}(x) +evaluate(f::EqX, y) = (y == f.x) + +in(x, itr) = any(EqX(x), itr) + +const ∈ = in +∉(x, itr)=!∈(x, itr) +∋(itr, x)= ∈(x, itr) +∌(itr, x)=!∋(itr, x) + +function contains(itr, x) + depwarn("contains(collection, item) is deprecated, use in(item, collection) instead", :contains) + in(x, itr) +end + +function contains(eq::Function, itr, x) + for y in itr + if eq(y, x) + return true end - return result end + return false end -function mapreduce(f::Callable, op::Callable, v0, itr) - v = v0 + +## countnz & count + +function count(pred::Union(Function,Func{1}), itr) + n = 0 for x in itr - v = op(v,f(x)) + if evaluate(pred, x) + n += 1 + end end - return v + return n end -# pairwise reduction, requires n > 1 (to allow type-stable loop) -function mr_pairwise(f::Callable, op::Callable, A::AbstractArray, i1,n) - if n < 128 - @inbounds v = op(f(A[i1]), f(A[i1+1])) - for i = i1+2:i1+n-1 - @inbounds v = op(v,f(A[i])) +function count(pred::Union(Function,Func{1}), a::AbstractArray) + n = 0 + for i = 1:length(a) + @inbounds if evaluate(pred, a[i]) + n += 1 end - return v - else - n2 = div(n,2) - return op(mr_pairwise(f,op,A, i1,n2), mr_pairwise(f,op,A, i1+n2,n-n2)) end + return n end -function mapreduce(f::Callable, op::Callable, A::AbstractArray) - n = length(A) - n == 0 ? error("argument is empty") : n == 1 ? f(A[1]) : mr_pairwise(f,op,A, 1,n) -end -function mapreduce(f::Callable, op::Callable, v0, A::AbstractArray) - n = length(A) - n == 0 ? v0 : n == 1 ? op(v0, f(A[1])) : op(v0, mr_pairwise(f,op,A, 1,n)) -end + +type NotEqZero <: Func{1} end +evaluate(NotEqZero, x) = (x != 0) + +countnz(a) = count(NotEqZero(), a) diff -Nru julia-0.3.0/base/regex.jl julia-0.3.0/base/regex.jl --- julia-0.3.0/base/regex.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/regex.jl 2014-06-03 19:33:18.000000000 +0000 @@ -9,6 +9,8 @@ options::Uint32 regex::Ptr{Void} extra::Ptr{Void} + ovec::Vector{Int32} + function Regex(pattern::String, options::Integer) pattern = bytestring(pattern) @@ -16,7 +18,7 @@ if (options & ~PCRE.OPTIONS_MASK) != 0 error("invalid regex options: $options") end - re = compile(new(pattern, options, C_NULL, C_NULL)) + re = compile(new(pattern, options, C_NULL, C_NULL, Array(Int32, 0))) finalizer(re, function(re::Regex) re.extra != C_NULL && PCRE.free_study(re.extra) @@ -43,6 +45,9 @@ if regex.regex == C_NULL regex.regex = PCRE.compile(regex.pattern, regex.options & PCRE.COMPILE_MASK) regex.extra = PCRE.study(regex.regex, PCRE.STUDY_JIT_COMPILE) + ncap = PCRE.info(regex.regex, regex.extra, + PCRE.INFO_CAPTURECOUNT, Int32) + resize!(regex.ovec, 3(ncap+1)) end regex end @@ -98,23 +103,28 @@ function ismatch(r::Regex, s::String, offset::Integer=0) compile(r) - PCRE.exec(r.regex, r.extra, bytestring(s), offset, r.options & PCRE.EXECUTE_MASK, false) + return PCRE.exec(r.regex, r.extra, bytestring(s), offset, r.options & PCRE.EXECUTE_MASK, + r.ovec) end + function ismatch(r::Regex, s::SubString, offset::Integer=0) compile(r) - PCRE.exec(r.regex, r.extra, s, offset, r.options & PCRE.EXECUTE_MASK, false) + return PCRE.exec(r.regex, r.extra, s, offset, r.options & PCRE.EXECUTE_MASK, + r.ovec) end function match(re::Regex, str::UTF8String, idx::Integer, add_opts::Uint32=uint32(0)) opts = re.options & PCRE.EXECUTE_MASK | add_opts compile(re) - m, n = PCRE.exec(re.regex, re.extra, str, idx-1, opts, true) - if isempty(m); return nothing; end - mat = SubString(str, m[1]+1, m[2]) + if !PCRE.exec(re.regex, re.extra, str, idx-1, opts, re.ovec) + return nothing + end + n = length(re.ovec)/3 - 1 + mat = SubString(str, re.ovec[1]+1, re.ovec[2]) cap = Union(Nothing,SubString{UTF8String})[ - m[2i+1] < 0 ? nothing : SubString(str, m[2i+1]+1, m[2i+2]) for i=1:n ] - off = Int[ m[2i+1]::Int32+1 for i=1:n ] - RegexMatch(mat, cap, m[1]+1, off) + re.ovec[2i+1] < 0 ? nothing : SubString(str, re.ovec[2i+1]+1, re.ovec[2i+2]) for i=1:n ] + off = Int[ re.ovec[2i+1]::Int32+1 for i=1:n ] + RegexMatch(mat, cap, re.ovec[1]+1, off) end match(re::Regex, str::Union(ByteString,SubString), idx::Integer, add_opts::Uint32=uint32(0)) = @@ -173,8 +183,8 @@ end opts = re.options & PCRE.EXECUTE_MASK compile(re) - m, n = PCRE.exec(re.regex, re.extra, str, idx-1, opts, true) - isempty(m) ? (0:-1) : ((m[1]+1):prevind(str,m[2]+1)) + PCRE.exec(re.regex, re.extra, str, idx-1, opts, re.ovec) ? + ((re.ovec[1]+1):prevind(str,re.ovec[2]+1)) : (0:-1) end search(s::String, r::Regex, idx::Integer) = error("regex search is only available for bytestrings; use bytestring(s) to convert") diff -Nru julia-0.3.0/base/socket.jl julia-0.3.0/base/socket.jl --- julia-0.3.0/base/socket.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/socket.jl 2014-06-03 19:33:18.000000000 +0000 @@ -310,7 +310,7 @@ uv.handle = 0 end -function uvfinalize(uv::Union(TTY,Pipe,TcpServer,TcpSocket)) +function uvfinalize(uv::Union(TTY,Pipe,PipeServer,TcpServer,TcpSocket)) if (uv.status != StatusUninit && uv.status != StatusInit) close(uv) end diff -Nru julia-0.3.0/base/sort.jl julia-0.3.0/base/sort.jl --- julia-0.3.0/base/sort.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/base/sort.jl 2014-06-03 19:33:18.000000000 +0000 @@ -125,6 +125,7 @@ # index of the first value of vector a that is greater than or equal to x; # returns length(v)+1 if x is greater than all values in v. function searchsortedfirst(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering) + 1 <= lo <= hi <= length(v) || throw(BoundsError()) lo = lo-1 hi = hi+1 @inbounds while lo < hi-1 @@ -141,6 +142,7 @@ # index of the last value of vector a that is less than or equal to x; # returns 0 if x is less than all values of v. function searchsortedlast(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering) + 1 <= lo <= hi <= length(v) || throw(BoundsError()) lo = lo-1 hi = hi+1 @inbounds while lo < hi-1 @@ -158,6 +160,7 @@ # if v does not contain x, returns a 0-length range # indicating the insertion point of x function searchsorted(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering) + 1 <= lo <= hi <= length(v) || throw(BoundsError()) lo = lo-1 hi = hi+1 @inbounds while lo < hi-1 @@ -176,37 +179,21 @@ end function searchsortedlast{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering) - if step(a) == 0 - lt(o, x, first(a)) ? 0 : length(a) - else - n = max(min(iround((x-first(a))/step(a))+1,length(a)),1) - lt(o, x, a[n]) ? n-1 : n - end + n = max(min(iround((x-first(a))/step(a))+1,length(a)),1) + lt(o, x, a[n]) ? n-1 : n end function searchsortedfirst{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering) - if step(a) == 0 - lt(o, first(a), x) ? length(a)+1 : 1 - else - n = max(min(iround((x-first(a))/step(a))+1,length(a)),1) - lt(o, a[n] ,x) ? n+1 : n - end + n = max(min(iround((x-first(a))/step(a))+1,length(a)),1) + lt(o, a[n] ,x) ? n+1 : n end function searchsortedlast{T<:Integer}(a::Range{T}, x::Real, o::DirectOrdering) - if step(a) == 0 - lt(o, x, first(a)) ? 0 : length(a) - else - max(min(fld(ifloor(x)-first(a),step(a))+1,length(a)),0) - end + max(min(fld(ifloor(x)-first(a),step(a))+1,length(a)),0) end function searchsortedfirst{T<:Integer}(a::Range{T}, x::Real, o::DirectOrdering) - if step(a) == 0 - lt(o, first(a), x) ? length(a)+1 : 1 - else - max(min(-fld(ifloor(-x)+first(a),step(a))+1,length(a)+1),1) - end + max(min(-fld(ifloor(-x)+first(a),step(a))+1,length(a)+1),1) end searchsorted{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering) = diff -Nru julia-0.3.0/base/special/bessel.jl julia-0.3.0/base/special/bessel.jl --- julia-0.3.0/base/special/bessel.jl 1970-01-01 00:00:00.000000000 +0000 +++ julia-0.3.0/base/special/bessel.jl 2014-06-03 19:33:18.000000000 +0000 @@ -0,0 +1,280 @@ +for jy in ("j","y"), nu in (0,1) + jynu = Expr(:quote, symbol(string(jy,nu))) + jynuf = Expr(:quote, symbol(string(jy,nu,"f"))) + bjynu = symbol(string("bessel",jy,nu)) + if jy == "y" + @eval begin + $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm), Float64, (Float64,), x), x) + $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x) + end + else + @eval begin + $bjynu(x::Float64) = ccall(($jynu,libm), Float64, (Float64,), x) + $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x) + end + end + @eval begin + $bjynu(x::Real) = $bjynu(float(x)) + $bjynu(x::Complex) = $(symbol(string("bessel",jy)))($nu,x) + @vectorize_1arg Number $bjynu + end +end + +type AmosException <: Exception + info::Int32 +end + +let + const ai::Array{Float64,1} = Array(Float64,2) + const ae::Array{Int32,1} = Array(Int32,2) +global airy +function airy(k::Int, z::Complex128) + id = int32(k==1 || k==3) + if k == 0 || k == 1 + ccall((:zairy_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + &real(z), &imag(z), + &id, &1, + pointer(ai,1), pointer(ai,2), + pointer(ae,1), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(ai[1],ai[2]) + else + throw(AmosException(ae[2])) + end + elseif k == 2 || k == 3 + ccall((:zbiry_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), + &real(z), &imag(z), + &id, &1, + pointer(ai,1), pointer(ai,2), + pointer(ae,1)) + if ae[1] == 0 || ae[1] == 3 # ignore underflow + return complex(ai[1],ai[2]) + else + throw(AmosException(ae[2])) + end + else + error("invalid argument") + end +end +end + +airy(z) = airy(0,z) +@vectorize_1arg Number airy +airyprime(z) = airy(1,z) +@vectorize_1arg Number airyprime +airyai(z) = airy(0,z) +@vectorize_1arg Number airyai +airyaiprime(z) = airy(1,z) +@vectorize_1arg Number airyaiprime +airybi(z) = airy(2,z) +@vectorize_1arg Number airybi +airybiprime(z) = airy(3,z) +@vectorize_1arg Number airybiprime + +airy(k::Number, x::FloatingPoint) = oftype(x, real(airy(k, complex(x)))) +airy(k::Number, x::Real) = airy(k, float(x)) +airy(k::Number, z::Complex64) = complex64(airy(k, complex128(z))) +airy(k::Number, z::Complex) = airy(convert(Int,k), complex128(z)) +@vectorize_2arg Number airy + +const cy = Array(Float64,2) +const ae = Array(Int32,2) +const wrk = Array(Float64,2) + +function _besselh(nu::Float64, k::Integer, z::Complex128) + ccall((:zbesh_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + &real(z), &imag(z), &nu, &1, &k, &1, + pointer(cy,1), pointer(cy,2), + pointer(ae,1), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end +end + +function _besseli(nu::Float64, z::Complex128) + ccall((:zbesi_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + &real(z), &imag(z), &nu, &1, &1, + pointer(cy,1), pointer(cy,2), + pointer(ae,1), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end +end + +function _besselj(nu::Float64, z::Complex128) + ccall((:zbesj_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + &real(z), &imag(z), &nu, &1, &1, + pointer(cy,1), pointer(cy,2), + pointer(ae,1), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end +end + +function _besselk(nu::Float64, z::Complex128) + ccall((:zbesk_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + &real(z), &imag(z), &nu, &1, &1, + pointer(cy,1), pointer(cy,2), + pointer(ae,1), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end +end + +function _bessely(nu::Float64, z::Complex128) + ccall((:zbesy_,openspecfun), Void, + (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, + Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), + &real(z), &imag(z), &nu, &1, &1, + pointer(cy,1), pointer(cy,2), + pointer(ae,1), pointer(wrk,1), + pointer(wrk,2), pointer(ae,2)) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end +end + +function besselh(nu::Float64, k::Integer, z::Complex128) + if nu < 0 + s = (k == 1) ? 1 : -1 + return _besselh(-nu, k, z) * complex(cospi(nu),-s*sinpi(nu)) + end + return _besselh(nu, k, z) +end + +function besseli(nu::Float64, z::Complex128) + if nu < 0 + return _besseli(-nu,z) - 2_besselk(-nu,z)*sinpi(nu)/pi + else + return _besseli(nu, z) + end +end + +function besselj(nu::Float64, z::Complex128) + if nu < 0 + return _besselj(-nu,z)cos(pi*nu) + _bessely(-nu,z)*sinpi(nu) + else + return _besselj(nu, z) + end +end + +besselj(nu::Integer, x::FloatingPoint) = typemin(Int32) <= nu <= typemax(Int32) ? + oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x)) : + besselj(float64(nu), x) + +besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ? + ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) : + besselj(float64(nu), x) + +besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z) + +function bessely(nu::Float64, z::Complex128) + if nu < 0 + return _bessely(-nu,z)*cospi(nu) - _besselj(-nu,z)*sinpi(nu) + else + return _bessely(nu, z) + end +end + +besselh(nu, z) = besselh(nu, 1, z) +besselh(nu::Real, k::Integer, z::Complex64) = complex64(besselh(float64(nu), k, complex128(z))) +besselh(nu::Real, k::Integer, z::Complex) = besselh(float64(nu), k, complex128(z)) +besselh(nu::Real, k::Integer, x::Real) = besselh(float64(nu), k, complex128(x)) +@vectorize_2arg Number besselh + +besseli(nu::Real, z::Complex64) = complex64(besseli(float64(nu), complex128(z))) +besseli(nu::Real, z::Complex) = besseli(float64(nu), complex128(z)) +besseli(nu::Real, x::Integer) = besseli(nu, float64(x)) +function besseli(nu::Real, x::FloatingPoint) + if x < 0 && !isinteger(nu) + throw(DomainError()) + end + oftype(x, real(besseli(float64(nu), complex128(x)))) +end +@vectorize_2arg Number besseli + +function besselj(nu::FloatingPoint, x::FloatingPoint) + if isinteger(nu) + if typemin(Int32) <= nu <= typemax(Int32) + return besselj(int(nu), x) + end + elseif x < 0 + throw(DomainError()) + end + oftype(x, real(besselj(float64(nu), complex128(x)))) +end + +besselj(nu::Real, z::Complex64) = complex64(besselj(float64(nu), complex128(z))) +besselj(nu::Real, z::Complex) = besselj(float64(nu), complex128(z)) +besselj(nu::Real, x::Integer) = besselj(nu, float(x)) +@vectorize_2arg Number besselj + +besselk(nu::Real, z::Complex64) = complex64(besselk(float64(nu), complex128(z))) +besselk(nu::Real, z::Complex) = besselk(float64(nu), complex128(z)) +besselk(nu::Real, x::Integer) = besselk(nu, float64(x)) +function besselk(nu::Real, x::FloatingPoint) + if x < 0 + throw(DomainError()) + end + if x == 0 + return oftype(x, Inf) + end + oftype(x, real(besselk(float64(nu), complex128(x)))) +end +@vectorize_2arg Number besselk + +bessely(nu::Real, z::Complex64) = complex64(bessely(float64(nu), complex128(z))) +bessely(nu::Real, z::Complex) = bessely(float64(nu), complex128(z)) +bessely(nu::Real, x::Integer) = bessely(nu, float64(x)) +function bessely(nu::Real, x::FloatingPoint) + if x < 0 + throw(DomainError()) + end + if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32) + return bessely(int(nu), x) + end + oftype(x, real(bessely(float64(nu), complex128(x)))) +end +function bessely(nu::Integer, x::FloatingPoint) + if x < 0 + throw(DomainError()) + end + return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x)) +end +function bessely(nu::Integer, x::Float32) + if x < 0 + throw(DomainError()) + end + return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x) +end +@vectorize_2arg Number bessely + +hankelh1(nu, z) = besselh(nu, 1, z) +@vectorize_2arg Number hankelh1 + +hankelh2(nu, z) = besselh(nu, 2, z) +@vectorize_2arg Number hankelh2 + diff -Nru julia-0.3.0/base/special/erf.jl julia-0.3.0/base/special/erf.jl --- julia-0.3.0/base/special/erf.jl 1970-01-01 00:00:00.000000000 +0000 +++ julia-0.3.0/base/special/erf.jl 2014-06-03 19:33:18.000000000 +0000 @@ -0,0 +1,222 @@ +@unix_only if WORD_SIZE == 64 +# TODO: complex return only on 64-bit for now +for f in (:erf, :erfc, :erfcx, :erfi, :Dawson) + fname = (f === :Dawson) ? :dawson : f + @eval begin + ($fname)(z::Complex128) = complex128(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) + ($fname)(z::Complex64) = complex64(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), complex128(z), float64(eps(Float32)))) + ($fname)(z::Complex) = ($fname)(complex128(z)) + end +end +end +for f in (:erfcx, :erfi, :Dawson) + fname = (f === :Dawson) ? :dawson : f + @eval begin + ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), x) + ($fname)(x::Float32) = float32(ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), float64(x))) + ($fname)(x::Integer) = ($fname)(float(x)) + @vectorize_1arg Number $fname + end +end + +# Compute the inverse of the error function: erf(erfinv(x)) == x, +# using the rational approximants tabulated in: +# J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev +# approximations for the inverse of the error function," Math. Comp. 30, +# pp. 827--830 (1976). +# http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7 +# http://www.jstor.org/stable/2005402 +function erfinv(x::Float64) + a = abs(x) + if a >= 1.0 + if x == 1.0 + return inf(Float64) + elseif x == -1.0 + return -inf(Float64) + end + throw(DomainError()) + elseif a <= 0.75 # Table 17 in Blair et al. + t = x*x - 0.5625 + return x * @horner(t, 0.16030_49558_44066_229311e2, + -0.90784_95926_29603_26650e2, + 0.18644_91486_16209_87391e3, + -0.16900_14273_46423_82420e3, + 0.65454_66284_79448_7048e2, + -0.86421_30115_87247_794e1, + 0.17605_87821_39059_0) / + @horner(t, 0.14780_64707_15138_316110e2, + -0.91374_16702_42603_13936e2, + 0.21015_79048_62053_17714e3, + -0.22210_25412_18551_32366e3, + 0.10760_45391_60551_23830e3, + -0.20601_07303_28265_443e2, + 0.1e1) + elseif a <= 0.9375 # Table 37 in Blair et al. + t = x*x - 0.87890625 + return x * @horner(t, -0.15238_92634_40726_128e-1, + 0.34445_56924_13612_5216, + -0.29344_39867_25424_78687e1, + 0.11763_50570_52178_27302e2, + -0.22655_29282_31011_04193e2, + 0.19121_33439_65803_30163e2, + -0.54789_27619_59831_8769e1, + 0.23751_66890_24448) / + @horner(t, -0.10846_51696_02059_954e-1, + 0.26106_28885_84307_8511, + -0.24068_31810_43937_57995e1, + 0.10695_12997_33870_14469e2, + -0.23716_71552_15965_81025e2, + 0.24640_15894_39172_84883e2, + -0.10014_37634_97830_70835e2, + 0.1e1) + else # Table 57 in Blair et al. + t = 1.0 / sqrt(-log(1.0 - a)) + return @horner(t, 0.10501_31152_37334_38116e-3, + 0.10532_61131_42333_38164_25e-1, + 0.26987_80273_62432_83544_516, + 0.23268_69578_89196_90806_414e1, + 0.71678_54794_91079_96810_001e1, + 0.85475_61182_21678_27825_185e1, + 0.68738_08807_35438_39802_913e1, + 0.36270_02483_09587_08930_02e1, + 0.88606_27392_96515_46814_9) / + (copysign(t, x) * + @horner(t, 0.10501_26668_70303_37690e-3, + 0.10532_86230_09333_27531_11e-1, + 0.27019_86237_37515_54845_553, + 0.23501_43639_79702_53259_123e1, + 0.76078_02878_58012_77064_351e1, + 0.11181_58610_40569_07827_3451e2, + 0.11948_78791_84353_96667_8438e2, + 0.81922_40974_72699_07893_913e1, + 0.40993_87907_63680_15361_45e1, + 0.1e1)) + end +end + +function erfinv(x::Float32) + a = abs(x) + if a >= 1.0f0 + if x == 1.0f0 + return inf(Float32) + elseif x == -1.0f0 + return -inf(Float32) + end + throw(DomainError()) + elseif a <= 0.75f0 # Table 10 in Blair et al. + t = x*x - 0.5625f0 + return x * @horner(t, -0.13095_99674_22f2, + 0.26785_22576_0f2, + -0.92890_57365f1) / + @horner(t, -0.12074_94262_97f2, + 0.30960_61452_9f2, + -0.17149_97799_1f2, + 0.1f1) + elseif a <= 0.9375f0 # Table 29 in Blair et al. + t = x*x - 0.87890625f0 + return x * @horner(t, -0.12402_56522_1f0, + 0.10688_05957_4f1, + -0.19594_55607_8f1, + 0.42305_81357f0) / + @horner(t, -0.88276_97997f-1, + 0.89007_43359f0, + -0.21757_03119_6f1, + 0.1f1) + else # Table 50 in Blair et al. + t = 1.0f0 / sqrt(-log(1.0f0 - a)) + return @horner(t, 0.15504_70003_116f0, + 0.13827_19649_631f1, + 0.69096_93488_87f0, + -0.11280_81391_617f1, + 0.68054_42468_25f0, + -0.16444_15679_1f0) / + (copysign(t, x) * + @horner(t, 0.15502_48498_22f0, + 0.13852_28141_995f1, + 0.1f1)) + end +end + +erfinv(x::Integer) = erfinv(float(x)) +@vectorize_1arg Real erfinv + +# Inverse complementary error function: use Blair tables for y = 1-x, +# exploiting the greater accuracy of y (vs. x) when y is small. +function erfcinv(y::Float64) + if y > 0.0625 + return erfinv(1.0 - y) + elseif y <= 0.0 + if y == 0.0 + return inf(Float64) + end + throw(DomainError()) + elseif y >= 1e-100 # Table 57 in Blair et al. + t = 1.0 / sqrt(-log(y)) + return @horner(t, 0.10501_31152_37334_38116e-3, + 0.10532_61131_42333_38164_25e-1, + 0.26987_80273_62432_83544_516, + 0.23268_69578_89196_90806_414e1, + 0.71678_54794_91079_96810_001e1, + 0.85475_61182_21678_27825_185e1, + 0.68738_08807_35438_39802_913e1, + 0.36270_02483_09587_08930_02e1, + 0.88606_27392_96515_46814_9) / + (t * + @horner(t, 0.10501_26668_70303_37690e-3, + 0.10532_86230_09333_27531_11e-1, + 0.27019_86237_37515_54845_553, + 0.23501_43639_79702_53259_123e1, + 0.76078_02878_58012_77064_351e1, + 0.11181_58610_40569_07827_3451e2, + 0.11948_78791_84353_96667_8438e2, + 0.81922_40974_72699_07893_913e1, + 0.40993_87907_63680_15361_45e1, + 0.1e1)) + else # Table 80 in Blair et al. + t = 1.0 / sqrt(-log(y)) + return @horner(t, 0.34654_29858_80863_50177e-9, + 0.25084_67920_24075_70520_55e-6, + 0.47378_13196_37286_02986_534e-4, + 0.31312_60375_97786_96408_3388e-2, + 0.77948_76454_41435_36994_854e-1, + 0.70045_68123_35816_43868_271e0, + 0.18710_42034_21679_31668_683e1, + 0.71452_54774_31351_45428_3e0) / + (t * @horner(t, 0.34654_29567_31595_11156e-9, + 0.25084_69079_75880_27114_87e-6, + 0.47379_53129_59749_13536_339e-4, + 0.31320_63536_46177_68848_0813e-2, + 0.78073_48906_27648_97214_733e-1, + 0.70715_04479_95337_58619_993e0, + 0.19998_51543_49112_15105_214e1, + 0.15072_90269_27316_80008_56e1, + 0.1e1)) + end +end + +function erfcinv(y::Float32) + if y > 0.0625f0 + return erfinv(1.0f0 - y) + elseif y <= 0.0f0 + if y == 0.0f0 + return inf(Float32) + end + throw(DomainError()) + else # Table 50 in Blair et al. + t = 1.0f0 / sqrt(-log(y)) + return @horner(t, 0.15504_70003_116f0, + 0.13827_19649_631f1, + 0.69096_93488_87f0, + -0.11280_81391_617f1, + 0.68054_42468_25f0, + -0.16444_15679_1f0) / + (t * + @horner(t, 0.15502_48498_22f0, + 0.13852_28141_995f1, + 0.1f1)) + end +end + +erfcinv(x::Integer) = erfcinv(float(x)) +@vectorize_1arg Real erfcinv + diff -Nru julia-0.3.0/base/special/gamma.jl julia-0.3.0/base/special/gamma.jl --- julia-0.3.0/base/special/gamma.jl 1970-01-01 00:00:00.000000000 +0000 +++ julia-0.3.0/base/special/gamma.jl 2014-06-03 19:33:18.000000000 +0000 @@ -0,0 +1,533 @@ +gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm), Float64, (Float64,), x), x) +gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm), Float32, (Float32,), x), x) +gamma(x::Real) = gamma(float(x)) +@vectorize_1arg Number gamma + +function lgamma_r(x::Float64) + signp = Array(Int32, 1) + y = ccall((:lgamma_r,libm), Float64, (Float64, Ptr{Int32}), x, signp) + return y, signp[1] +end +function lgamma_r(x::Float32) + signp = Array(Int32, 1) + y = ccall((:lgammaf_r,libm), Float32, (Float32, Ptr{Int32}), x, signp) + return y, signp[1] +end +lgamma_r(x::Real) = lgamma_r(float(x)) + +lfact(x::Real) = (x<=1 ? zero(float(x)) : lgamma(x+one(x))) +@vectorize_1arg Number lfact + +const clg_coeff = [76.18009172947146, + -86.50532032941677, + 24.01409824083091, + -1.231739572450155, + 0.1208650973866179e-2, + -0.5395239384953e-5] + +function clgamma_lanczos(z) + const sqrt2pi = 2.5066282746310005 + + y = x = z + temp = x + 5.5 + zz = log(temp) + zz = zz * (x+0.5) + temp -= zz + ser = complex(1.000000000190015, 0) + for j=1:6 + y += 1.0 + zz = clg_coeff[j]/y + ser += zz + end + zz = sqrt2pi*ser / x + return log(zz) - temp +end + +function lgamma(z::Complex) + if real(z) <= 0.5 + a = clgamma_lanczos(1-z) + b = log(sinpi(z)) + const logpi = 1.14472988584940017 + z = logpi - b - a + else + z = clgamma_lanczos(z) + end + complex(real(z), angle_restrict_symm(imag(z))) +end + +gamma(z::Complex) = exp(lgamma(z)) + +# Derivatives of the digamma function +function psifn(x::Float64, n::Int, kode::Int, m::Int) +# Translated from http://www.netlib.org/slatec/src/dpsifn.f +# Note: Underflow handling at 380 in original is skipped + const nmax = 100 + ans = Array(Float64, m) +#----------------------------------------------------------------------- +# bernoulli numbers +#----------------------------------------------------------------------- + const b = [1.00000000000000000e+00, + -5.00000000000000000e-01,1.66666666666666667e-01, + -3.33333333333333333e-02,2.38095238095238095e-02, + -3.33333333333333333e-02,7.57575757575757576e-02, + -2.53113553113553114e-01,1.16666666666666667e+00, + -7.09215686274509804e+00,5.49711779448621554e+01, + -5.29124242424242424e+02,6.19212318840579710e+03, + -8.65802531135531136e+04,1.42551716666666667e+06, + -2.72982310678160920e+07,6.01580873900642368e+08, + -1.51163157670921569e+10,4.29614643061166667e+11, + -1.37116552050883328e+13,4.88332318973593167e+14, + -1.92965793419400681e+16] + trm = Array(Float64, 22) + trmr = Array(Float64, 100) +#***first executable statement dpsifn + if x <= 0.0 throw(DomainError()) end + if n < 0 error("n must be non-negative") end + if kode < 1 | kode > 2 error("kode must be one or two") end + if m < 1 error("m must be larger than one") end + mm = m + const nx = min(-exponent(realmin(Float64)) + 1, exponent(realmax(Float64))) + const r1m5 = log10(2) + const r1m4 = Base.eps(Float64) * 0.5 + const wdtol = max(r1m4, 0.5e-18) +#----------------------------------------------------------------------- +# elim = approximate exponential over and underflow limit +#----------------------------------------------------------------------- + const elim = 2.302*(nx*r1m5 - 3.0) + xln = log(x) + nn = n + mm - 1 + fn = nn + t = (fn + 1)*xln +#----------------------------------------------------------------------- +# overflow and underflow test for small and large x +#----------------------------------------------------------------------- + if abs(t) > elim + if t <= 0.0 error("n too large") end + error("overflow; x too small or n+m-1 too large or both") + end + if x < wdtol + ans[1] = x^(-n - 1) + if mm != 1 + k = 1 + for i = 2:mm + ans[k + 1] = ans[k]/x + k += 1 + end + end + if n != 0 return ans end + if kode == 2 ans[1] = ans[1] + xln end + return ans + end +#----------------------------------------------------------------------- +# compute xmin and the number of terms of the series, fln+1 +#----------------------------------------------------------------------- + rln = r1m5 * precision(x) + rln = min(rln, 18.06) + fln = max(rln, 3.0) - 3.0 + yint = 3.50 + 0.40*fln + slope = 0.21 + fln*(0.0006038*fln + 0.008677) + xm = yint + slope*fn + mx = itrunc(xm) + 1 + xmin = mx + if n != 0 + xm = -2.302*rln - min(0.0,xln) + arg = xm/n + arg = min(0.0,arg) + eps = exp(arg) + xm = 1.0 - eps + if abs(arg) < 1.0e-3 xm = -arg end + fln = x*xm/eps + xm = xmin - x + if (xm > 7.0) & (fln < 15.0) + nn = itrunc(fln) + 1 + np = n + 1 + t1 = (n + 1)*xln + t = exp(-t1) + s = t + den = x + for i = 1:nn + den += 1.0 + trm[i] = den^(-np) + s += trm[i] + end + ans[1] = s + if n == 0 + if kode == 2 ans[1] = s + xln end + end + if mm == 1 return ans end +#----------------------------------------------------------------------- +# generate higher derivatives, j.gt.n +#----------------------------------------------------------------------- + tol = wdtol/5.0 + for j = 2:mm + t = t/x + s = t + tols = t*tol + den = x + for i = 1:nn + den += 1.0 + trm[i] = trm[i]/den + s += trm[i] + if trm[i] < tols break end + end + ans[j] = s + end + return ans + end + end + + xdmy = x + xdmln = xln + xinc = 0.0 + if x < xmin + nx = itrunc(x) + xinc = xmin - nx + xdmy = x + xinc + xdmln = log(xdmy) + end +#----------------------------------------------------------------------- +# generate w(n+mm-1,x) by the asymptotic expansion +#----------------------------------------------------------------------- + t = fn*xdmln + t1 = xdmln + xdmln + t2 = t + xdmln + tk = max(abs(t), abs(t1), abs(t2)) + if tk > elim error("underflow") end + tss = exp(-t) + tt = 0.5/xdmy + t1 = tt + tst = wdtol*tt + if nn != 0 t1 = tt + 1.0/fn end + rxsq = 1.0/(xdmy*xdmy) + ta = 0.5*rxsq + t = (fn + 1)*ta + s = t*b[3] + if abs(s) >= tst + tk = 2.0 + for k = 4:22 + t = t*((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0))*rxsq + trm[k] = t*b[k] + if abs(trm[k]) < tst break end + s += trm[k] + tk += 2.0 + end + end + s = (s + t1)*tss + while true + if xinc != 0.0 +#----------------------------------------------------------------------- +# backward recur from xdmy to x +#----------------------------------------------------------------------- + nx = itrunc(xinc) + np = nn + 1 + if nx > nmax error("n too large") end + if nn == 0 break end + xm = xinc - 1.0 + fx = x + xm +#----------------------------------------------------------------------- +# this loop should not be changed. fx is accurate when x is small +#----------------------------------------------------------------------- + for i = 1:nx + trmr[i] = fx^(-np) + s += trmr[i] + xm -= 1.0 + fx = x + xm + end + end + ans[mm] = s + if fn == 0 + if kode != 2 + ans[1] = s - xdmln + return ans + end + if xdmy == x return ans end + xq = xdmy/x + ans[1] = s - log(xq) + return ans + end +#----------------------------------------------------------------------- +# generate lower derivatives, j.lt.n+mm-1 +#----------------------------------------------------------------------- + if mm == 1 return ans end + for j = 2:mm + fn -= 1 + tss *= xdmy + t1 = tt + if fn != 0 t1 = tt + 1.0/fn end + t = (fn + 1)*ta + s = t*b[3] + if abs(s) >= tst + tk = 4 + fn + for k = 4:22 #110 + trm[k] = trm[k]*(fn + 1)/tk + if abs(trm[k]) < tst break end + s += trm[k] + tk += 2.0 + end + end + s = (s + t1)*tss + if xinc != 0.0 + if fn == 0 break end + xm = xinc - 1.0 + fx = x + xm + for i = 1:nx + trmr[i] = trmr[i]*fx + s += trmr[i] + xm -= 1.0 + fx = x + xm + end + end + mx = mm - j + 1 + ans[mx] = s + if fn == 0 + if kode != 2 + ans[1] = s - xdmln + return ans + end + if xdmy == x return ans end + xq = xdmy/x + ans[1] = s - log(xq) + return ans + end + end + if fn == 0 break end + return ans + end +#----------------------------------------------------------------------- +# recursion for n = 0 +#----------------------------------------------------------------------- + for i = 1:nx + s += 1.0/(x + nx - i) + end + if kode != 2 + ans[1] = s - xdmln + return ans + end + if xdmy == x return ans end + xq = xdmy/x + ans[1] = s - log(xq) + return ans +end +polygamma(k::Int, x::Float64) = (2rem(k,2) - 1)*psifn(x, k, 1, 1)[1]/gamma(k + 1) +polygamma(k::Int, x::Float32) = float32(polygamma(k, float64(x))) +polygamma(k::Int, x::Real) = polygamma(k, float64(x)) + +# Translation of psi.c from cephes +function digamma(x::Float64) + negative = false + nz = 0.0 + + if x <= 0.0 + negative = true + q = x + p = floor(q) + if p == q + return NaN + end + + nz = q - p + if nz != 0.5 + if nz > 0.5 + p += 1.0 + nz = q - p + end + nz = pi / tan(pi * nz) + else + nz = 0.0 + end + x = 1.0 - x + end + + if x <= 10.0 && x == floor(x) + y = 0.0 + for i = 1:x-1 + y += 1.0 / i + end + y -= γ # γ == -digamma(1) == 0.5772156649015328606065121; + + if negative + y -= nz + end + return y + end + + w = 0.0 + while x < 10.0 + w += 1.0 / x + x += 1.0 + end + + if x < 1.0e17 + z = 1.0 / (x*x) + y = @horner(z, 8.33333333333333333333e-2, -8.33333333333333333333e-3, 3.96825396825396825397e-3, + -4.16666666666666666667e-3, 7.57575757575757575758e-3,-2.10927960927960927961e-2, + 8.33333333333333333333e-2) + y *= z + else + y = 0.0 + end + + y = log(x) - 0.5/x - y - w + + if negative + y -= nz + end + + return y +end +digamma(x::Float32) = float32(digamma(float64(x))) +digamma(x::Real) = digamma(float64(x)) +@vectorize_1arg Real digamma + +trigamma(x::Real) = polygamma(1, x) +@vectorize_1arg Real trigamma + +# Inverse digamma function +# +# Implementation of fixed point algorithm described in +# "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000 +function invdigamma(y::Float64) + # Closed form initial estimates + if y >= -2.22 + x_old = exp(y) + 0.5 + x_new = x_old + else + x_old = -1.0 / (y - digamma(1.0)) + x_new = x_old + end + + # Fixed point algorithm + delta = Inf + iteration = 0 + while delta > 1e-12 && iteration < 25 + iteration += 1 + x_new = x_old - (digamma(x_old) - y) / trigamma(x_old) + delta = abs(x_new - x_old) + x_old = x_new + end + + return x_new +end +invdigamma(x::Float32) = float32(invdigamma(float64(x))) +invdigamma(x::Real) = invdigamma(float64(x)) +@vectorize_1arg Real invdigamma + +function beta(x::Number, w::Number) + yx, sx = lgamma_r(x) + yw, sw = lgamma_r(w) + yxw, sxw = lgamma_r(x+w) + return copysign(exp(yx + yw - yxw), sx*sw*sxw) +end +lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w) +@vectorize_2arg Number beta +@vectorize_2arg Number lbeta + +const eta_coeffs = + [.99999999999999999997, + -.99999999999999999821, + .99999999999999994183, + -.99999999999999875788, + .99999999999998040668, + -.99999999999975652196, + .99999999999751767484, + -.99999999997864739190, + .99999999984183784058, + -.99999999897537734890, + .99999999412319859549, + -.99999996986230482845, + .99999986068828287678, + -.99999941559419338151, + .99999776238757525623, + -.99999214148507363026, + .99997457616475604912, + -.99992394671207596228, + .99978893483826239739, + -.99945495809777621055, + .99868681159465798081, + -.99704078337369034566, + .99374872693175507536, + -.98759401271422391785, + .97682326283354439220, + -.95915923302922997013, + .93198380256105393618, + -.89273040299591077603, + .83945793215750220154, + -.77148960729470505477, + .68992761745934847866, + -.59784149990330073143, + .50000000000000000000, + -.40215850009669926857, + .31007238254065152134, + -.22851039270529494523, + .16054206784249779846, + -.10726959700408922397, + .68016197438946063823e-1, + -.40840766970770029873e-1, + .23176737166455607805e-1, + -.12405987285776082154e-1, + .62512730682449246388e-2, + -.29592166263096543401e-2, + .13131884053420191908e-2, + -.54504190222378945440e-3, + .21106516173760261250e-3, + -.76053287924037718971e-4, + .25423835243950883896e-4, + -.78585149263697370338e-5, + .22376124247437700378e-5, + -.58440580661848562719e-6, + .13931171712321674741e-6, + -.30137695171547022183e-7, + .58768014045093054654e-8, + -.10246226511017621219e-8, + .15816215942184366772e-9, + -.21352608103961806529e-10, + .24823251635643084345e-11, + -.24347803504257137241e-12, + .19593322190397666205e-13, + -.12421162189080181548e-14, + .58167446553847312884e-16, + -.17889335846010823161e-17, + .27105054312137610850e-19] + +function eta(z::Union(Float64,Complex128)) + if z == 0 + return oftype(z, 0.5) + end + re, im = reim(z) + if im==0 && re < 0 && re==round(re/2)*2 + return zero(z) + end + reflect = false + if re < 0.5 + z = 1-z + reflect = true + end + s = zero(z) + for n = length(eta_coeffs):-1:1 + c = eta_coeffs[n] + p = n^-z + s += c * p + end + if reflect + z2 = 2.0^z + b = 2.0 - (2.0*z2) + f = z2 - 2 + piz = pi^z + + b = b/f/piz + + return s * gamma(z) * b * cospi(z/2) + end + return s +end + +eta(x::Integer) = eta(float64(x)) +eta(x::Real) = oftype(x,eta(float64(x))) +eta(z::Complex) = oftype(z,eta(complex128(z))) +@vectorize_1arg Number eta + +function zeta(z::Number) + zz = 2^z + eta(z) * zz/(zz-2) +end +@vectorize_1arg Number zeta diff -Nru julia-0.3.0/base/special/trig.jl julia-0.3.0/base/special/trig.jl --- julia-0.3.0/base/special/trig.jl 1970-01-01 00:00:00.000000000 +0000 +++ julia-0.3.0/base/special/trig.jl 2014-06-03 19:33:18.000000000 +0000 @@ -0,0 +1,184 @@ +function sinpi(x::Real) + if isinf(x) + return throw(DomainError()) + elseif isnan(x) + return nan(x) + end + + rx = copysign(float(rem(x,2)),x) + arx = abs(rx) + + if arx < oftype(rx,0.25) + return sin(pi*rx) + elseif arx <= oftype(rx,0.75) + arx = oftype(rx,0.5) - arx + return copysign(cos(pi*arx),rx) + elseif arx < oftype(rx,1.25) + rx = (one(rx) - arx)*sign(rx) + return sin(pi*rx) + elseif arx <= oftype(rx,1.75) + arx = oftype(rx,1.5) - arx + return -copysign(cos(pi*arx),rx) + else + rx = rx - copysign(oftype(rx,2.0),rx) + return sin(pi*rx) + end +end + +function cospi(x::Real) + if isinf(x) + return throw(DomainError()) + elseif isnan(x) + return nan(x) + end + + rx = abs(float(rem(x,2))) + + if rx <= oftype(rx,0.25) + return cos(pi*rx) + elseif rx < oftype(rx,0.75) + rx = oftype(rx,0.5) - rx + return sin(pi*rx) + elseif rx <= oftype(rx,1.25) + rx = one(rx) - rx + return -cos(pi*rx) + elseif rx < oftype(rx,1.75) + rx = rx - oftype(rx,1.5) + return sin(pi*rx) + else + rx = oftype(rx,2.0) - rx + return cos(pi*rx) + end +end + +sinpi(x::Integer) = zero(x) +cospi(x::Integer) = isodd(x) ? -one(x) : one(x) + +function sinpi(z::Complex) + zr, zi = reim(z) + if !isfinite(zi) && zr == 0 return complex(zr, zi) end + if isnan(zr) && !isfinite(zi) return complex(zr, zi) end + if !isfinite(zr) && zi == 0 return complex(oftype(zr, NaN), zi) end + if !isfinite(zr) && isfinite(zi) return complex(oftype(zr, NaN), oftype(zi, NaN)) end + if !isfinite(zr) && !isfinite(zi) return complex(zr, oftype(zi, NaN)) end + pizi = pi*zi + complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi)) +end + +function cospi(z::Complex) + zr, zi = reim(z) + if !isfinite(zi) && zr == 0 + return complex(isnan(zi) ? zi : oftype(zi, Inf), + isnan(zi) ? zr : zr*-sign(zi)) + end + if !isfinite(zr) && isinf(zi) + return complex(oftype(zr, Inf), oftype(zi, NaN)) + end + if isinf(zr) + return complex(oftype(zr, NaN), zi==0 ? -copysign(zi, zr) : oftype(zi, NaN)) + end + if isnan(zr) && zi==0 return complex(zr, abs(zi)) end + pizi = pi*zi + complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi)) +end +@vectorize_1arg Number sinpi +@vectorize_1arg Number cospi + + +sinc(x::Number) = x==0 ? one(x) : oftype(x,sinpi(x)/(pi*x)) +sinc(x::Integer) = x==0 ? one(x) : zero(x) +sinc{T<:Integer}(x::Complex{T}) = sinc(float(x)) +@vectorize_1arg Number sinc +cosc(x::Number) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x) +cosc(x::Integer) = cosc(float(x)) +cosc{T<:Integer}(x::Complex{T}) = cosc(float(x)) +@vectorize_1arg Number cosc + +for (finv, f) in ((:sec, :cos), (:csc, :sin), (:cot, :tan), + (:sech, :cosh), (:csch, :sinh), (:coth, :tanh), + (:secd, :cosd), (:cscd, :sind), (:cotd, :tand)) + @eval begin + ($finv){T<:Number}(z::T) = one(T) / (($f)(z)) + ($finv){T<:Number}(z::AbstractArray{T}) = one(T) ./ (($f)(z)) + end +end + +for (fa, fainv) in ((:asec, :acos), (:acsc, :asin), (:acot, :atan), + (:asech, :acosh), (:acsch, :asinh), (:acoth, :atanh)) + @eval begin + ($fa){T<:Number}(y::T) = ($fainv)(one(T) / y) + ($fa){T<:Number}(y::AbstractArray{T}) = ($fainv)(one(T) ./ y) + end +end + +function sind(x::Real) + if isinf(x) + return throw(DomainError()) + elseif isnan(x) + return nan(x) + end + + rx = copysign(float(rem(x,360)),x) + arx = abs(rx) + + if arx < oftype(rx,45.0) + return sin(deg2rad(rx)) + elseif arx <= oftype(rx,135.0) + arx = oftype(rx,90.0) - arx + return copysign(cos(deg2rad(arx)),rx) + elseif arx < oftype(rx,225.0) + rx = (oftype(rx,180.0) - arx)*sign(rx) + return sin(deg2rad(rx)) + elseif arx <= 315.0 + arx = oftype(rx,270.0) - arx + return -copysign(cos(deg2rad(arx)),rx) + else + rx = rx - copysign(oftype(rx,360.0),rx) + return sin(deg2rad(rx)) + end +end +@vectorize_1arg Real sind + +function cosd(x::Real) + if isinf(x) + return throw(DomainError()) + elseif isnan(x) + return nan(x) + end + + rx = abs(float(rem(x,360))) + + if rx <= oftype(rx,45.0) + return cos(deg2rad(rx)) + elseif rx < oftype(rx,135.0) + rx = oftype(rx,90.0) - rx + return sin(deg2rad(rx)) + elseif rx <= oftype(rx,225.0) + rx = oftype(rx,180.0) - rx + return -cos(deg2rad(rx)) + elseif rx < oftype(rx,315.0) + rx = rx - oftype(rx,270.0) + return sin(deg2rad(rx)) + else + rx = oftype(rx,360.0) - rx + return cos(deg2rad(rx)) + end +end +@vectorize_1arg Real cosd + +tand(x::Real) = sind(x) / cosd(x) +@vectorize_1arg Real tand + +for (fd, f) in ((:sind, :sin), (:cosd, :cos), (:tand, :tan)) + @eval begin + ($fd)(z) = ($f)(deg2rad(z)) + end +end + +for (fd, f) in ((:asind, :asin), (:acosd, :acos), (:atand, :atan), + (:asecd, :asec), (:acscd, :acsc), (:acotd, :acot)) + @eval begin + ($fd)(y) = rad2deg(($f)(y)) + @vectorize_1arg Real $fd + end +end diff -Nru julia-0.3.0/base/version_git.jl julia-0.3.0/base/version_git.jl --- julia-0.3.0/base/version_git.jl 2014-06-02 19:54:52.000000000 +0000 +++ julia-0.3.0/base/version_git.jl 2014-06-03 19:34:20.000000000 +0000 @@ -11,12 +11,12 @@ end const GIT_VERSION_INFO = GitVersionInfo( - "17b6262464fb935378ba67ca25b8bf96897a16c2", - "17b6262*", + "06b809c848343402a9a591b65a076ec0a7f9f30a", + "06b809c*", "master", - 3363, - "2014-06-02 04:12 UTC", + 3424, + "2014-06-03 04:42 UTC", false, 0, - 1401682354. + 1401770534. ) diff -Nru julia-0.3.0/base/version_git.jl.backup julia-0.3.0/base/version_git.jl.backup --- julia-0.3.0/base/version_git.jl.backup 2014-06-02 19:54:51.000000000 +0000 +++ julia-0.3.0/base/version_git.jl.backup 2014-06-03 19:34:20.000000000 +0000 @@ -11,12 +11,12 @@ end const GIT_VERSION_INFO = GitVersionInfo( - "17b6262464fb935378ba67ca25b8bf96897a16c2", - "17b6262*", + "06b809c848343402a9a591b65a076ec0a7f9f30a", + "06b809c*", "master", - 3363, - "2014-06-02 04:12 UTC", + 3424, + "2014-06-03 04:42 UTC", false, 0, - 1401682354. + 1401770534. ) diff -Nru julia-0.3.0/debian/bzr-builder.manifest julia-0.3.0/debian/bzr-builder.manifest --- julia-0.3.0/debian/bzr-builder.manifest 2014-06-02 19:54:05.000000000 +0000 +++ julia-0.3.0/debian/bzr-builder.manifest 2014-06-03 19:33:21.000000000 +0000 @@ -1,2 +1,2 @@ -# bzr-builder format 0.3 deb-version {debupstream}-1082 -lp:~staticfloat/julianightlies/trunk revid:staticfloat@gmail.com-20140602120013-u6frsa6yzqip4h1y +# bzr-builder format 0.3 deb-version {debupstream}-1085 +lp:~staticfloat/julianightlies/trunk revid:staticfloat@gmail.com-20140603120012-pgg2fsj009z9oo23 diff -Nru julia-0.3.0/debian/changelog julia-0.3.0/debian/changelog --- julia-0.3.0/debian/changelog 2014-06-02 19:54:05.000000000 +0000 +++ julia-0.3.0/debian/changelog 2014-06-03 19:33:21.000000000 +0000 @@ -1,14 +1,14 @@ -julia (0.3.0-1082~ubuntu12.10.1) quantal; urgency=low +julia (0.3.0-1085~ubuntu12.10.1) quantal; urgency=low * Auto build. - -- Elliot Saba Mon, 02 Jun 2014 19:54:05 +0000 + -- Elliot Saba Tue, 03 Jun 2014 19:33:21 +0000 -julia (0.3.0-prerelease+20140602) precise; urgency=low +julia (0.3.0-prerelease+20140603) precise; urgency=low * nightly git build - -- Elliot Saba Mon, 02 Jun 2014 08:00:13 -0400 + -- Elliot Saba Tue, 03 Jun 2014 08:00:11 -0400 julia (0.2.0-pre) unstable; urgency=low diff -Nru julia-0.3.0/doc/helpdb.jl julia-0.3.0/doc/helpdb.jl --- julia-0.3.0/doc/helpdb.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/doc/helpdb.jl 2014-06-03 19:33:18.000000000 +0000 @@ -2864,7 +2864,7 @@ "), -("Base","readdlm","readdlm(source, delim::Char, T::Type, eol::Char; has_header=false, use_mmap, ignore_invalid_chars=false, quotes=true, dims, comments=true, comment_char='#') +("Base","readdlm","readdlm(source, delim::Char, T::Type, eol::Char; header=false, skipstart=0, use_mmap, ignore_invalid_chars=false, quotes=true, dims, comments=true, comment_char='#') Read a matrix from the source where each line (separated by \"eol\") gives one row, with elements separated by the given @@ -2877,10 +2877,13 @@ or zero. Other useful values of \"T\" include \"ASCIIString\", \"String\", and \"Any\". - If \"has_header\" is \"true\", the first row of data would be read - as headers and the tuple \"(data_cells, header_cells)\" is returned + If \"header\" is \"true\", the first row of data will be read + as header and the tuple \"(data_cells, header_cells)\" is returned instead of only \"data_cells\". + Specifying \"skipstart\" will ignore the corresponding number + of initial lines from the input. + If \"use_mmap\" is \"true\", the file specified by \"source\" is memory mapped for potential speedups. Default is \"true\" except on Windows. On Windows, you may want to specify \"true\" if the file diff -Nru julia-0.3.0/doc/manual/variables-and-scoping.rst julia-0.3.0/doc/manual/variables-and-scoping.rst --- julia-0.3.0/doc/manual/variables-and-scoping.rst 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/doc/manual/variables-and-scoping.rst 2014-06-03 19:33:18.000000000 +0000 @@ -167,6 +167,24 @@ prompt this declaration would be necessary in order to modify a global variable. +Multiple variables can be declared global using the following syntax:: + + function foo() + global x=1, y="bar", z=3 + end + + julia> foo() + 3 + + julia> x + 1 + + julia> y + "bar" + + julia> z + 3 + The ``let`` statement provides a different way to introduce variables. Unlike assignments to local variables, ``let`` statements allocate new variable bindings each time they run. An assignment modifies an existing diff -Nru julia-0.3.0/doc/stdlib/base.rst julia-0.3.0/doc/stdlib/base.rst --- julia-0.3.0/doc/stdlib/base.rst 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/doc/stdlib/base.rst 2014-06-03 19:33:18.000000000 +0000 @@ -1899,13 +1899,15 @@ Create an iterable object that will yield each line from a stream. -.. function:: readdlm(source, delim::Char, T::Type, eol::Char; has_header=false, use_mmap, ignore_invalid_chars=false, quotes=true, dims, comments=true, comment_char='#') +.. function:: readdlm(source, delim::Char, T::Type, eol::Char; header=false, skipstart=0, use_mmap, ignore_invalid_chars=false, quotes=true, dims, comments=true, comment_char='#') Read a matrix from the source where each line (separated by ``eol``) gives one row, with elements separated by the given delimeter. The source can be a text file, stream or byte array. Memory mapped files can be used by passing the byte array representation of the mapped segment as source. If ``T`` is a numeric type, the result is an array of that type, with any non-numeric elements as ``NaN`` for floating-point types, or zero. Other useful values of ``T`` include ``ASCIIString``, ``String``, and ``Any``. - If ``has_header`` is ``true``, the first row of data would be read as headers and the tuple ``(data_cells, header_cells)`` is returned instead of only ``data_cells``. + If ``header`` is ``true``, the first row of data will be read as header and the tuple ``(data_cells, header_cells)`` is returned instead of only ``data_cells``. + + Specifying ``skipstart`` will ignore the corresponding number of initial lines from the input. If ``use_mmap`` is ``true``, the file specified by ``source`` is memory mapped for potential speedups. Default is ``true`` except on Windows. On Windows, you may want to specify ``true`` if the file is large, and is only read once and not written to. diff -Nru julia-0.3.0/doc/stdlib/sort.rst julia-0.3.0/doc/stdlib/sort.rst --- julia-0.3.0/doc/stdlib/sort.rst 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/doc/stdlib/sort.rst 2014-06-03 19:33:18.000000000 +0000 @@ -96,7 +96,7 @@ Sorting Functions ----------------- -.. function:: sort!(v, [dim,] [alg=,] [by=,] [lt=,] [rev=false]) +.. function:: sort!(v, [alg=,] [by=,] [lt=,] [rev=false]) Sort the vector ``v`` in place. ``QuickSort`` is used by default for numeric arrays while ``MergeSort`` is used for other arrays. You can specify an algorithm to use via diff -Nru julia-0.3.0/src/ccall.cpp julia-0.3.0/src/ccall.cpp --- julia-0.3.0/src/ccall.cpp 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/ccall.cpp 2014-06-03 19:33:18.000000000 +0000 @@ -752,6 +752,18 @@ builder.CreateBitCast(emit_nthptr_addr(ary, addressOf?1:0), lrt), rt); } + if (fptr == &jl_is_leaf_type || + (f_lib==NULL && f_name && !strcmp(f_name, "jl_is_leaf_type"))) { + jl_value_t *arg = args[4]; + if (jl_is_expr(arg) || jl_is_symbol(arg)) { + jl_value_t* ty = expr_type(arg, ctx); + if (jl_is_type_type(ty) && !jl_is_typevar(jl_tparam0(ty))) { + int isleaf = jl_is_leaf_type(jl_tparam0(ty)); + JL_GC_POP(); + return ConstantInt::get(T_int32, isleaf); + } + } + } // save place before arguments, for possible insertion of temp arg // area saving code. diff -Nru julia-0.3.0/src/codegen.cpp julia-0.3.0/src/codegen.cpp --- julia-0.3.0/src/codegen.cpp 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/codegen.cpp 2014-06-03 19:33:18.000000000 +0000 @@ -1616,6 +1616,16 @@ } } } + else if (f->fptr == &jl_f_subtype && nargs == 2) { + rt1 = expr_type(args[1], ctx); + rt2 = expr_type(args[2], ctx); + if (jl_is_type_type(rt1) && !jl_is_typevar(jl_tparam0(rt1)) && + jl_is_type_type(rt2) && !jl_is_typevar(jl_tparam0(rt2))) { + int issub = jl_subtype(jl_tparam0(rt1), jl_tparam0(rt2), 0); + JL_GC_POP(); + return ConstantInt::get(T_int1, issub); + } + } else if (f->fptr == &jl_f_tuplelen && nargs==1) { jl_value_t *aty = expr_type(args[1], ctx); rt1 = aty; if (jl_is_tuple(aty)) { @@ -2042,6 +2052,10 @@ &jl_tupleref(ctx->sp,0), jl_tuple_len(ctx->sp)/2); if (jl_is_leaf_type(ty)) { + if (jl_has_typevars(ty)) { + // add root for types not cached. issue #7065 + jl_add_linfo_root(ctx->linfo, ty); + } JL_GC_POP(); return literal_pointer_val(ty); } diff -Nru julia-0.3.0/src/gc.c julia-0.3.0/src/gc.c --- julia-0.3.0/src/gc.c 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/gc.c 2014-06-03 19:33:18.000000000 +0000 @@ -10,9 +10,9 @@ #include "julia_internal.h" #ifdef _P64 -#define GC_PAGE_SZ (1536*sizeof(void*))//bytes +#define GC_PAGE_SZ 12288//bytes #else -#define GC_PAGE_SZ (2048*sizeof(void*))//bytes +#define GC_PAGE_SZ 8192//bytes #endif #ifdef __cplusplus @@ -65,13 +65,14 @@ static size_t allocd_bytes = 0; static int64_t total_allocd_bytes = 0; static size_t freed_bytes = 0; -#define default_collect_interval (3200*1024*sizeof(void*)) -static size_t collect_interval = default_collect_interval; #ifdef _P64 +#define default_collect_interval (5600*1024*sizeof(void*)) static size_t max_collect_interval = 1250000000UL; #else -static size_t max_collect_interval = 500000000UL; +#define default_collect_interval (3200*1024*sizeof(void*)) +static size_t max_collect_interval = 500000000UL; #endif +static size_t collect_interval = default_collect_interval; int jl_in_gc; // referenced from switchto task.c #ifdef OBJPROFILE @@ -275,7 +276,7 @@ } JL_CATCH { JL_PRINTF(JL_STDERR, "error in running finalizer: "); - jl_show(jl_stderr_obj(), jl_exception_in_transit); + jl_static_show(JL_STDERR, jl_exception_in_transit); JL_PUTC('\n',JL_STDERR); } ff = jl_t1(ff); @@ -287,7 +288,7 @@ } JL_CATCH { JL_PRINTF(JL_STDERR, "error in running finalizer: "); - jl_show(jl_stderr_obj(), jl_exception_in_transit); + jl_static_show(JL_STDERR, jl_exception_in_transit); JL_PUTC('\n',JL_STDERR); } } diff -Nru julia-0.3.0/src/intrinsics.cpp julia-0.3.0/src/intrinsics.cpp --- julia-0.3.0/src/intrinsics.cpp 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/intrinsics.cpp 2014-06-03 19:33:18.000000000 +0000 @@ -473,26 +473,40 @@ return builder.CreateZExt(JL_INT(auto_unbox(x,ctx)), to); } -static Value *emit_eqfsi64(Value *x, Value *y) +static Value *emit_eqfsi(Value *x, Value *y) { x = FP(x); Value *fy = JL_INT(y); + + // using all 64-bit is slightly faster than using mixed sizes + Value *xx = x, *vv = fy; + if (x->getType() == T_float32) + xx = builder.CreateFPExt(xx, T_float64); + if (vv->getType()->getPrimitiveSizeInBits() < 64) + vv = builder.CreateSExt(vv, T_int64); + + Value *back = builder.CreateSIToFP(vv, xx->getType()); return builder.CreateAnd - (builder.CreateFCmpOEQ(x, builder.CreateSIToFP(fy, T_float64)), - builder.CreateICmpEQ(fy, builder.CreateFPToSI - (builder.CreateSIToFP(fy, T_float64), - T_int64))); + (builder.CreateFCmpOEQ(xx, back), + builder.CreateICmpEQ(vv, builder.CreateFPToSI(back, vv->getType()))); } -static Value *emit_eqfui64(Value *x, Value *y) +static Value *emit_eqfui(Value *x, Value *y) { x = FP(x); Value *fy = JL_INT(y); + + // using all 64-bit is slightly faster than using mixed sizes + Value *xx = x, *vv = fy; + if (x->getType() == T_float32) + xx = builder.CreateFPExt(xx, T_float64); + if (vv->getType()->getPrimitiveSizeInBits() < 64) + vv = builder.CreateZExt(vv, T_int64); + + Value *back = builder.CreateUIToFP(vv, xx->getType()); return builder.CreateAnd - (builder.CreateFCmpOEQ(x, builder.CreateUIToFP(fy, T_float64)), - builder.CreateICmpEQ(fy, builder.CreateFPToUI - (builder.CreateUIToFP(fy, T_float64), - T_int64))); + (builder.CreateFCmpOEQ(xx, back), + builder.CreateICmpEQ(vv, builder.CreateFPToUI(back, vv->getType()))); } static Value *emit_checked_fptosi(Type *to, Value *x, jl_codectx_t *ctx) @@ -506,12 +520,7 @@ prepare_global(jlinexacterr_var), ctx); } else { - Value *xx = x, *vv = v; - if (x->getType() == T_float32) - xx = builder.CreateFPExt(x, T_float64); - if (to->getPrimitiveSizeInBits() < 64) - vv = builder.CreateSExt(v, T_int64); - raise_exception_unless(emit_eqfsi64(xx, vv), prepare_global(jlinexacterr_var), ctx); + raise_exception_unless(emit_eqfsi(x, v), prepare_global(jlinexacterr_var), ctx); } return v; } @@ -532,12 +541,7 @@ prepare_global(jlinexacterr_var), ctx); } else { - Value *xx = x, *vv = v; - if (x->getType() == T_float32) - xx = builder.CreateFPExt(x, T_float64); - if (to->getPrimitiveSizeInBits() < 64) - vv = builder.CreateZExt(v, T_int64); - raise_exception_unless(emit_eqfui64(xx, vv), prepare_global(jlinexacterr_var), ctx); + raise_exception_unless(emit_eqfui(x, v), prepare_global(jlinexacterr_var), ctx); } return v; } @@ -796,6 +800,32 @@ } HANDLE(uitofp,2) return builder.CreateUIToFP(JL_INT(auto_unbox(args[2],ctx)), FTnbits(try_to_determine_bitstype_nbits(args[1],ctx))); HANDLE(sitofp,2) return builder.CreateSIToFP(JL_INT(auto_unbox(args[2],ctx)), FTnbits(try_to_determine_bitstype_nbits(args[1],ctx))); + + case fptoui: + if (nargs == 1) { + Value *x = FP(auto_unbox(args[1], ctx)); + return builder.CreateFPToUI(FP(x), JL_INTT(x->getType())); + } + else if (nargs == 2) { + return builder.CreateFPToUI(FP(auto_unbox(args[2],ctx)), + Type::getIntNTy(jl_LLVMContext, try_to_determine_bitstype_nbits(args[1],ctx))); + } + else { + jl_error("fptoui: wrong number of arguments"); + } + case fptosi: + if (nargs == 1) { + Value *x = FP(auto_unbox(args[1], ctx)); + return builder.CreateFPToSI(FP(x), JL_INTT(x->getType())); + } + else if (nargs == 2) { + return builder.CreateFPToSI(FP(auto_unbox(args[2],ctx)), + Type::getIntNTy(jl_LLVMContext, try_to_determine_bitstype_nbits(args[1],ctx))); + } + else { + jl_error("fptosi: wrong number of arguments"); + } + HANDLE(fptrunc,2) return builder.CreateFPTrunc(FP(auto_unbox(args[2],ctx)), FTnbits(try_to_determine_bitstype_nbits(args[1],ctx))); HANDLE(fpext,2) { // when extending a float32 to a float64, we need to force @@ -946,8 +976,8 @@ HANDLE(lt_float,2) return builder.CreateFCmpOLT(FP(x), FP(y)); HANDLE(le_float,2) return builder.CreateFCmpOLE(FP(x), FP(y)); - HANDLE(eqfsi64,2) return emit_eqfsi64(x, y); - HANDLE(eqfui64,2) return emit_eqfui64(x, y); + HANDLE(eqfsi64,2) return emit_eqfsi(x, y); + HANDLE(eqfui64,2) return emit_eqfui(x, y); HANDLE(ltfsi64,2) { x = FP(x); fy = JL_INT(y); @@ -1175,8 +1205,6 @@ } #endif - HANDLE(fptoui,1) return builder.CreateFPToUI(FP(x), JL_INTT(x->getType())); - HANDLE(fptosi,1) return builder.CreateFPToSI(FP(x), JL_INTT(x->getType())); HANDLE(fpsiround,1) HANDLE(fpuiround,1) { diff -Nru julia-0.3.0/src/jl_uv.c julia-0.3.0/src/jl_uv.c --- julia-0.3.0/src/jl_uv.c 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/jl_uv.c 2014-06-03 19:33:18.000000000 +0000 @@ -309,18 +309,6 @@ if ((handle->type == UV_NAMED_PIPE || handle->type == UV_TCP) && uv_is_writable((uv_stream_t*)handle)) { - // Make sure that the stream has not already been marked closed in Julia. - // A double shutdown would cause the process to hang on exit. - if (handle->data) { - JULIA_CB(isopen, handle->data, 0); - if (!jl_is_int32(ret)) { - jl_error("jl_close_uv: _uv_hook_isopen must return an int32."); - } - if (!jl_unbox_int32(ret)){ - return; - } - } - uv_shutdown_t *req = (uv_shutdown_t*)malloc(sizeof(uv_shutdown_t)); req->data = 0; /* diff -Nru julia-0.3.0/src/llvm-simdloop.cpp julia-0.3.0/src/llvm-simdloop.cpp --- julia-0.3.0/src/llvm-simdloop.cpp 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/src/llvm-simdloop.cpp 2014-06-03 19:33:18.000000000 +0000 @@ -82,8 +82,13 @@ for (Instruction *I = Phi; ; I=J) { J = NULL; // Find the user of instruction I that is within loop L. +#if defined(LLVM_VERSION_MAJOR) && LLVM_VERSION_MAJOR == 3 && LLVM_VERSION_MINOR >= 5 + for (User *UI : I->users()) { /*}*/ + Instruction *U = cast(UI); +#else for (Value::use_iterator UI = I->use_begin(), UE = I->use_end(); UI != UE; ++UI) { Instruction *U = cast(*UI); +#endif if (L->contains(U)) { if (J) { DEBUG(dbgs() << "LSL: not a reduction var because op has two internal uses: " << *I << "\n"); diff -Nru julia-0.3.0/test/collections.jl julia-0.3.0/test/collections.jl --- julia-0.3.0/test/collections.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/collections.jl 2014-06-03 19:33:18.000000000 +0000 @@ -171,6 +171,39 @@ @test d == {8=>19, 19=>2, 42=>4} end +# show +for d in (["\n" => "\n", "1" => "\n", "\n" => "2"], + [string(i) => i for i = 1:30], + [reshape(1:i^2,i,i) => reshape(1:i^2,i,i) for i = 1:24], + [utf8(Char['α':'α'+i]) => utf8(Char['α':'α'+i]) for i = (1:10)*10]) + for cols in (12, 40, 80), rows in (2, 10, 24) + # Ensure output is limited as requested + s = IOBuffer() + Base.showdict(s, d, true, rows, cols) + out = split(takebuf_string(s),'\n') + for line in out[2:end] + @test strwidth(line) <= cols + end + @test length(out) <= rows + + for f in (keys, values) + s = IOBuffer() + Base.showkv(s, f(d), true, rows, cols) + out = split(takebuf_string(s),'\n') + for line in out[2:end] + @test strwidth(line) <= cols + end + @test length(out) <= rows + end + end + # Simply ensure these do not throw errors + Base.showdict(IOBuffer(), d, false) + @test !isempty(summary(d)) + @test !isempty(summary(keys(d))) + @test !isempty(summary(values(d))) +end + + # issue #2540 d = {x => 1 for x in ['a', 'b', 'c']} diff -Nru julia-0.3.0/test/combinatorics.jl julia-0.3.0/test/combinatorics.jl --- julia-0.3.0/test/combinatorics.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/combinatorics.jl 2014-06-03 19:33:18.000000000 +0000 @@ -15,8 +15,12 @@ @test collect(filter(x->(iseven(x[1])),combinations([1,2,3],2))) == {[2,3]} @test collect(partitions(4)) == {[4], [3,1], [2,2], [2,1,1], [1,1,1,1]} @test collect(partitions(8,3)) == {[6,1,1], [5,2,1], [4,3,1], [4,2,2], [3,3,2]} +@test collect(partitions(8, 1)) == {[8]} +@test collect(partitions(8, 9)) == {} @test collect(partitions([1,2,3])) == {{[1,2,3]}, {[1,2],[3]}, {[1,3],[2]}, {[1],[2,3]}, {[1],[2],[3]}} @test collect(partitions([1,2,3,4],3)) == {{[1,2],[3],[4]}, {[1,3],[2],[4]}, {[1],[2,3],[4]}, {[1,4],[2],[3]}, {[1],[2,4],[3]},{[1],[2],[3,4]}} +@test collect(partitions([1,2,3,4],1)) == {{[1, 2, 3, 4]}} +@test collect(partitions([1,2,3,4],5)) == {} @test length(collect(partitions(30))) == length(partitions(30)) @test length(collect(partitions(90,4))) == length(partitions(90,4)) @@ -32,7 +36,7 @@ #Issue 6154 @test binomial(int32(34), int32(15)) == binomial(BigInt(34), BigInt(15)) == 1855967520 @test binomial(int64(67), int64(29)) == binomial(BigInt(67), BigInt(29)) == 7886597962249166160 -@test binomial(int128(131), int128(62)) == binomial(BigInt(131), BigInt(62)) == 157311720980559117816198361912717812000 +@test binomial(int128(131), int128(62)) == binomial(BigInt(131), BigInt(62)) == 157311720980559117816198361912717812000 # issue #6579 @test factorial(0) == 1 diff -Nru julia-0.3.0/test/linalg2.jl julia-0.3.0/test/linalg2.jl --- julia-0.3.0/test/linalg2.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/linalg2.jl 2014-06-03 19:33:18.000000000 +0000 @@ -3,7 +3,6 @@ import Base.LinAlg import Base.LinAlg: BlasComplex, BlasFloat, BlasReal - # basic tridiagonal operations n = 5 @@ -191,248 +190,6 @@ @test_approx_eq vals LinAlg.LAPACK.syev!('N','U',copy(Asym)) end -#Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences -function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) - n = size(a, 1) - @test n==size(b,1) && size(a,2)==size(b,2) - error==nothing && (error=n^3*(eps(S)+eps(T))) - for i=1:n - ev1, ev2 = a[:,i], b[:,i] - deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) - if !isnan(deviation) - @test_approx_eq_eps deviation 0.0 error - end - end -end - -############################## -# Tests for special matrices # -############################## - -#Triangular matrices -n=12 -for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) - A = convert(Matrix{elty}, randn(n, n)) - b = convert(Matrix{elty}, randn(n, 2)) - if elty <: Complex - A += im*convert(Matrix{elty}, randn(n, n)) - b += im*convert(Matrix{elty}, randn(n, 2)) - end - - for M in (triu(A), tril(A)) - TM = Triangular(M) - #Multiplication - @test_approx_eq A*M A*TM - @test_approx_eq M*A TM*A - @test_approx_eq M*M TM*TM - @test_approx_eq M*b TM*b - condM = elty <:BlasFloat ? cond(TM, Inf) : convert(relty, cond(complex128(M), Inf)) - #Linear solver - x = M \ b - tx = TM \ b - @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) - if elty <: BlasFloat #test naivesub! against LAPACK - tx = [LinAlg.naivesub!(TM, b[:,1]) LinAlg.naivesub!(TM, b[:,2])] - @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) - end - - #Eigensystems - vals1, vecs1 = eig(complex128(M)) - vals2, vecs2 = eig(TM) - res1=norm(complex128(vecs1*diagm(vals1)*inv(vecs1) - M)) - res2=norm(complex128(vecs2*diagm(vals2)*inv(vecs2) - full(TM))) - @test_approx_eq_eps res1 res2 res1+res2 - - if elty <:BlasFloat - #Condition number tests - can be VERY approximate - for p in [1.0, Inf] - @test_approx_eq_eps cond(TM, p) cond(M, p) (cond(TM,p)+cond(M,p))*0.2 - end - end - end -end - -#Tridiagonal matrices -for relty in (Float32, Float64), elty in (relty, Complex{relty}) - a = convert(Vector{elty}, randn(n-1)) - b = convert(Vector{elty}, randn(n)) - c = convert(Vector{elty}, randn(n-1)) - if elty <: Complex - a += im*convert(Vector{elty}, randn(n-1)) - b += im*convert(Vector{elty}, randn(n)) - c += im*convert(Vector{elty}, randn(n-1)) - end - - A=Tridiagonal(a, b, c) - fA=(elty<:Complex?complex128:float64)(full(A)) - for func in (det, inv) - @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) - end -end - -#SymTridiagonal (symmetric tridiagonal) matrices -for relty in (Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work - a = convert(Vector{elty}, randn(n)) - b = convert(Vector{elty}, randn(n-1)) - if elty <: Complex - a += im*convert(Vector{elty}, randn(n)) - b += im*convert(Vector{elty}, randn(n-1)) - end - - A=SymTridiagonal(a, b) - fA=(elty<:Complex?complex128:float64)(full(A)) - for func in (det, inv) - @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) - end -end - -Ainit = randn(n) -Binit = randn(n-1) -for elty in (Float32, Float64) - A = convert(Array{elty, 1}, Ainit) - B = convert(Array{elty, 1}, Binit) - zero, infinity = convert(elty, 0), convert(elty, Inf) - #This tests eigenvalue and eigenvector computations using stebz! and stein! - w, iblock, isplit = LinAlg.LAPACK.stebz!('V','B',-infinity,infinity,0,0,zero,A,B) - evecs = LinAlg.LAPACK.stein!(A,B,w) - - (e, v)=eig(SymTridiagonal(A,B)) - @test_approx_eq e w - #Take into account possible phase (sign) difference in eigenvectors - for i=1:n - ev1 = v[:,i] - ev2 = evecs[:,i] - deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) - @test_approx_eq_eps deviation 0.0 n^2*eps(abs(convert(elty, 2.0))) - end - - #Test stein! call using iblock and isplit - w, iblock, isplit = LinAlg.LAPACK.stebz!('V','B',-infinity,infinity,0,0,zero,A,B) - evecs = LinAlg.LAPACK.stein!(A, B, w, iblock, isplit) - test_approx_eq_vecs(v, evecs) -end - -#Bidiagonal matrices -for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) - dv = convert(Vector{elty}, randn(n)) - ev = convert(Vector{elty}, randn(n-1)) - b = convert(Matrix{elty}, randn(n, 2)) - if (elty <: Complex) - dv += im*convert(Vector{elty}, randn(n)) - ev += im*convert(Vector{elty}, randn(n-1)) - b += im*convert(Matrix{elty}, randn(n, 2)) - end - for isupper in (true, false) #Test upper and lower bidiagonal matrices - T = Bidiagonal(dv, ev, isupper) - - @test size(T, 1) == size(T, 2) == n - @test size(T) == (n, n) - @test full(T) == diagm(dv) + diagm(ev, isupper?1:-1) - @test Bidiagonal(full(T), isupper) == T - z = zeros(elty, n) - - # idempotent tests - for func in (conj, transpose, ctranspose) - @test func(func(T)) == T - end - - #Linear solver - Tfull = full(T) - condT = cond(complex128(Tfull)) - x = T \ b - tx = Tfull \ b - @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) - - #Test eigenvalues/vectors - d1, v1 = eig(T) - d2, v2 = eig((elty<:Complex?complex128:float64)(Tfull)) - @test_approx_eq isupper?d1:reverse(d1) d2 - if elty <: Real - test_approx_eq_vecs(v1, isupper?v2:v2[:,n:-1:1]) - end - - if (elty <: BlasReal) - #Test singular values/vectors - @test_approx_eq svdvals(Tfull) svdvals(T) - u1, d1, v1 = svd(Tfull) - u2, d2, v2 = svd(T) - @test_approx_eq d1 d2 - if elty <: Real - test_approx_eq_vecs(u1, u2) - test_approx_eq_vecs(v1, v2) - end - @test_approx_eq_eps 0 vecnorm(u2*diagm(d2)*v2'-Tfull) n*max(n^2*eps(relty), vecnorm(u1*diagm(d1)*v1'-Tfull)) - end - end -end - -#Diagonal matrices -n=12 -for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) - d=convert(Vector{elty}, randn(n)) - v=convert(Vector{elty}, randn(n)) - U=convert(Matrix{elty}, randn(n,n)) - if elty <: Complex - d+=im*convert(Vector{elty}, randn(n)) - v+=im*convert(Vector{elty}, randn(n)) - U+=im*convert(Matrix{elty}, randn(n,n)) - end - D = Diagonal(d) - DM = diagm(d) - @test_approx_eq_eps D*v DM*v n*eps(relty)*(elty<:Complex ? 2:1) - @test_approx_eq_eps D*U DM*U n^2*eps(relty)*(elty<:Complex ? 2:1) - if relty != BigFloat - @test_approx_eq_eps D\v DM\v 2n^2*eps(relty)*(elty<:Complex ? 2:1) - @test_approx_eq_eps D\U DM\U 2n^3*eps(relty)*(elty<:Complex ? 2:1) - end - for func in (det, trace) - @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) - end - if relty <: BlasFloat - for func in (expm,) - @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) - end - end - if elty <: BlasComplex - for func in (logdet, sqrtm) - @test_approx_eq_eps func(D) func(DM) n^2*eps(relty) - end - end -end - -#Test interconversion between special matrix types -N=12 -A=Diagonal([1:N]*1.0) -for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] - @test full(convert(newtype, A)) == full(A) -end - -for isupper in (true, false) - A=Bidiagonal([1:N]*1.0, [1:N-1]*1.0, isupper) - for newtype in [Bidiagonal, Tridiagonal, Triangular, Matrix] - @test full(convert(newtype, A)) == full(A) - end - A=Bidiagonal([1:N]*1.0, [1:N-1]*0.0, isupper) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] - @test full(convert(newtype, A)) == full(A) - end -end - -A=SymTridiagonal([1:N]*1.0, [1:N-1]*1.0) -for newtype in [Tridiagonal, Matrix] - @test full(convert(newtype, A)) == full(A) -end - -A=Tridiagonal([1:N-1]*0.0, [1:N]*1.0, [1:N-1]*0.0) #morally Diagonal -for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] - @test full(convert(newtype, A)) == full(A) -end - -A=Triangular(full(Diagonal([1:N]*1.0))) #morally Diagonal -for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] - @test full(convert(newtype, A)) == full(A) -end - # Test gglse for elty in (Float32, Float64, Complex64, Complex128) A = convert(Array{elty, 2}, [1 1 1 1; 1 3 1 1; 1 -1 3 1; 1 1 1 3; 1 1 1 -1]) diff -Nru julia-0.3.0/test/linalg4.jl julia-0.3.0/test/linalg4.jl --- julia-0.3.0/test/linalg4.jl 1970-01-01 00:00:00.000000000 +0000 +++ julia-0.3.0/test/linalg4.jl 2014-06-03 19:33:18.000000000 +0000 @@ -0,0 +1,313 @@ +using Base.Test +import Base.LinAlg +import Base.LinAlg: BlasComplex, BlasFloat, BlasReal + +#Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences +function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) + n = size(a, 1) + @test n==size(b,1) && size(a,2)==size(b,2) + error==nothing && (error=n^3*(eps(S)+eps(T))) + for i=1:n + ev1, ev2 = a[:,i], b[:,i] + deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) + if !isnan(deviation) + @test_approx_eq_eps deviation 0.0 error + end + end +end + +############################## +# Tests for special matrices # +############################## + +n=12 #Size of matrix problem to test + +#Triangular matrices +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) + A = convert(Matrix{elty}, randn(n, n)) + b = convert(Matrix{elty}, randn(n, 2)) + if elty <: Complex + A += im*convert(Matrix{elty}, randn(n, n)) + b += im*convert(Matrix{elty}, randn(n, 2)) + end + + for M in (triu(A), tril(A)) + TM = Triangular(M) + + ##Idempotent tests #XXX - not implemented + #for func in (conj, transpose, ctranspose) + # @test full(func(func(TM))) == M + #end + + #Linear solver + x = M \ b + tx = TM \ b + condM = elty <:BlasFloat ? cond(TM, Inf) : convert(relty, cond(complex128(M), Inf)) + @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) + if elty <: BlasFloat #test naivesub! against LAPACK + tx = [LinAlg.naivesub!(TM, b[:,1]) LinAlg.naivesub!(TM, b[:,2])] + @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) + end + + #Eigensystems + vals1, vecs1 = eig(complex128(M)) + vals2, vecs2 = eig(TM) + res1=norm(complex128(vecs1*diagm(vals1)*inv(vecs1) - M)) + res2=norm(complex128(vecs2*diagm(vals2)*inv(vecs2) - full(TM))) + @test_approx_eq_eps res1 res2 res1+res2 + + if elty <:BlasFloat + #Condition number tests - can be VERY approximate + for p in [1.0, Inf] + @test_approx_eq_eps cond(TM, p) cond(M, p) (cond(TM,p)+cond(M,p)) + end + end + + #Binary operations + B = convert(Matrix{elty}, randn(n, n)) + for M2 in (triu(B), tril(B)) + TM2 = Triangular(M2) + for op in (*,) #+, - not implemented + @test_approx_eq full(op(TM, TM2)) op(M, M2) + @test_approx_eq full(op(TM, M2)) op(M, M2) + @test_approx_eq full(op(M, TM2)) op(M, M2) + end + end + end +end + +#Tridiagonal matrices +for relty in (Float32, Float64), elty in (relty, Complex{relty}) + a = convert(Vector{elty}, randn(n-1)) + b = convert(Vector{elty}, randn(n)) + c = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Vector{elty}, randn(n)) + c += im*convert(Vector{elty}, randn(n-1)) + end + + A=Tridiagonal(a, b, c) + fA=(elty<:Complex?complex128:float64)(full(A)) + + #Simple unary functions + for func in (det, inv) + @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) + end + + #Binary operations + a = convert(Vector{elty}, randn(n-1)) + b = convert(Vector{elty}, randn(n)) + c = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Vector{elty}, randn(n)) + c += im*convert(Vector{elty}, randn(n-1)) + end + + B=Tridiagonal(a, b, c) + fB=(elty<:Complex?complex128:float64)(full(B)) + + for op in (+, -, *) + @test_approx_eq full(op(A, B)) op(fA, fB) + end +end + +#SymTridiagonal (symmetric tridiagonal) matrices +for relty in (Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work + a = convert(Vector{elty}, randn(n)) + b = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n)) + b += im*convert(Vector{elty}, randn(n-1)) + end + + A=SymTridiagonal(a, b) + fA=(elty<:Complex?complex128:float64)(full(A)) + + #Idempotent tests + for func in (conj, transpose, ctranspose) + @test func(func(A)) == A + end + + #Simple unary functions + for func in (det, inv) + @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) + end + + #Eigensystems + zero, infinity = convert(elty, 0), convert(elty, Inf) + #This tests eigenvalue and eigenvector computations using stebz! and stein! + w, iblock, isplit = LinAlg.LAPACK.stebz!('V','B',-infinity,infinity,0,0,zero,a,b) + evecs = LinAlg.LAPACK.stein!(a,b,w) + + (e, v)=eig(SymTridiagonal(a,b)) + @test_approx_eq e w + test_approx_eq_vecs(v, evecs) + + #stein! call using iblock and isplit + w, iblock, isplit = LinAlg.LAPACK.stebz!('V','B',-infinity,infinity,0,0,zero,a,b) + evecs = LinAlg.LAPACK.stein!(a,b,w,iblock,isplit) + test_approx_eq_vecs(v, evecs) + + #Binary operations + a = convert(Vector{elty}, randn(n)) + b = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Vector{elty}, randn(n)) + end + + B=SymTridiagonal(a, b) + fB=(elty<:Complex?complex128:float64)(full(B)) + + for op in (+, -, *) + @test_approx_eq full(op(A, B)) op(fA, fB) + end +end + +#Bidiagonal matrices +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) + dv = convert(Vector{elty}, randn(n)) + ev = convert(Vector{elty}, randn(n-1)) + b = convert(Matrix{elty}, randn(n, 2)) + if (elty <: Complex) + dv += im*convert(Vector{elty}, randn(n)) + ev += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Matrix{elty}, randn(n, 2)) + end + + for isupper in (true, false) #Test upper and lower bidiagonal matrices + T = Bidiagonal(dv, ev, isupper) + + @test size(T, 1) == size(T, 2) == n + @test size(T) == (n, n) + @test full(T) == diagm(dv) + diagm(ev, isupper?1:-1) + @test Bidiagonal(full(T), isupper) == T + z = zeros(elty, n) + + #Idempotent tests + for func in (conj, transpose, ctranspose) + @test func(func(T)) == T + end + + #Linear solver + Tfull = full(T) + condT = cond(complex128(Tfull)) + x = T \ b + tx = Tfull \ b + @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) + + #Eigensystems + d1, v1 = eig(T) + d2, v2 = eig((elty<:Complex?complex128:float64)(Tfull)) + @test_approx_eq isupper?d1:reverse(d1) d2 + if elty <: Real + test_approx_eq_vecs(v1, isupper?v2:v2[:,n:-1:1]) + end + + #Singular systems + if (elty <: BlasReal) + @test_approx_eq svdvals(Tfull) svdvals(T) + u1, d1, v1 = svd(Tfull) + u2, d2, v2 = svd(T) + @test_approx_eq d1 d2 + if elty <: Real + test_approx_eq_vecs(u1, u2) + test_approx_eq_vecs(v1, v2) + end + @test_approx_eq_eps 0 vecnorm(u2*diagm(d2)*v2'-Tfull) n*max(n^2*eps(relty), vecnorm(u1*diagm(d1)*v1'-Tfull)) + end + + #Binary operations + for isupper2 in (true, false) + dv = convert(Vector{elty}, randn(n)) + ev = convert(Vector{elty}, randn(n-1)) + T2 = Bidiagonal(dv, ev, isupper2) + Tfull2 = full(T2) + for op in (+, -, *) + @test_approx_eq full(op(T, T2)) op(Tfull, Tfull2) + end + end + end +end + +#Diagonal matrices +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) + d=convert(Vector{elty}, randn(n)) + v=convert(Vector{elty}, randn(n)) + U=convert(Matrix{elty}, randn(n,n)) + if elty <: Complex + d+=im*convert(Vector{elty}, randn(n)) + v+=im*convert(Vector{elty}, randn(n)) + U+=im*convert(Matrix{elty}, randn(n,n)) + end + D = Diagonal(d) + DM = diagm(d) + + #Linear solve + @test_approx_eq_eps D*v DM*v n*eps(relty)*(elty<:Complex ? 2:1) + @test_approx_eq_eps D*U DM*U n^2*eps(relty)*(elty<:Complex ? 2:1) + if relty != BigFloat + @test_approx_eq_eps D\v DM\v 2n^2*eps(relty)*(elty<:Complex ? 2:1) + @test_approx_eq_eps D\U DM\U 2n^3*eps(relty)*(elty<:Complex ? 2:1) + end + + #Simple unary functions + for func in (det, trace) + @test_approx_eq_eps func(D) func(DM) n^2*eps(relty)*(elty<:Complex ? 2:1) + end + if relty <: BlasFloat + for func in (expm,) + @test_approx_eq_eps func(D) func(DM) n^3*eps(relty) + end + end + if elty <: BlasComplex + for func in (logdet, sqrtm) + @test_approx_eq_eps func(D) func(DM) n^2*eps(relty)*2 + end + end + #Binary operations + d = convert(Vector{elty}, randn(n)) + D2 = Diagonal(d) + DM2= diagm(d) + for op in (+, -, *) + @test_approx_eq full(op(D, D2)) op(DM, DM2) + end +end + + +#Test interconversion between special matrix types +a=[1.0:n] +A=Diagonal(a) +for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] + @test full(convert(newtype, A)) == full(A) +end + +for isupper in (true, false) + A=Bidiagonal(a, [1.0:n-1], isupper) + for newtype in [Bidiagonal, Tridiagonal, Triangular, Matrix] + @test full(convert(newtype, A)) == full(A) + end + A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] + @test full(convert(newtype, A)) == full(A) + end +end + +A=SymTridiagonal(a, [1.0:n-1]) +for newtype in [Tridiagonal, Matrix] + @test full(convert(newtype, A)) == full(A) +end + +A=Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal +for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] + @test full(convert(newtype, A)) == full(A) +end + +A=Triangular(full(Diagonal(a))) #morally Diagonal +for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] + @test full(convert(newtype, A)) == full(A) +end + + diff -Nru julia-0.3.0/test/Makefile julia-0.3.0/test/Makefile --- julia-0.3.0/test/Makefile 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/Makefile 2014-06-03 19:33:18.000000000 +0000 @@ -7,7 +7,7 @@ math functional bigint sorting statistics spawn parallel arpack file \ git pkg resolve suitesparse complex version pollfd mpfr broadcast \ socket floatapprox priorityqueue readdlm regex float16 combinatorics \ - sysinfo rounding ranges mod2pi euler show linalg1 linalg2 linalg3 lineedit \ + sysinfo rounding ranges mod2pi euler show lineedit \ replcompletions backtrace repl default: all diff -Nru julia-0.3.0/test/math.jl julia-0.3.0/test/math.jl --- julia-0.3.0/test/math.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/math.jl 2014-06-03 19:33:18.000000000 +0000 @@ -149,6 +149,7 @@ @test_approx_eq lbeta(-1/2, 3) log(16/3) # gamma, lgamma (complex argument) +@test gamma(1:25) == gamma(Float64[1:25]) @test_approx_eq gamma(1/2) sqrt(π) @test_approx_eq gamma(-1/2) -2sqrt(π) @test_approx_eq lgamma(-1/2) log(abs(gamma(-1/2))) diff -Nru julia-0.3.0/test/readdlm.jl julia-0.3.0/test/readdlm.jl --- julia-0.3.0/test/readdlm.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/readdlm.jl 2014-06-03 19:33:18.000000000 +0000 @@ -102,6 +102,30 @@ @test isequaldlm(readcsv(IOBuffer("1,2,\"#3\"\n4,5,6")), [1. 2. "#3";4. 5. 6.], Any) @test isequaldlm(readcsv(IOBuffer("1,2,3\n #with leading whitespace\n4,5,6")), [1. 2. 3.;" " "" "";4. 5. 6.], Any) +# test skipstart +let x = ["a" "b" "c"; "d" "e" "f"; "g" "h" "i"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9], io = IOBuffer() + writedlm(io, x, quotes=false) + seek(io, 0) + (data, hdr) = readdlm(io, header=true, skipstart=3) + @test data == [1 2 3; 4 5 6; 7 8 9] + @test hdr == ["A" "B" "C"] +end +let x = ["a" "b" "\nc"; "d" "\ne" "f"; "g" "h" "i\n"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9] + io = IOBuffer() + writedlm(io, x, quotes=true) + seek(io, 0) + (data, hdr) = readdlm(io, header=true, skipstart=6) + @test data == [1 2 3; 4 5 6; 7 8 9] + @test hdr == ["A" "B" "C"] + + io = IOBuffer() + writedlm(io, x, quotes=false) + seek(io, 0) + (data, hdr) = readdlm(io, header=true, skipstart=6) + @test data == [1 2 3; 4 5 6; 7 8 9] + @test hdr == ["A" "B" "C"] +end + # source: http://www.i18nguy.com/unicode/unicode-example-utf8.zip let i18n_data = ["Origin (English)", "Name (English)", "Origin (Native)", "Name (Native)", "Australia", "Nicole Kidman", "Australia", "Nicole Kidman", @@ -161,4 +185,9 @@ i18n_buff = PipeBuffer() writedlm(i18n_buff, i18n_arr, ',') @test i18n_arr == readcsv(i18n_buff) + + hdr = i18n_arr[1, :] + data = i18n_arr[2:end, :] + writedlm(i18n_buff, i18n_arr, ',') + @test (data, hdr) == readcsv(i18n_buff, header=true) end diff -Nru julia-0.3.0/test/reduce.jl julia-0.3.0/test/reduce.jl --- julia-0.3.0/test/reduce.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/reduce.jl 2014-06-03 19:33:18.000000000 +0000 @@ -1,28 +1,60 @@ -## foldl & foldr -# folds -- reduce.jl -@test foldl(-,[1:5]) == -13 -@test foldl(-,10,[1:5]) == foldl(-,[10,1:5]) +# fold(l|r) & mapfold(l|r) +@test foldl(-, 1:5) == -13 +@test foldl(-, 10, 1:5) == -5 -@test foldr(-,[1:5]) == 3 -@test foldr(-,10,[1:5]) == foldr(-,[1:5,10]) +@test Base.mapfoldl(abs2, -, 2:5) == -46 +@test Base.mapfoldl(abs2, -, 10, 2:5) == -44 -# reduce -- reduce.jl +@test foldr(-, 1:5) == 3 +@test foldr(-, 10, 1:5) == -7 + +@test Base.mapfoldr(abs2, -, 2:5) == -14 +@test Base.mapfoldr(abs2, -, 10, 2:5) == -4 + +# reduce & mapreduce @test reduce((x,y)->"($x+$y)", [9:11]) == "((9+10)+11)" @test reduce(max, [8 6 7 5 3 0 9]) == 9 @test reduce(+, 1000, [1:5]) == (1000 + 1 + 2 + 3 + 4 + 5) -# mapreduce -- reduce.jl @test mapreduce(-, +, [-10 -9 -3]) == ((10 + 9) + 3) @test mapreduce((x)->x[1:3], (x,y)->"($x+$y)", ["abcd", "efgh", "01234"]) == "((abc+efg)+012)" +# sum -z = zeros(2,2,2,2) -for i=1:16 - z[i] = i -end +@test sum(Int8[]) === 0 +@test sum(Int[]) === 0 +@test sum(Float64[]) === 0.0 + +@test sum(int8(3)) === int8(3) +@test sum(3) === 3 +@test sum(3.0) === 3.0 + +@test sum([int8(3)]) === 3 +@test sum([3]) === 3 +@test sum([3.0]) === 3.0 + +z = reshape(1:16, (2,2,2,2)) +fz = float(z) +@test sum(z) === 136 +@test sum(fz) === 136.0 + +@test_throws ErrorException sum(sin, Int[]) +@test sum(sin, 3) == sin(3.0) +@test sum(sin, [3]) == sin(3.0) +@test sum(sin, z) == sum(sin, fz) == sum(sin(fz)) -@test sum(z) == sum(z,(1,2,3,4))[1] == 136 +z = [-4, -3, 2, 5] +fz = float(z) +@test Base.sumabs(Float64[]) === 0.0 +@test Base.sumabs([int8(-2)]) === 2 +@test Base.sumabs(z) === 14 +@test Base.sumabs(fz) === 14.0 + +@test Base.sumabs2(Float64[]) === 0.0 +@test Base.sumabs2([int8(-2)]) === 4 +@test Base.sumabs2(z) === 54 +@test Base.sumabs2(fz) === 54.0 # check variants of summation for type-stability and other issues (#6069) sum2(itr) = invoke(sum, (Any,), itr) @@ -48,48 +80,36 @@ end @test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7])) -prod2(itr) = invoke(prod, (Any,), itr) -@test prod(Int[]) == prod2(Int[]) == 1 -@test prod(Int[7]) == prod2(Int[7]) == 7 -@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1 7])) == typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7])) - -v = cell(2,2,1,1) -v[1,1,1,1] = 28.0 -v[1,2,1,1] = 36.0 -v[2,1,1,1] = 32.0 -v[2,2,1,1] = 40.0 +@test sum_kbn([1,1e100,1,-1e100]) == 2 -@test isequal(v,sum(z,(3,4))) +# prod -@test sum_kbn([1,1e100,1,-1e100]) == 2 +prod(Int[]) === 0 +prod(Int8[]) === 0 +prod(Float64[]) === 0.0 -z = rand(10^6) -let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5]) - @test (es - sum(z)) < es * 1e-13 - cs = cumsum(z) - @test (es - cs[end]) < es * 1e-13 - @test (es2 - cs[10^5]) < es2 * 1e-13 -end +prod([3]) === 0 +prod([int8(3)]) === 0 +prod([3.0]) === 0.0 -@test_throws ErrorException sum(sin, Int[]) -@test Base.sumabs(Float64[]) === 0.0 -@test Base.sumabs2(Float64[]) === 0.0 +prod(z) === 120 +prod(fz) === 120.0 -@test sum(sin, [1]) == sin(1) -@test Base.sumabs([int8(-2)]) === 2 -@test Base.sumabs2([int8(-2)]) === 4 +# check type-stability +prod2(itr) = invoke(prod, (Any,), itr) +@test prod(Int[]) === prod2(Int[]) === 1 +@test prod(Int[7]) === prod2(Int[7]) === 7 +@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1, 7])) == Int +@test typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7])) == Int + +# maximum & minimum & extrema -x = -2:3 -@test sum(sin, x) == sum(sin(x)) -@test Base.sumabs(x) === 9 -@test Base.sumabs2(x) === 19 - -@test_approx_eq sum(sin, z) sum(sin(z)) -@test_approx_eq Base.sumabs(z) sum(abs(z)) -@test_approx_eq Base.sumabs2(z) sum(abs2(z)) +@test_throws ErrorException maximum(Int[]) +@test_throws ErrorException minimum(Int[]) @test maximum(5) == 5 @test minimum(5) == 5 +@test extrema(5) == (5, 5) @test maximum([4, 3, 5, 2]) == 5 @test minimum([4, 3, 5, 2]) == 2 @@ -99,29 +119,87 @@ @test isnan(minimum([NaN])) @test isequal(extrema([NaN]), (NaN, NaN)) +@test maximum([NaN, 2., 3.]) == 3. +@test minimum([NaN, 2., 3.]) == 2. +@test extrema([NaN, 2., 3.]) == (2., 3.) + @test maximum([4., 3., NaN, 5., 2.]) == 5. @test minimum([4., 3., NaN, 5., 2.]) == 2. @test extrema([4., 3., NaN, 5., 2.]) == (2., 5.) -@test extrema(1:5) == (1,5) +@test Base.maxabs(Int[]) == 0 +@test_throws ErrorException Base.minabs(Int[]) @test Base.maxabs(-2) == 2 @test Base.minabs(-2) == 2 @test Base.maxabs([1, -2, 3, -4]) == 4 @test Base.minabs([-1, 2, -3, 4]) == 1 -@test maximum(abs2, 3:7) == 49 -@test minimum(abs2, 3:7) == 9 +@test maximum(x->abs2(x), 3:7) == 49 +@test minimum(x->abs2(x), 3:7) == 9 -@test any([true false; false false], 2) == [true false]' -@test any([true false; false false], 1) == [true false] +# any & all -@test all([true true; false true], 2) == [true false]' -@test all([true false; false true], 1) == [false false] +@test any(Bool[]) == false +@test any([true]) == true +@test any([false, false]) == false +@test any([false, true]) == true +@test any([true, false]) == true +@test any([true, true]) == true +@test any([true, true, true]) == true +@test any([true, false, true]) == true +@test any([false, false, false]) == false + +@test all(Bool[]) == true +@test all([true]) == true +@test all([false, false]) == false +@test all([false, true]) == false +@test all([true, false]) == false +@test all([true, true]) == true +@test all([true, true, true]) == true +@test all([true, false, true]) == false +@test all([false, false, false]) == false + +@test any(x->x>0, Int[]) == false +@test any(x->x>0, [-3]) == false +@test any(x->x>0, [4]) == true +@test any(x->x>0, [-3, 4, 5]) == true + +@test all(x->x>0, Int[]) == true +@test all(x->x>0, [-3]) == false +@test all(x->x>0, [4]) == true +@test all(x->x>0, [-3, 4, 5]) == false + +# in + +@test in(1, Int[]) == false +@test in(1, Int[1]) == true +@test in(1, Int[2]) == false +@test in(0, 1:3) == false +@test in(1, 1:3) == true +@test in(2, 1:3) == true + +# count & countnz + +@test count(x->x>0, Int[]) == 0 +@test count(x->x>0, -3:5) == 5 + +@test countnz(Int[]) == 0 +@test countnz(Int[0]) == 0 +@test countnz(Int[1]) == 1 +@test countnz([1, 0, 2, 0, 3, 0, 4]) == 4 ## cumsum, cummin, cummax +z = rand(10^6) +let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5]) + @test (es - sum(z)) < es * 1e-13 + cs = cumsum(z) + @test (es - cs[end]) < es * 1e-13 + @test (es2 - cs[10^5]) < es2 * 1e-13 +end + @test isequal(cummin([1, 2, 5, -1, 3, -2]), [1, 1, 1, -1, -1, -2]) @test isequal(cummax([1, 2, 5, -1, 3, -2]), [1, 2, 5, 5, 5, 5]) diff -Nru julia-0.3.0/test/runtests.jl julia-0.3.0/test/runtests.jl --- julia-0.3.0/test/runtests.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/runtests.jl 2014-06-03 19:33:18.000000000 +0000 @@ -20,7 +20,7 @@ if "linalg" in tests # specifically selected case filter!(x -> x != "linalg", tests) - prepend!(tests, ["linalg1", "linalg2", "linalg3"]) + prepend!(tests, ["linalg1", "linalg2", "linalg3", "linalg4"]) end net_required_for = ["socket", "parallel"] diff -Nru julia-0.3.0/test/sorting.jl julia-0.3.0/test/sorting.jl --- julia-0.3.0/test/sorting.jl 2014-06-02 19:54:00.000000000 +0000 +++ julia-0.3.0/test/sorting.jl 2014-06-03 19:33:18.000000000 +0000 @@ -50,6 +50,10 @@ @test searchsorted(1:10, 1, by=(x -> x >= 5)) == searchsorted([1:10], 1, by=(x -> x >= 5)) @test searchsorted(1:10, 10, by=(x -> x >= 5)) == searchsorted([1:10], 10, by=(x -> x >= 5)) +@test_throws BoundsError searchsortedfirst(1:10, 5, 0, 7, Base.Order.Forward) +@test_throws BoundsError searchsortedfirst(1:10, 5, 1, 11, Base.Order.Forward) +@test_throws BoundsError searchsortedfirst(1:10, 5, 5, 1, Base.Order.Forward) + a = rand(1:10000, 1000) for alg in [InsertionSort, MergeSort]