see the discussion at https://github.com/JuliaLang/julia/issues/5857 (search for NTuple). it seems that NTuple currently is not a contiguous block of memory, but a collection of pointers to its contents. it also seems that this will change...
andrew On Saturday, 3 May 2014 18:58:09 UTC-4, Andrea Pagnani wrote: > > Strange behavior > > I define Type Polynomial according to Stefan definition and > > julia> a = tuple(rand(500,1)...); b = tuple(rand(500,1)...); > > julia> @time P1=Polynomial(a); > elapsed time: 24.441438322 seconds (614490296 bytes allocated) > > julia> @time P2=Polynomial(b); > elapsed time: 2.4303e-5 seconds (7296 bytes allocated) > > julia> @time P1 + P2; > elapsed time: 0.001270669 seconds (1037376 bytes allocated) > > > > I tried with a tuple of 1000 random elements on my macBook Pro 4 Giga and > it takes forever (and starts swapping). It seems that the compiler does not > handle very well this Type (after the first round type definition for P1, > creating P2 is super fast, and uses a reasonable amount of memory). > > Of course one can use other strategies for polynomials, but the behaviour > is somehow surprising. > > > > > > > On Friday, May 2, 2014 5:55:50 PM UTC+2, Stefan Karpinski wrote: >> >> There isn't a built-in polynomial type, but it's pretty easy to make one. >> If the type of the polynomial depends on its degree, then you can dispatch >> on it. For example, you could do something like this: >> >> immutable Polynomial{n,T<:Number} >> a::NTuple{n,T} >> end >> Polynomial(a::Number...) = Polynomial(promote(a...)) >> >> function +{d1,d2}(p1::Polynomial{d1}, p2::Polynomial{d2}) >> d1 <= d2 || return p2 + p1 >> a = ntuple(d2) do k >> k <= d1 ? p1.a[k] + p2.a[k] : p2.a[k] >> end >> Polynomial(promote(a...)) >> end >> >> julia> p1 = Polynomial(1,2,3) >> Polynomial{3,Int64}((1,2,3)) >> >> julia> p2 = Polynomial(1,2) >> Polynomial{2,Int64}((1,2)) >> >> julia> p3 = Polynomial(1,1.5) >> Polynomial{2,Float64}((1.0,1.5)) >> >> julia> p1 + p2 >> Polynomial{3,Int64}((2,4,3)) >> >> julia> p1 + p3 >> Polynomial{3,Float64}((2.0,3.5,3.0)) >> >> julia> p3 + p1 >> Polynomial{3,Float64}((2.0,3.5,3.0)) >> >> julia> p3 + p2 >> Polynomial{2,Float64}((2.0,3.5)) >> >> >> Since the degree of the polynomial is part of the type, you can dispatch >> based on it. I'm not sure that's necessary though. You could always just >> check what the degree of a polynomial is and choose between algorithms >> explicitly. >> >> >> On Fri, May 2, 2014 at 11:25 AM, Albert Díaz <albee...@gmail.com> wrote: >> >>> We are trying to develop a program to find the different roots of a nth >>> degree polynomial using a few quite simple methods (fixed-point, bisection, >>> newton..) and we would like to know how, when given a function, to >>> recognise it's degree in order to choose between different root-finding >>> methods. >>> >>> >>> Thanks >>> >> >>