Discussion that should be archived in sage-support. (It's not a question at
this point.)
------- Forwarded message -------
From: "Kiran S. Kedlaya" <[EMAIL PROTECTED]>
To: "William Stein" <[EMAIL PROTECTED]>
Cc: "Joshua Kantor" <[EMAIL PROTECTED]>
Subject: Re: Fwd: GSL
Date: Sat, 27 Jan 2007 05:55:23 -0800
Typo: those should be result[2*i] and result[2*i+1] instead of result[i]
and result[i+1] in line -3. But anyway, this seems to run much faster
than numpy, which is what I was expecting. (Did I mention that the
ability to effortlessly swap in different software components is what I
really love about SAGE?)
There is a section in the SAGE Reference Manual (section 9) about
library interfaces; it might be good to include a subsection about GSL,
maybe with this example included (as corrected above).
Thanks,
Kiran
William Stein wrote:
> Here's a plain text ("edit") version of the worksheet...
>
> {{{
> %sagex
> include 'gsl.pxi'
> import sage.rings.complex_double
> def poly_root_finder(f):
> cdef double* coeff_list
> cdef double* result
> cdef int n,i
> cdef gsl_poly_complex_workspace* w
> l = f.coeffs()
> n = len(l)
> coeff_list = <double *>malloc(sizeof(double)*n)
> result = <double *> malloc(2*(n-1)*sizeof(double))
> for i from 0<=i<n:
> coeff_list[i] = float(l[i])
> w = <gsl_poly_complex_workspace*> gsl_poly_complex_workspace_alloc(n)
> gsl_poly_complex_solve(coeff_list,n,w,result)
> gsl_poly_complex_workspace_free(w)
> free(coeff_list)
> root_list= [ sage.rings.complex_double.CDF(result[i],result[i+1]) for i
> in range(n-1)]
> free(result)
> return root_list
> }}}
>
> {{{
> f=x^5 -1
> }}}
>
> {{{
> roots = poly_root_finder(f)
> }}}
>
> {{{
> roots
> ///
> [-0.809016994375 + 0.587785243988*I, 0.587785252292 - 0.809017002583*I,
> -0.809016994375 - 0.587785243988*I, -0.587785252292 + 0.309017002583*I,
> 0.309016994375 + 0.951056540012*I]
> }}}
>
> {{{
> f(roots[0])
> ///
> 2.22044604925e-16 + 8.98261499904e-17*I
> }}}
>
> {{{
>
> }}}
>
> On Sat, 27 Jan 2007 00:10:24 -0800, Joshua Kantor <[EMAIL PROTECTED]> wrote:
>
>
>> Hello Kiran,
>>
>> Sage doesn't really have an interface to GSL, as GSL is
>> just a c library that can be used from within pyrex.
>> I wrote a simple function poly_root_finder that
>> uses GSL's polynomial root finder to find all
>> the roots of a polynomial, they are returned
>> as a list of complex numbers. There is an example
>> in the worksheet of its usage. Let me know what the results
>> of your comparison are, or if you have any questions.
>>
>> Josh
>>
>>
>> On Fri, 26 Jan 2007, William Stein wrote:
>>
>>
>>> You might be able to say something helpful...
>>>
>>> ------- Forwarded message -------
>>> From: "Kiran S. Kedlaya" <[EMAIL PROTECTED]>
>>> To: "William Stein" <[EMAIL PROTECTED]>
>>> Cc:
>>> Subject: GSL
>>> Date: Fri, 26 Jan 2007 16:00:09 -0800
>>>
>>> Is there any documentation available about SAGE's interface to GSL?
>>> Specifically, I'm trying to figure out how to access the polynomial root
>>> finder in GSL, so I can compare it to numpy.
>>>
>>> Kiran
>>>
>>>
>>>
>
>
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://sage.math.washington.edu/sage/ and http://sage.scipy.org/sage/
-~----------~----~----~----~------~----~------~--~---
