Re: [julia-users] How can Julia recognise the degree of a polynomial?

2014-05-03 Thread Andrea Pagnani
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.comjavascript:
  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 




Re: [julia-users] How can Julia recognise the degree of a polynomial?

2014-05-03 Thread Stefan Karpinski
Yes, sorry about that. Using tuples for large polynomials is not a good idea 
but it does allow the type system to reason about degrees. I'm not sure how 
useful having the compiler reason about degrees really is though.

 On May 3, 2014, at 8:15 PM, andrew cooke and...@acooke.org wrote:
 
 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 
 


Re: [julia-users] How can Julia recognise the degree of a polynomial?

2014-05-02 Thread Stefan Karpinski
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 albeert@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



Re: [julia-users] How can Julia recognise the degree of a polynomial?

2014-05-02 Thread Jameson Nash
See also loladiro's Polynomials.jl package, which would accept pull
requests for alternative root finding methods to be added (there's an open
issue on my original fork to add a realroots function, at minimum)

On Friday, May 2, 2014, Stefan Karpinski ste...@karpinski.org 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 
 albeert@gmail.comjavascript:_e(%7B%7D,'cvml','albeert@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