#4225: faster sqrt for complex numbers
----------------------+-----------------------------------------------------
Reporter: robertwb | Owner: tbd
Type: defect | Status: new
Priority: major | Milestone: sage-3.1.3
Component: algebra | Resolution:
Keywords: |
----------------------+-----------------------------------------------------
Comment (by robertwb):
Code:
{{{
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
# self = x + yi = (a+bi)^2
# expand, substitute, solve
# a^2 = (x + sqrt(x^2+y^2))/2
cdef bint avoid_branch = mpfr_sgn(self.__re) < 0 and
mpfr_cmpabs(self.__im, self.__re) < 0
cdef mpfr_t a2
mpfr_init2(a2, self._prec)
mpfr_hypot(a2, self.__re, self.__im, rnd)
if avoid_branch:
# x + sqrt(x^2+y^2) numerically unstable for x near negative real
axis
# so we compute sqrt of (-z) and shift by i at the end
mpfr_sub(a2, a2, self.__re, rnd)
else:
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 avoid_branch:
# Now we shift by i depending on what side of the branch we were
on.
mpfr_swap(z.__re, z.__im)
# Note that y was never negated, so b already has the opposite
sign.
if mpfr_sgn(self.__im) < 0:
mpfr_neg(z.__re, z.__re, rnd)
mpfr_neg(z.__im, z.__im, rnd)
if all:
return [z, -z]
else:
return z
}}}
I will post a patch (with more documentation) when I have a copy of sage
with #4132 included.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/4225#comment:1>
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
-~----------~----~----~----~------~----~------~--~---
