I wrote up some quick code to use GSL's polynomial solver--see below.
It looks like GSL is about 4 times as fast as numpy, but seem to give
about the same precision (on my test polynomial).
One may want to use newton-raphson to refine these roots, but I'm
pretty sure one would have to work with more than 53 bits of
precision internally to do so.
- Robert
On Sep 4, 2007, at 5:30 PM, Robert Miller wrote:
>>> The root finding is much much more accurate, and no doubt faster
>>> too.
>>
>> I really hope this is true. Have you tested it?
>
> I haven't. However, I remember testing the same example that was in
> GSL's documentation. It was something like x^3-1, something that had 1
> as a root. The GSL root ~ 1 was accurate to 10 or 12 places, and the
> numpy root only to about 5. I have no reason to assume GSL is faster
> except that numpy converts the problem to an eigenvalue problem first,
> and GSL doesn't, and numpy is python, GSL is C.
-------------------------------------------------------
{{{id=14|
%cython
from sage.rings.all import CDF
cdef extern from "gsl/gsl_poly.h":
ctypedef struct gsl_poly_complex_workspace:
size_t nc
double * matrix
gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(int n)
void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w)
int gsl_poly_complex_solve(double *a, size_t n,
gsl_poly_complex_workspace *w, double* z)
def roots(f):
f = list(f)
cdef Py_ssize_t i, n = len(f)-1
cdef double* roots = <double*>sage_malloc(sizeof(double)*n*2)
cdef double* coeffs = <double*>sage_malloc(sizeof(double)*(n+1))
if roots == NULL or coeffs == NULL:
raise MemoryError
cdef gsl_poly_complex_workspace *w =
gsl_poly_complex_workspace_alloc(n+1)
for i from 0 <= i <= n:
coeffs[i] = f[i]
gsl_poly_complex_solve(coeffs, n+1, w, roots)
gsl_poly_complex_workspace_free(w)
res = [(roots[2*i], roots[2*i+1]) for i from 0 <= i < n]
sage_free(roots)
sage_free(coeffs)
return res
}}}
{{{id=15|
x = ZZ['x'].0
f = prod([x-i for i in range(1,12)]); f
///
x^11 - 66*x^10 + 1925*x^9 - 32670*x^8 + 357423*x^7 - 2637558*x^6 +
13339535*x^5 - 45995730*x^4 + 105258076*x^3 - 150917976*x^2 +
120543840*x - 39916800
}}}
{{{id=16|
%time
for _ in range(10^4): z = roots(f)
///
CPU time: 0.82 s, Wall time: 0.87 s
}}}
{{{id=17|
z
///
[(0.9999999999995961, 0.0), (2.0000000000088503, 0.0),
(2.9999999998633893, 0.0), (4.0000000011247705, 0.0),
(4.9999999950367791, 0.0), (6.0000000126757103, 0.0),
(6.9999999802111796, 0.0), (8.0000000193450411, 0.0),
(8.9999999882746984, 0.0), (10.000000004115627, 0.0),
(10.999999999344341, 0.0)]
}}}
{{{id=19|
# plot(f, xmin=0, xmax=10).show()
}}}
{{{id=21|
g = RDF['x'](list(reversed(list(f)))) # reverse coeffs due to bug
(now fixed)
}}}
{{{id=23|
g.roots()
///
[10.9999999991, 10.0000000049, 8.9999999883, 8.00000001573,
6.99999998694, 6.00000000701, 4.99999999753, 4.00000000055,
2.99999999993, 2.0, 1.0]
}}}
{{{id=24|
%time
for _ in range(10^4):
z = g.roots()
///
CPU time: 3.65 s, Wall time: 3.73 s
}}}
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
