#4132: [with patch, nearly positive review] complex arithmetic passes via pari
------------------------------+---------------------------------------------
Reporter: robertwb | Owner: AlexGhitza
Type: enhancement | Status: assigned
Priority: major | Milestone: sage-3.1.3
Component: basic arithmetic | Resolution:
Keywords: |
------------------------------+---------------------------------------------
Changes (by robertwb):
* summary: [with patch, needs review] complex arithmetic passes via pari
=> [with patch, nearly positive review] complex
arithmetic passes via pari
Comment:
Nice work. You've probably noticed by now that the trig operations are
expensive, so you could possibly shave 25% of some of the calls where you
compute both sinh and cosh by using the identity
cosh(x) = sqrt(1+sinh(x)).
On the note of expensive trig, might I suggest a faster sqrt
{{{
def new_sqrt(ComplexNumber self, all=False):
cdef ComplexNumber z = self._new()
if mpfr_zero_p(self.__im):
if mpfr_sgn(self.__re) >= 0:
mpfr_set_ui(z.__im, 0, rnd)
mpfr_sqrt(z.__re, self.__re, rnd)
else:
mpfr_set_ui(z.__re, 0, rnd)
mpfr_neg(z.__im, self.__re, rnd)
mpfr_sqrt(z.__im, z.__im, rnd)
if all:
return [z, -z] if z else [z]
else:
return z
# z = x + yi = (a+bi)^2
# expand, substitute, solve (note that a is nonzero)
# a^2 = (x + sqrt(x^2+y^2))/2
cdef mpfr_t a2
mpfr_init2(a2, self._prec)
mpfr_hypot(a2, self.__re, self.__im, rnd)
mpfr_add(a2, a2, self.__re, rnd)
mpfr_mul_2si(a2, a2, -1, rnd)
# a = sqrt(a2)
mpfr_sqrt(z.__re, a2, rnd)
# b = y/(2a)
mpfr_div(z.__im, self.__im, z.__re, rnd)
mpfr_mul_2si(z.__im, z.__im, -1, rnd)
mpfr_clear(a2)
if all:
return [z, -z]
else:
return z
}}}
which should run in 3 or 4 µs with 53 bits, and an even more significant
improvement further out.
This patch is very good, exactly what I wanted when I wrote the ticket, so
thanks. The only thing that ''needs'' to change for a positive review is I
noticed in several case you write {{{long(n)}}} when you want an unsigned
long. What this does is create a python int, call the python type "long"
on it to get a python long, then extract that as a c long. In summary,
just write {{{n}}} to get a C [unsigned] long just as you would in C.
(This is a common confusion with Cython, as the C long and python long are
totally different things.)
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/4132#comment:8>
Sage <http://sagemath.org/>
Sage - Open Source Mathematical Software: Building the Car Instead of
Reinventing the Wheel
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
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-trac?hl=en
-~----------~----~----~----~------~----~------~--~---
