I added the following function to my polynomial_element.py file. The good news
is that tests on a small polynomial (order 7) ran 1k to 2k times faster than
__call__(). This thoroughly blows matlab out of the water for an inlined
function, and is about 10x faster than a hard-coded polynomial.
Sadly, whether or not those results scale became hard to test because pyrex
seems to barf on larger polynomials using this. (my guess is that it doesn't
like the 100+ nested parentheses)
So I guess my question is, should we expect pyrex to take an arbitrary number
of nested parentheses, or is this a job for Josh's inline c idea? (because
clearly, abandoning 3 orders of magnitude speedup isn't an option)
Test code:
{{{
%pyrex
def f(double x): #generated by p.pyrex()
return(((((((1.0)*x+0.0)*x+0.0)*x+0.0)*x+3.0)*x+1.0)*x+1.0)*x+2.0
def test_f(int n,x):
for i in range(n):
v = f(x)
def test_p(int n, p,x):
for i in range(n):
v = p(x)
}}}
{{{
a = x^7 + 3*x^3 + x^2 + x + 2
}}}
{{{
%time test_f(100000,RDF(100.123134651346))
///
CPU time: 0.04 s, Wall time: 0.05 s
}}}
{{{
%time test_p(1000,a,RDF(100.123134651346))
///
CPU time: 0.70 s, Wall time: 0.73 s
}}}
Code added to Polynomial class:
def pyrex(self,name='f'):
"""
Return a string of valid pyrex code which is a function
that quickly (~1k times faster than pure python) evaluates
the polynomial as a double, at a double value of x. This
borrows code from the above __call__ to "unroll" the loop
that uses Horner's rule into a long string:
def f(double x):
return ((...((p[d])*x+p[d-1])*x...)*x+p[0]
For example,
x^7 + 3x^4 + 2x+1
gets unrolled to (output cleaned up for aesthetic reasons):
def f(double x):
(((((((1.0)*x+0.0)*x+0.0)*x+3.0)*x+0.0)*x+0.0)*x+2.0)*x+1.0
INPUT:
name -- the name in the pyrex function that evaluates
this polynomial as a double
OUTPUT:
a valid pyrex string
AUTHORS:
-- David Joyner, William Stein, wrote __call__
-- Tom Boothby: refactored DJ & WS code to return pyrex
"""
d = self.degree()
result = "%.16e"%self[d]
i = d - 1
while i >= 0:
result = "(%s)*%s + %.16e"%(result,'x',self[i])
i -= 1
return """def %s(double x): return %s"""%(name,result)
--~--~---------~--~----~------------~-------~--~----~
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/
-~----------~----~----~----~------~----~------~--~---
