#7719: Improvements to complex AGM
--------------------------------+-------------------------------------------
Reporter: cremona | Owner: AlexGhitza
Type: enhancement | Status: new
Priority: major | Milestone: sage-4.3.1
Component: basic arithmetic | Keywords: complex agm
Work_issues: | Author:
Upstream: N/A | Reviewer:
Merged: |
--------------------------------+-------------------------------------------
This is related to #6021 but also of independent interest.
As of 4.3 we use pari to compute the complex agm (i.e. a.agm(b) where a, b
are complex). Now the complex agm is multi-valued, and for the
application to periods of elliptic curves it does matter which value is
used (see upcoming paper by John Cremona and Thotsaphon Thongjunthug).
The pari-computed value is not this "optimal" value. So I have
implemented a native Sage version, replacing the existing code in
sage/rings/complex_number.pyx.
Much to my surprise, the new code is in some cases 50 times faster than
calling the pari library -- despite the fact that I have not yet cython-
optimised the code!
{{{
sage: CC = ComplexField(200)
sage: a = CC(-0.95,-0.65)
sage: b = CC(0.683,0.747)
sage: %timeit t = a.agm(b, algorithm="pari")
100 loops, best of 3: 7.04 ms per loop
sage: %timeit t = a.agm(b, algorithm="principal")
10000 loops, best of 3: 136 mus per loop
sage: %timeit t = a.agm(b, algorithm="optimal")
10000 loops, best of 3: 146 mus per loop
}}}
Here "mus" means microseconds (this was run on a computer for which
displaying greek "mu" causes an error). "pari" is the old way calling the
pari library function; "principal" is a native implementation of
essentially the same; "optimal" is the native implementation returning the
so-called optimal value.
Some details: AGM(a,b) is the common limit of two sequences {{{a_n,b_n}}}
under the iteration {{{(a,b) -> ((a+b)/2, sqrt(a*b))}}} and the issue is
which square root to take. The complete story is a wonderful but quite
long one (which started with Gauss). Essentially, the "principal"
algorithm always takes the principal branch of the square root (with
positive real part); pari does the same after an initial step where
AGM(a,b) is replaced by a*AGM(1/b/a); the optimal sequence (which gives
the largest limit) is the one for which the sign is always chosen so that
sqrt(a*b) is closest to (a+b)/2. Note that the optimal sequence is
preserved under scaling, so gives AGM(z*a,z*b)=z*AGM(a,b), but this is not
true of the principal sequence.
I have a patch, but before posting it I'll ask for some help optimising it
cythonically.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/7719>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
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.
