That looks brilliant. I think it would be fine to have the faster version in unreadable mpfr provided that the "real" formula was included as comments.
Robert, can you post some version of what you have so far so I can see how you have made it faster? About principal vs. optimal. Neither are standard terms, I invented them yesterday. Don't be tricked into thinking that the "principal" value is in any way canonical. Only the "optimal" version has that claim. John 2009/12/18 Robert Bradshaw <rober...@math.washington.edu>: > > On Dec 17, 2009, at 9:58 PM, William Stein wrote: > >> On Thu, Dec 17, 2009 at 9:35 PM, David Roe <r...@math.harvard.edu> >> wrote: >>> Hey John, >>> I worked on it tonight, and I'm not sure how much you want to >>> optimize it. >>> Is a factor of 2 or 3 speedup worth making the code much less >>> readable (I'd >>> include comments, but...)? I could also probably improve a few >>> things and >>> get maybe 10% and leave it mostly as is. >>> David >> >> I posted some remarks on the ticket. My initial benchmarks suggest >> that John's implementation of AGM may be about 10 times slower than >> PARI's, and that Magma's (which is only in the real case) may be 10 >> times faster than PARI... making John's 100 times slower than Pari? >> Hopefully I'm wrong, but that's what I get. So I hope there is a way >> to speed it up by more than a factor of 2 or 3. > > I implemented the optimal and principle options for CDF, which are > very fast (15x pari): > > sage: a = CDF(-0.95,-0.65) > sage: b = CDF(0.683,0.747) > sage: a.agm(b) > 0.338175462986 - 0.0135326969565*I > sage: a.agm(b, algorithm='optimal') > -0.371591652352 + 0.319894660207*I > sage: a.agm(b, algorithm='pari') > 0.080689185076 + 0.239036532686*I > sage: %timeit a.agm(b) > 100000 loops, best of 3: 2.24 µs per loop > sage: %timeit a.agm(b, algorithm='optimal') > 100000 loops, best of 3: 2.42 µs per loop > sage: %timeit a.agm(b, algorithm='pari') > 10000 loops, best of 3: 32.2 µs per loop > sage: %timeit a.real() # nearly trivial call for > comparison > 1000000 loops, best of 3: 575 ns per loop > > I was also able to shave off about 30 (out of 180) microseconds for > your agm via a more efficient versions of > > one = self.parent().one_element() > two = one+one > half = one/two > eps = two**(2-self.prec()) > > One could write the algorithm directly against mpfr (I'm guessing for > much more than a 2-3x speedup), but that would make for really hard to > read code. Another option might be to use the CDF for <=53-bit (it > only performs basic arithmatic, so IEEE standards should be just as > good as mpfr), and the more generic code for larger precisions (where, > at the limit, the overhead won't matter at all). Wouldn't help much > for intermediate precisions though. > > Another note on the code--inplace operations aren't actually inplace > (due to immutability of real numbers). > > - Robert > > -- > To post to this group, send an email to sage-devel@googlegroups.com > To unsubscribe from this group, send an email to > sage-devel+unsubscr...@googlegroups.com > For more options, visit this group at > http://groups.google.com/group/sage-devel > URL: http://www.sagemath.org > -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
