Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Sheehan Olver]

But I agree large degree polynomials are rarely needed...so this package looks 
really great! 

Sent from my iPhone

> On 30 Apr 2016, at 22:48, Mosè Giordano  wrote:
> 
> 2016-04-30 14:22 GMT+02:00 Sheehan Olver:
>> 
>> Hi Mose,
>> 
>> I think AMVW is specifically designed for large degree polynomials: it's an 
>> O(n^2) algorithm that works via orthogonal transformations on the companion 
>> matrix
>> 
>> How do the timings compare on, say, a degree 200 polynomial?
> 
> Just replace 101 with 201 in my previous script.  For 200th-degree 
> polynomial, PolynomialRoots.jl is still a bit faster, AMVW becomes faster at 
> even larger  degrees:
> 
> INFO: 400th-order polynomial
> 
> PR maximum error:   3.1349840509961133
> AMVW maximum error: 8.091848829362318e30
> 
> PR benchmark:
>  Benchmark Results 
>  Time per evaluation: 133.95 ms [125.45 ms, 142.45 ms]
> Proportion of time in GC: 0.07% [0.00%, 0.27%]
> Memory allocated: 1.20 mb
>Number of allocations: 38837 allocations
>Number of samples: 75
>Number of evaluations: 75
>  Time spent benchmarking: 10.27 s
> 
> AMVW:
>  Benchmark Results 
>  Time per evaluation: 52.97 ms [52.36 ms, 53.57 ms]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 33.64 kb
>Number of allocations: 9 allocations
>Number of samples: 100
>Number of evaluations: 100
>  Time spent benchmarking: 5.45 s
> 
> I wouldn't say that speed is of any utility if one gets such inaccurate 
> results.  Hpwever, I feel that for large degree polynomials limited precision 
> is set to be always a problem, either in root finding or polynomial 
> evaluation.
> 
> Bye,
> Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Sheehan Olver]

I think AMVW is backward stable: so it's calculating the exact roots of a near 
by polynomial

But even when roots are inaccurate, certain functions of the roots are 
accurate: e.g.

f(lambda_1) + ... + f(lambda_n)

Where f is a polynomial.

Sent from my iPhone

> On 30 Apr 2016, at 22:48, Mosè Giordano  wrote:
> 
> 2016-04-30 14:22 GMT+02:00 Sheehan Olver:
>> 
>> Hi Mose,
>> 
>> I think AMVW is specifically designed for large degree polynomials: it's an 
>> O(n^2) algorithm that works via orthogonal transformations on the companion 
>> matrix
>> 
>> How do the timings compare on, say, a degree 200 polynomial?
> 
> Just replace 101 with 201 in my previous script.  For 200th-degree 
> polynomial, PolynomialRoots.jl is still a bit faster, AMVW becomes faster at 
> even larger  degrees:
> 
> INFO: 400th-order polynomial
> 
> PR maximum error:   3.1349840509961133
> AMVW maximum error: 8.091848829362318e30
> 
> PR benchmark:
>  Benchmark Results 
>  Time per evaluation: 133.95 ms [125.45 ms, 142.45 ms]
> Proportion of time in GC: 0.07% [0.00%, 0.27%]
> Memory allocated: 1.20 mb
>Number of allocations: 38837 allocations
>Number of samples: 75
>Number of evaluations: 75
>  Time spent benchmarking: 10.27 s
> 
> AMVW:
>  Benchmark Results 
>  Time per evaluation: 52.97 ms [52.36 ms, 53.57 ms]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 33.64 kb
>Number of allocations: 9 allocations
>Number of samples: 100
>Number of evaluations: 100
>  Time spent benchmarking: 5.45 s
> 
> I wouldn't say that speed is of any utility if one gets such inaccurate 
> results.  Hpwever, I feel that for large degree polynomials limited precision 
> is set to be always a problem, either in root finding or polynomial 
> evaluation.
> 
> Bye,
> Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

2016-04-30 14:22 GMT+02:00 Sheehan Olver:
>
> Hi Mose,
>
> I think AMVW is specifically designed for large degree polynomials: it's 
> an O(n^2) algorithm that works via orthogonal transformations on the 
> companion matrix
>
> How do the timings compare on, say, a degree 200 polynomial?
>

Just replace 101 with 201 in my previous script.  For 200th-degree 
polynomial, PolynomialRoots.jl is still a bit faster, AMVW becomes faster 
at even larger  degrees:

INFO: 400th-order polynomial

PR maximum error:   3.1349840509961133
AMVW maximum error: 8.091848829362318e30

PR benchmark:
 Benchmark Results 
 Time per evaluation: 133.95 ms [125.45 ms, 142.45 ms]
Proportion of time in GC: 0.07% [0.00%, 0.27%]
Memory allocated: 1.20 mb
   Number of allocations: 38837 allocations
   Number of samples: 75
   Number of evaluations: 75
 Time spent benchmarking: 10.27 s

AMVW:
 Benchmark Results 
 Time per evaluation: 52.97 ms [52.36 ms, 53.57 ms]
Proportion of time in GC: 0.00% [0.00%, 0.00%]
Memory allocated: 33.64 kb
   Number of allocations: 9 allocations
   Number of samples: 100
   Number of evaluations: 100
 Time spent benchmarking: 5.45 s

I wouldn't say that speed is of any utility if one gets such inaccurate 
results.  Hpwever, I feel that for large degree polynomials limited 
precision is set to be always a problem, either in root finding or 
polynomial evaluation.

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Sheehan Olver]

Hi Mose,

I think AMVW is specifically designed for large degree polynomials: it's an 
O(n^2) algorithm that works via orthogonal transformations on the companion 
matrix

How do the timings compare on, say, a degree 200 polynomial?

Sent from my iPhone

> On 30 Apr 2016, at 20:44, Mosè Giordano  wrote:
> 
> Hi Sheehan,
> 
> 2016-04-30 8:12 GMT+02:00 Sheehan Olver:
>> 
>> How does this compare to AMVW.jl?
>> 
>> https://github.com/andreasnoack/AMVW.jl
> 
> I compared both programs with this script (requires development version of 
> PolynomialRoots):
> 
> using Benchmarks, PolynomialRoots, AMVW
> 
> for polynomial in ([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 
> 9im), 1.],
>rand(Complex{Float64}, 101))
> info("$(length(polynomial)-1)th-order polynomial")
> PRroots  = PolynomialRoots.roots(polynomial)
> AMVWroots = AMVW.rootsAMVW(polynomial)
> println("\nPR maximum error:   ",
> maximum(abs(PolynomialRoots.evalpoly(PRroots, polynomial
> println("AMVW maximum error: ",
> maximum(abs(PolynomialRoots.evalpoly(AMVWroots, polynomial
> println("\nPR benchmark:\n", @benchmark PolynomialRoots.roots(polynomial))
> println("AMVW:\n", @benchmark AMVW.rootsAMVW(polynomial))
> end
> 
>  This is the result:
> 
> INFO: 5th-order polynomial
> 
> PR maximum error:   8.526512829121202e-14
> AMVW maximum error: 3.0505379675102492e-12
> 
> PR benchmark:
>  Benchmark Results 
>  Time per evaluation: 2.82 μs [2.74 μs, 2.91 μs]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 656.00 bytes
>Number of allocations: 9 allocations
>Number of samples: 2001
>Number of evaluations: 12501
>  R² of OLS model: 0.952
>  Time spent benchmarking: 0.78 s
> 
> AMVW:
>  Benchmark Results 
>  Time per evaluation: 26.72 μs [24.77 μs, 28.67 μs]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 1.03 kb
>Number of allocations: 9 allocations
>Number of samples: 100
>Number of evaluations: 100
>  Time spent benchmarking: 0.12 s
> 
> INFO: 100th-order polynomial
> 
> PR maximum error:   2.2724270820617676e-7
> AMVW maximum error: 1.6100194395356297e-6
> 
> PR benchmark:
>  Benchmark Results 
>  Time per evaluation: 904.96 μs [887.77 μs, 922.15 μs]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 18.45 kb
>Number of allocations: 424 allocations
>Number of samples: 100
>Number of evaluations: 100
>  Time spent benchmarking: 0.19 s
> 
> AMVW:
>  Benchmark Results 
>  Time per evaluation: 4.03 ms [3.97 ms, 4.09 ms]
> Proportion of time in GC: 0.00% [0.00%, 0.00%]
> Memory allocated: 9.41 kb
>Number of allocations: 9 allocations
>Number of samples: 100
>Number of evaluations: 100
>  Time spent benchmarking: 0.51 s
> 
> So PolynomialRoots.jl is both faster and more precise than AMVW.jl.  Note 
> that this package employs an algorithm based on finding eigenvalues of 
> companion matrix, like Polynomials.jl.  Note that both algorithms often fail 
> to converge to the actual roots for large degree polynomials, but next 
> version of PolynomialRoots will feature support for multiple precision 
> calculations.  I'll hopefully release it in the next few days.
> 
> Bye,
> Mosè

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi Sheehan,

2016-04-30 8:12 GMT+02:00 Sheehan Olver:
>
> How does this compare to AMVW.jl?
>
> https://github.com/andreasnoack/AMVW.jl
>

I compared both programs with this script (requires development version of 
PolynomialRoots):

using Benchmarks, PolynomialRoots, AMVW

for polynomial in ([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 
9im), 1.],
   rand(Complex{Float64}, 101))
info("$(length(polynomial)-1)th-order polynomial")
PRroots  = PolynomialRoots.roots(polynomial)
AMVWroots = AMVW.rootsAMVW(polynomial)
println("\nPR maximum error:   ",
maximum(abs(PolynomialRoots.evalpoly(PRroots, polynomial
println("AMVW maximum error: ",
maximum(abs(PolynomialRoots.evalpoly(AMVWroots, polynomial
println("\nPR benchmark:\n", @benchmark PolynomialRoots.roots(polynomial
))
println("AMVW:\n", @benchmark AMVW.rootsAMVW(polynomial))
end

 This is the result:

INFO: 5th-order polynomial

PR maximum error:   8.526512829121202e-14
AMVW maximum error: 3.0505379675102492e-12

PR benchmark:
 Benchmark Results 
 Time per evaluation: 2.82 μs [2.74 μs, 2.91 μs]
Proportion of time in GC: 0.00% [0.00%, 0.00%]
Memory allocated: 656.00 bytes
   Number of allocations: 9 allocations
   Number of samples: 2001
   Number of evaluations: 12501
 R² of OLS model: 0.952
 Time spent benchmarking: 0.78 s

AMVW:
 Benchmark Results 
 Time per evaluation: 26.72 μs [24.77 μs, 28.67 μs]
Proportion of time in GC: 0.00% [0.00%, 0.00%]
Memory allocated: 1.03 kb
   Number of allocations: 9 allocations
   Number of samples: 100
   Number of evaluations: 100
 Time spent benchmarking: 0.12 s

INFO: 100th-order polynomial

PR maximum error:   2.2724270820617676e-7
AMVW maximum error: 1.6100194395356297e-6

PR benchmark:
 Benchmark Results 
 Time per evaluation: 904.96 μs [887.77 μs, 922.15 μs]
Proportion of time in GC: 0.00% [0.00%, 0.00%]
Memory allocated: 18.45 kb
   Number of allocations: 424 allocations
   Number of samples: 100
   Number of evaluations: 100
 Time spent benchmarking: 0.19 s

AMVW:
 Benchmark Results 
 Time per evaluation: 4.03 ms [3.97 ms, 4.09 ms]
Proportion of time in GC: 0.00% [0.00%, 0.00%]
Memory allocated: 9.41 kb
   Number of allocations: 9 allocations
   Number of samples: 100
   Number of evaluations: 100
 Time spent benchmarking: 0.51 s

So PolynomialRoots.jl is both faster and more precise than AMVW.jl.  Note 
that this package employs an algorithm based on finding eigenvalues of 
companion matrix, like Polynomials.jl.  Note that both algorithms often 
fail to converge to the actual roots for large degree polynomials, but next 
version of PolynomialRoots will feature support for multiple precision 
calculations.  I'll hopefully release it in the next few days.

Bye,
Mosè

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Sheehan Olver]

How does this compare to AMVW.jl?

https://github.com/andreasnoack/AMVW.jl

On Wednesday, April 27, 2016 at 9:15:35 PM UTC+10, Mosè Giordano wrote:
>
> Dear all,
>
> I am happy to announce the first version of CmplxRoots.jl 
> , a fast root finder of real 
> and complex polynomials.
>
> This is a Julia implementation of the algorithm described in 
>
>- J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver 
>and Its Further Optimization for Binary Microlenses", arXiv:1203.1034 
>
>
> I first wrote the package as a wrapper around the Fortran implementation 
> available at http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ but then 
> I rewrote it entirely in Julia, finding that it's even faster than 
> ccall'ing functions in the shared object.
>
> CmplxRoots.jl is a registered Julia package, so you can install it with 
> the package manager:
>
> Pkg.add("CmplxRoots")
>
> Two functions are exposed to the users:
>
> roots(polynomial[, roots, polish=true])
> roots5(polynomial[, roots])
>
> The first function can be used to solve polynomials of any order, the 
> second one is specialized and optimized for (and works only for) 
> polynomials of fifth-order, but I found it to be a bit slower than roots()
> .
>
> The only mandatory argument is the list of coefficients of the polynomial, 
> in ascending order.  So, if you want to find all the roots of polynomial
>
> (x - i)*(x - 2)*(x - 3*i)*(x - 4)*(x - 5*i) =
>   x^5 - (9*i + 6)*x^4 + (54*i - 15)*x^3 + (138 - 57*i)*x^2 - (184 + 90*i)*x + 
> 120*i
>
>
> you can use
>
> julia> roots([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im), 
> 1])
> 5-element Array{Complex{Float64},1}:
>  -1.55431e-15+5.0im
>  4.0+1.77636e-16im
>   1.55028e-15+3.0im
>  -1.04196e-16+1.0im
>  2.0-2.00595e-16im
>
> julia> roots5([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im
> ), 1])
> 5-element Array{Complex{Float64},1}:
>   5.92119e-16+5.0im
>   4.0-1.4803e-16im
>  2.0+1.88202e-16im
>  -1.88738e-15+3.0im
>   2.10942e-15+1.0im
>
> Optional roots argument is the array of starting values that can be used 
> to find the final roots (it won't be changed in-place).  Note that even if 
> the roots are all real, the output array is always of type Complex{Float64}.
>
> Also Polynomials.jl package has a root finding function, but according to 
> some random tests performed with @benchmark (from Benchmarks.jl package) 
> it's roughly 20 times slower than roots() function provided by 
> CmplxRoots.jl.  Of course do not expect results with machine precision (you 
> can have catastrophic cancellation even for simple quadratic polynomials), 
> but some random tests indicate that CmplxRoots.jl is slightly more precise 
> than Polynomials.jl (just compute abs2() of polynomials at the root).
>
> The package CmplxRoots.jl is licensed under the GNU Lesser General Public 
> License version 3 or any later version, as well as under a "customary 
> scientific license", which implies that if this code was important in the 
> scientific process or for the results of your scientific work, you are 
> asked for the appropriate citation of the paper Skowron & Gould 2012 (
> http://arxiv.org/abs/1203.1034).
>
> Happy root-finding,
> Mosè
>

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi David,

2016-04-30 0:17 GMT+02:00 David P. Sanders :
>
>
> El jueves, 28 de abril de 2016, 16:09:54 (UTC-4), Mosè Giordano escribió:
>>
>> Hi David,
>>
>> 2016-04-28 14:10 GMT+02:00 David P. Sanders :
>> > As far as I could see, the original library has an Apache license, so
>> > you should be able to use MIT.
>> >
>> > I believe that you need to include a copy of the original license in
>> > your package?
>>
>> Original code is dual-licensed, in that case the author of derived
>> work can choose one of the two licenses or both.[1]  I went for LGPL,
>> that is weakly protective
>
>
> What is being protected, your own work?

And of the authors of original library.

>> but perfectly compatible with MIT "Expat"
>
>
> No, it is only half compatible with MIT, i.e. only in one direction: you can
> take code
> from my MIT-licensed package (e.g. ValidatedNumerics) and use it in yours,
> but I *cannot*
> take a piece of code from your package, modify it, and re-use it in mine,
> without "infecting" my package with the LGPL code, which will automatically
> mean that
> I will have to change my license to LGPL.
>
> (This is all only my understanding. Please do correct me if I am wrong.)

No, this isn't correct.  You can take an LPGL library and use it in
your (even proprietary) program without the need to change the license
of your program.  What you can't do is to change the license of the
LGPLed library itself into a more permissive one, but the library
doesn't affect the license of programs making use of it.  This is the
main difference with GLP.  Remember that the first "L" stands for
"Lesser" ;-)  For this reason, an LGPL library can be used in a MIT
program.

>> license, which is the most common in Julia ecosystem, so that most of
>> Julia package can use this package.
>
>
> Again, they can only *use* it.
> I cannot even *look* at your code (which I would very much like to do),
> since if I find something useful, I can't reuse it in my (MIT) package.

You can have the core of your program with MIT license and still
incorporate LGPL code, provided that that part is released under the
terms of LGPL.  Julia itself has some files with separate licenses,
including LGPL:
https://github.com/JuliaLang/julia/blob/master/LICENSE.md  This
doesn't affect the rest of the program.

>> For the time being I don't plan
>> to change the license.
>
>
> I would urge you to reconsider that. You could, for example, license it
> with the same dual license as in the original package -- otherwise, you
> are actually restricting your users more than the original package does
> (for the reasons I discussed above).
>
> The end result otherwise will be that someone
> who needs to use the functionality will end up going back to the original
> Fortran code
> and rewriting what you have already done (which I certainly do not want to
> do),
> or using less-good code from somewhere else.

Ok, I'll consider using the same dual-license as in the original
library (no guarantee), but please note my above comments, they should
alleviate your worries ;-)

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Steven G. Johnson]

On Friday, April 29, 2016 at 3:58:56 AM UTC-4, 2n wrote:
>
> So translating header files for bindings is considered a derived work also?
>

It depends on what's in the header file.   Of course, if you are calling 
someone else's library, then you generally are a derived work anyway.
 

> And also the GNU guys and the porting of unix sofware to linux, that can't 
> be considered derived work.
>

GNU re-implemented the Unix userland from scratch.  They didn't look at 
anyone's code, only the specifications.
 

> I guess it comes down to the definition of the word "expressive" like you 
> say, but programming languages are so precise, I wouldn't have thought 
> there is much room for expressiveness.
>
 
Try looking at two people's programs that do "the same thing" for anything 
moderately complicated.   Programming is quite a creative exercise.

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[David P. Sanders]

El jueves, 28 de abril de 2016, 16:09:54 (UTC-4), Mosè Giordano escribió:
>
> Hi David, 
>
> 2016-04-28 14:10 GMT+02:00 David P. Sanders  >: 
> > As far as I could see, the original library has an Apache license, so 
> you should be able to use MIT. 
> > 
> > I believe that you need to include a copy of the original license in 
> your package? 
>
> Original code is dual-licensed, in that case the author of derived 
> work can choose one of the two licenses or both.[1]  I went for LGPL, 
> that is weakly protective 


What is being protected, your own work? 

but perfectly compatible with MIT "Expat" 
>

No, it is only half compatible with MIT, i.e. only in one direction: you 
can take code
from my MIT-licensed package (e.g. ValidatedNumerics) and use it in yours, 
but I *cannot*
take a piece of code from your package, modify it, and re-use it in mine,
without "infecting" my package with the LGPL code, which will automatically 
mean that 
I will have to change my license to LGPL. 

(This is all only my understanding. Please do correct me if I am wrong.)

 

> license, which is the most common in Julia ecosystem, so that most of 
> Julia package can use this package.  


Again, they can only *use* it.
I cannot even *look* at your code (which I would very much like to do),
since if I find something useful, I can't reuse it in my (MIT) package.

 

> For the time being I don't plan 
> to change the license. 
>

I would urge you to reconsider that. You could, for example, license it 
with the same dual license as in the original package -- otherwise, you 
are actually restricting your users more than the original package does 
(for the reasons I discussed above). 

The end result otherwise will be that someone
who needs to use the functionality will end up going back to the original 
Fortran code
and rewriting what you have already done (which I certainly do not want to 
do), 
or using less-good code from somewhere else.
 

>
> > +1 for ComplexRoots.jl 


Steven makes a good point that PolynomialRoots.jl would be a better name.

 

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

2016-04-28 22:48 GMT+02:00 Mosè Giordano :
> Hi Steven,
>
> 2016-04-28 14:58 GMT+02:00 Steven G. Johnson :
>> On Thursday, April 28, 2016 at 8:46:06 AM UTC-4, Steven G. Johnson wrote:
>>>
>>> Out of curiosity, how does the performance and accuracy compare to the
>>> "roots" function in the Polynomials.jl package, which uses the standard
>>> companion-matrix approach?
>>
>>
>> Just checked: it is much faster than Polynomials.roots, which I guess is not
>> surprising for a Newton-like method compared to an eigenvalue method,
>
> Yes, as I wrote in my first mail I found my package to be ~20 times
> faster than Polynomials.jl and usually more precise (**very** roughly,
> evaluating abs2(evaluation of polynomial at root) result in ~10^-29
> for Polynomials.jl and ~10^-30 for CmplxRoots.jl).  Probably I
> wouldn't have published the package if it had been both slower and less 
> precise.
>
>> although perhaps the latter is more reliable in some cases?
>>
>> I'm not sure what the state of the art here is, though.
>
> Most of the methods I found are based on finding eigenvalues of
> companion matrix.  A somewhat popular method is Jenkins–Traub
> algorithm (https://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_algorithm),
> but I didn't find a FLOSS implementation of that algorithm in any
> language, to start with.
>
> In my test suite I added a couple of tests of random polynomials, to
> catch major flaws of root finding algorithm.

I hit a limitation of the algorithm used in PolynomailRoots.jl (new
name of CmplxRoots.jl, but it is hasn't been merged in METADATA.jl
yet): it gives highly inaccurate results for large-order polynomials
(above ~30).  This isn't a systematic behavior, it happens from time
to time, but the larger the degree of polynomial, the more failures
you'll get.  Fortunately, this can be mitigated by using new multiple
precision feature.  For the record, Polynomials.jl has the same
problem as PolynomialRoots.jl, with the former still performing worse
than the latter also at high orders.

Regarding interesting algorithms, Aberth–Ehrlich method
(https://en.wikipedia.org/wiki/Aberth_method) looks promising and it
has already been implemented in multiple precision (see
https://en.wikipedia.org/wiki/MPSolve and
http://numpi.dm.unipi.it/mpsolve-2.2/, this isn't FLOSS though).

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

2016-04-29 9:58 GMT+02:00 2n :
> So translating header files for bindings is considered a derived work also?
> And also the GNU guys and the porting of unix sofware to linux, that can't
> be considered derived work. I guess it comes down to the definition of the
> word "expressive" like you say, but programming languages are so precise, I
> wouldn't have thought there is much room for expressiveness.

It depends on what you translate.  If it's really needed to
communicate with the library you link to, that may not be considered a
derived work.  If no creativity is involved you may say that there is
no copyright on that part of the code, but I guess this strongly
depends also on local jurisdictions.

LPGL addresses this question for object code, see section 3:

--8<---cut here---start->8---
  3. Object Code Incorporating Material from Library Header Files.

  The object code form of an Application may incorporate material from
a header file that is part of the Library.  You may convey such object
code under terms of your choice, provided that, if the incorporated
material is not limited to numerical parameters, data structure
layouts and accessors, or small macros, inline functions and templates
(ten or fewer lines in length), you do both of the following:

   a) Give prominent notice with each copy of the object code that the
   Library is used in it and that the Library and its use are
   covered by this License.

   b) Accompany the object code with a copy of the GNU GPL and this license
   document.
--8<---cut here---end--->8---

So, if you limit to something really essential, no copyright is claimed on it.

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[2n]

So translating header files for bindings is considered a derived work also? 
And also the GNU guys and the porting of unix sofware to linux, that can't 
be considered derived work. I guess it comes down to the definition of the 
word "expressive" like you say, but programming languages are so precise, I 
wouldn't have thought there is much room for expressiveness.

On Friday, 29 April 2016 10:22:33 UTC+8, Steven G. Johnson wrote:
>
> On Thursday, April 28, 2016 at 9:21:03 PM UTC-4, 2n wrote:
>>
>> If you re-write or translate something doesn't it necessarily mean that 
>> you are no longer using licensed code? that you have full copyrights to the 
>> new code? 
>>
>
> No, typically a translation would be considered a "derived work" of the 
> original work.   For a derived work, the copyright is jointly held by both 
> you and the original author, i.e. you need permission from both to 
> distribute (which, in the free/open-source world, you typically get via a 
> license, in this case the LGPL).
>
> A re-write is a tricker issue.   A "clean room" rewrite, in which the 
> algorithm is re-implemented by someone who hasn't looked at the original 
> code, but has only read the mathematical description of the algorithm, is 
> typically only copyrighted by the new author.   However, if you've 
> rewritten after studying the original code, it may be difficult to avoid 
> copying "expressive" elements of the original code that are copyrighted.
>
> (Ultimately, there is no deterministic way to know whether something is a 
> derived work of something else -- it can only be determined on a 
> case-by-case basis in a court of law.   In practice, people usually assume 
> that if you've studied someone's code and implement something similar then 
> it is probably a derived work.) 
>

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Steven G. Johnson]

On Thursday, April 28, 2016 at 9:21:03 PM UTC-4, 2n wrote:
>
> If you re-write or translate something doesn't it necessarily mean that 
> you are no longer using licensed code? that you have full copyrights to the 
> new code? 
>

No, typically a translation would be considered a "derived work" of the 
original work.   For a derived work, the copyright is jointly held by both 
you and the original author, i.e. you need permission from both to 
distribute (which, in the free/open-source world, you typically get via a 
license, in this case the LGPL).

A re-write is a tricker issue.   A "clean room" rewrite, in which the 
algorithm is re-implemented by someone who hasn't looked at the original 
code, but has only read the mathematical description of the algorithm, is 
typically only copyrighted by the new author.   However, if you've 
rewritten after studying the original code, it may be difficult to avoid 
copying "expressive" elements of the original code that are copyrighted.

(Ultimately, there is no deterministic way to know whether something is a 
derived work of something else -- it can only be determined on a 
case-by-case basis in a court of law.   In practice, people usually assume 
that if you've studied someone's code and implement something similar then 
it is probably a derived work.) 

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[2n]

On Wednesday, 27 April 2016 21:10:47 UTC+8, Mosè Giordano wrote:
>
> Hi Kristoffer, 
>
> 2016-04-27 14:32 GMT+02:00 Kristoffer Carlsson  >: 
> > This looks interesting. 
>
> Thank you! 
>
> > It is in my opinion a bit unfortunate that you chose 
> > to use the GPL license for this since it means that it will not be 
> usable by 
> > most packages in Julia which are under MIT license. 
>
> Actually I used LGPL, not GPL, because the original Fortran library I 
> translated to Julia is released with that license, and, as far as I 
> know, you can you an LGPL program into a MIT one. 
>
> Bye, 
> Mosè 
>

If you re-write or translate something doesn't it necessarily mean that you 
are no longer using licensed code? that you have full copyrights to the new 
code? 

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

2016-04-28 22:34 GMT+02:00 Steven G. Johnson :
>
>
> On Thursday, April 28, 2016 at 4:09:54 PM UTC-4, Mosè Giordano wrote:
>>
>> > +1 for ComplexRoots.jl
>>
>> Ok, then I'll see how to change the name of the package in METADATA.jl
>
>
> I would suggest something like PolynomialRoots.jl; the names "ComplexRoots"
> are "CmplxRoots" are not very descriptive.

Good point, thanks!

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi Steven,

2016-04-28 14:58 GMT+02:00 Steven G. Johnson :
> On Thursday, April 28, 2016 at 8:46:06 AM UTC-4, Steven G. Johnson wrote:
>>
>> Out of curiosity, how does the performance and accuracy compare to the
>> "roots" function in the Polynomials.jl package, which uses the standard
>> companion-matrix approach?
>
>
> Just checked: it is much faster than Polynomials.roots, which I guess is not
> surprising for a Newton-like method compared to an eigenvalue method,

Yes, as I wrote in my first mail I found my package to be ~20 times
faster than Polynomials.jl and usually more precise (**very** roughly,
evaluating abs2(evaluation of polynomial at root) result in ~10^-29
for Polynomials.jl and ~10^-30 for CmplxRoots.jl).  Probably I
wouldn't have published the package if it had been both slower and less precise.

> although perhaps the latter is more reliable in some cases?
>
> I'm not sure what the state of the art here is, though.

Most of the methods I found are based on finding eigenvalues of
companion matrix.  A somewhat popular method is Jenkins–Traub
algorithm (https://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_algorithm),
but I didn't find a FLOSS implementation of that algorithm in any
language, to start with.

In my test suite I added a couple of tests of random polynomials, to
catch major flaws of root finding algorithm.

> (I find it a little bit worrying that the Skowron and Gould paper justify
> their method on account of it being "1.6 to 3 times" faster than the roots
> function in Numerical Recipes (NR), and that NR was the inspiration for
> their algorithm.  NR routines typically have deeply suboptimal performance,
> and that book is terrible starting point for developing state-of-the-art
> algorithms because NR's approaches are typically antiquated.)

I know that algorithm showed in Numerical Recipes are often far from
being the state of art (in addition, their routines aren't FLOSS), but
that series of book is fairly common in astrophysics even today, maybe
because at least the first two authors are well known astrophysicists.

I started from Skowron & Gould package because their package is FLOSS,
it's in Fortran and it's easy to be called from Julia (indeed I first
wrote a wrapper around the Fortran shared object, that took me a few
minutes), and I used it often, without major problems.

Bye,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Steven G. Johnson]

On Thursday, April 28, 2016 at 4:09:54 PM UTC-4, Mosè Giordano wrote:
>
> > +1 for ComplexRoots.jl 
>
> Ok, then I'll see how to change the name of the package in METADATA.jl 
>

I would suggest something like PolynomialRoots.jl; the names "ComplexRoots" 
are "CmplxRoots" are not very descriptive.

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi David,

2016-04-28 14:10 GMT+02:00 David P. Sanders :
> As far as I could see, the original library has an Apache license, so you 
> should be able to use MIT.
>
> I believe that you need to include a copy of the original license in your 
> package?

Original code is dual-licensed, in that case the author of derived
work can choose one of the two licenses or both.[1]  I went for LGPL,
that is weakly protective but perfectly compatible with MIT "Expat"
license, which is the most common in Julia ecosystem, so that most of
Julia package can use this package.  For the time being I don't plan
to change the license.

> +1 for ComplexRoots.jl

Ok, then I'll see how to change the name of the package in METADATA.jl

In a separate branch of current repository I'm also working on adding
support for arbitrary precision calculations, so the renamed package
will feature this new capability.

Bye,
Mosè


Note
[1] See also http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/: «The
authors release the source codes associated with the Paper under terms
of the GNU Lesser General Public License version 2 or any later
version, **or** under the Apache License, Version 2.0»

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Steven G. Johnson]

On Thursday, April 28, 2016 at 8:46:06 AM UTC-4, Steven G. Johnson wrote:
>
> Out of curiosity, how does the performance and accuracy compare to the 
> "roots" function in the Polynomials.jl package, which uses the standard 
> companion-matrix approach?
>

Just checked: it is much faster than Polynomials.roots, which I guess is 
not surprising for a Newton-like method compared to an eigenvalue method, 
although perhaps the latter is more reliable in some cases?

I'm not sure what the state of the art here is, though. 

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Steven G. Johnson]

On Wednesday, April 27, 2016 at 7:15:35 AM UTC-4, Mosè Giordano wrote:
>
> Dear all,
>
> I am happy to announce the first version of CmplxRoots.jl 
> , a fast root finder of real 
> and complex polynomials.
>
> This is a Julia implementation of the algorithm described in 
>
>- J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver 
>and Its Further Optimization for Binary Microlenses", arXiv:1203.1034 
>
>
> Out of curiosity, how does the performance and accuracy compare to the 
"roots" function in the Polynomials.jl package, which uses the standard 
companion-matrix approach?

(I find it a little bit worrying that the Skowron and Gould paper justify 
their method on account of it being "1.6 to 3 times" faster than the roots 
function in *Numerical Recipes* (NR), and that NR was the inspiration for 
their algorithm.  NR routines typically have deeply suboptimal performance, 
and that book is terrible starting point for developing state-of-the-art 
algorithms because NR's approaches are typically antiquated.) 

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[David P. Sanders]

As far as I could see, the original library has an Apache license, so you 
should be able to use MIT. 

I believe that you need to include a copy of the original license in your 
package? 

+1 for ComplexRoots.jl

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[DNF]

There is a reason, then ;)

But I see you have already changed the name from cmplx_roots to 
CmplxRoots.jl to keep up with Julia naming conventions, so going all the 
way would be even better! Who knows what odd reasons or compromises may be 
behind the original name.

Just to refer to the Julia style guide 

 :

>
>- conciseness is valued, but avoid abbreviation (indexin() 
>
>  
> rather 
>than indxin()) as it becomes difficult to remember whether and how 
>particular words are abbreviated.
>
> It doesn't say anything about packages, but I don't see why it should not 
apply. 

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi,

2016-04-28 12:57 GMT+02:00 DNF :
> Just a question: Is there some particular reason you chose to to skip the
> letters 'o' and 'e' in the name?
>
> It seems like the kind of thing that would make the package hard to find,
> and cause confusion and frustration. In fact, while I immediately saw that
> you skipped the 'e', it read it many times before realizing that the 'o' was
> also gone. Doing Pkg.add() I would have been pretty perplexed, and every
> time >> using Com it would fail.

It stems from the name of the original Fortran code the package is
based on: http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/

If also other people will suggest to change the name I can consider
this possibility.

Cheers,
Mosè

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[DNF]

Just a question: Is there some particular reason you chose to to skip the 
letters 'o' and 'e' in the name?

It seems like the kind of thing that would make the package hard to find, 
and cause confusion and frustration. In fact, while I immediately saw that 
you skipped the 'e', it read it many times before realizing that the 'o' 
was also gone. Doing Pkg.add() I would have been pretty perplexed, and 
every time >> using Com it would fail.


On Wednesday, April 27, 2016 at 1:15:35 PM UTC+2, Mosè Giordano wrote:
>
> Dear all,
>
> I am happy to announce the first version of CmplxRoots.jl 
> , a fast root finder of real 
> and complex polynomials.
>
> This is a Julia implementation of the algorithm described in 
>
>- J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver 
>and Its Further Optimization for Binary Microlenses", arXiv:1203.1034 
>
>
> I first wrote the package as a wrapper around the Fortran implementation 
> available at http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ but then 
> I rewrote it entirely in Julia, finding that it's even faster than 
> ccall'ing functions in the shared object.
>
> CmplxRoots.jl is a registered Julia package, so you can install it with 
> the package manager:
>
> Pkg.add("CmplxRoots")
>
> Two functions are exposed to the users:
>
> roots(polynomial[, roots, polish=true])
> roots5(polynomial[, roots])
>
> The first function can be used to solve polynomials of any order, the 
> second one is specialized and optimized for (and works only for) 
> polynomials of fifth-order, but I found it to be a bit slower than roots()
> .
>
> The only mandatory argument is the list of coefficients of the polynomial, 
> in ascending order.  So, if you want to find all the roots of polynomial
>
> (x - i)*(x - 2)*(x - 3*i)*(x - 4)*(x - 5*i) =
>   x^5 - (9*i + 6)*x^4 + (54*i - 15)*x^3 + (138 - 57*i)*x^2 - (184 + 90*i)*x + 
> 120*i
>
>
> you can use
>
> julia> roots([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im), 
> 1])
> 5-element Array{Complex{Float64},1}:
>  -1.55431e-15+5.0im
>  4.0+1.77636e-16im
>   1.55028e-15+3.0im
>  -1.04196e-16+1.0im
>  2.0-2.00595e-16im
>
> julia> roots5([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im
> ), 1])
> 5-element Array{Complex{Float64},1}:
>   5.92119e-16+5.0im
>   4.0-1.4803e-16im
>  2.0+1.88202e-16im
>  -1.88738e-15+3.0im
>   2.10942e-15+1.0im
>
> Optional roots argument is the array of starting values that can be used 
> to find the final roots (it won't be changed in-place).  Note that even if 
> the roots are all real, the output array is always of type Complex{Float64}.
>
> Also Polynomials.jl package has a root finding function, but according to 
> some random tests performed with @benchmark (from Benchmarks.jl package) 
> it's roughly 20 times slower than roots() function provided by 
> CmplxRoots.jl.  Of course do not expect results with machine precision (you 
> can have catastrophic cancellation even for simple quadratic polynomials), 
> but some random tests indicate that CmplxRoots.jl is slightly more precise 
> than Polynomials.jl (just compute abs2() of polynomials at the root).
>
> The package CmplxRoots.jl is licensed under the GNU Lesser General Public 
> License version 3 or any later version, as well as under a "customary 
> scientific license", which implies that if this code was important in the 
> scientific process or for the results of your scientific work, you are 
> asked for the appropriate citation of the paper Skowron & Gould 2012 (
> http://arxiv.org/abs/1203.1034).
>
> Happy root-finding,
> Mosè
>

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Avik Sengupta]

Of course, if you started with a LGPL library, then this much be LGPL. 

Just for some context, a majority of Julia packages, including the core 
language are under MIT license. So there is a certain preference for that, 
but other licenses are certainly not precluded. There is a lot of GPL/LGPL 
code within the Julia ecosystem. 

Regards
-
Avik

On Thursday, 28 April 2016 08:54:53 UTC+1, Mosè Giordano wrote:
>
> 2016-04-27 15:10 GMT+02:00 Mosè Giordano  >: 
> > Hi Kristoffer, 
> > 
> > 2016-04-27 14:32 GMT+02:00 Kristoffer Carlsson  >: 
> >> This looks interesting. 
> > 
> > Thank you! 
> > 
> >> It is in my opinion a bit unfortunate that you chose 
> >> to use the GPL license for this since it means that it will not be 
> usable by 
> >> most packages in Julia which are under MIT license. 
> > 
> > Actually I used LGPL, not GPL, because the original Fortran library I 
> > translated to Julia is released with that license, and, as far as I 
> > know, you can you an LGPL program into a MIT one. 
> ^^^ 
>
> Oops, while writing the email I skipped the verb: in place of the 
> second wrong "you" read "use".  To avoid further confusion, a MIT 
> program can link to an LPGL one, provided that certain conditions are 
> met (but no license change is required), see section 4 of LGPLv3 text. 
>
> Cheers, 
> Mosè 
>

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

2016-04-27 15:10 GMT+02:00 Mosè Giordano :
> Hi Kristoffer,
>
> 2016-04-27 14:32 GMT+02:00 Kristoffer Carlsson :
>> This looks interesting.
>
> Thank you!
>
>> It is in my opinion a bit unfortunate that you chose
>> to use the GPL license for this since it means that it will not be usable by
>> most packages in Julia which are under MIT license.
>
> Actually I used LGPL, not GPL, because the original Fortran library I
> translated to Julia is released with that license, and, as far as I
> know, you can you an LGPL program into a MIT one.
^^^

Oops, while writing the email I skipped the verb: in place of the
second wrong "you" read "use".  To avoid further confusion, a MIT
program can link to an LPGL one, provided that certain conditions are
met (but no license change is required), see section 4 of LGPLv3 text.

Cheers,
Mosè

Re: [julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Mosè Giordano]

Hi Kristoffer,

2016-04-27 14:32 GMT+02:00 Kristoffer Carlsson :
> This looks interesting.

Thank you!

> It is in my opinion a bit unfortunate that you chose
> to use the GPL license for this since it means that it will not be usable by
> most packages in Julia which are under MIT license.

Actually I used LGPL, not GPL, because the original Fortran library I
translated to Julia is released with that license, and, as far as I
know, you can you an LGPL program into a MIT one.

Bye,
Mosè

[julia-users] Re: [ANN] CmplxRoots.jl: Fast Complex Polynomial Root Finder

[Kristoffer Carlsson]

This looks interesting. It is in my opinion a bit unfortunate that you 
chose to use the GPL license for this since it means that it will not be 
usable by most packages in Julia which are under MIT license.

On Wednesday, April 27, 2016 at 1:15:35 PM UTC+2, Mosè Giordano wrote:
>
> Dear all,
>
> I am happy to announce the first version of CmplxRoots.jl 
> , a fast root finder of real 
> and complex polynomials.
>
> This is a Julia implementation of the algorithm described in 
>
>- J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver 
>and Its Further Optimization for Binary Microlenses", arXiv:1203.1034 
>
>
> I first wrote the package as a wrapper around the Fortran implementation 
> available at http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ but then 
> I rewrote it entirely in Julia, finding that it's even faster than 
> ccall'ing functions in the shared object.
>
> CmplxRoots.jl is a registered Julia package, so you can install it with 
> the package manager:
>
> Pkg.add("CmplxRoots")
>
> Two functions are exposed to the users:
>
> roots(polynomial[, roots, polish=true])
> roots5(polynomial[, roots])
>
> The first function can be used to solve polynomials of any order, the 
> second one is specialized and optimized for (and works only for) 
> polynomials of fifth-order, but I found it to be a bit slower than roots()
> .
>
> The only mandatory argument is the list of coefficients of the polynomial, 
> in ascending order.  So, if you want to find all the roots of polynomial
>
> (x - i)*(x - 2)*(x - 3*i)*(x - 4)*(x - 5*i) =
>   x^5 - (9*i + 6)*x^4 + (54*i - 15)*x^3 + (138 - 57*i)*x^2 - (184 + 90*i)*x + 
> 120*i
>
>
> you can use
>
> julia> roots([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im), 
> 1])
> 5-element Array{Complex{Float64},1}:
>  -1.55431e-15+5.0im
>  4.0+1.77636e-16im
>   1.55028e-15+3.0im
>  -1.04196e-16+1.0im
>  2.0-2.00595e-16im
>
> julia> roots5([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im
> ), 1])
> 5-element Array{Complex{Float64},1}:
>   5.92119e-16+5.0im
>   4.0-1.4803e-16im
>  2.0+1.88202e-16im
>  -1.88738e-15+3.0im
>   2.10942e-15+1.0im
>
> Optional roots argument is the array of starting values that can be used 
> to find the final roots (it won't be changed in-place).  Note that even if 
> the roots are all real, the output array is always of type Complex{Float64}.
>
> Also Polynomials.jl package has a root finding function, but according to 
> some random tests performed with @benchmark (from Benchmarks.jl package) 
> it's roughly 20 times slower than roots() function provided by 
> CmplxRoots.jl.  Of course do not expect results with machine precision (you 
> can have catastrophic cancellation even for simple quadratic polynomials), 
> but some random tests indicate that CmplxRoots.jl is slightly more precise 
> than Polynomials.jl (just compute abs2() of polynomials at the root).
>
> The package CmplxRoots.jl is licensed under the GNU Lesser General Public 
> License version 3 or any later version, as well as under a "customary 
> scientific license", which implies that if this code was important in the 
> scientific process or for the results of your scientific work, you are 
> asked for the appropriate citation of the paper Skowron & Gould 2012 (
> http://arxiv.org/abs/1203.1034).
>
> Happy root-finding,
> Mosè
>