Alternatively you could just use the implementation here:
http://www.ark.in-berlin.de/part.c
which seems to only rely on mpfr and GMP. However, since the Pari
version was based on it, I suspect it may also need correction.
He does seem to adjust the precision as he goes, but I've no idea how
far above or below the required precision this is. However there is no
need to use a precision of 55 from a certain point on. That probably
forces it to use two doubles instead of one in the floating point
computations, which I think would be unnecesary. You can look on
mathworld for how well the Rademacher series converges.
It would probably be quicker to do a fresh implementation from scratch
given the application in mind. The code there may not be GPL'd also.
To check it works, one should at least look at the congruences modulo
5, 7, 11 and 13. The author of the mpfr version does seem to check
them mod 11 but for very small n and only for a very few values of n.
If the gcds prove to be a bottleneck, I have some code that is roughly
twice as fast as GMP for computing these. Certainly the code in the
mpfr version is suboptimal and needs some serious optimisation when
the terms in the dedekind sum fit into a single limb, which on a 64
bit machine is probably always, when n is of a reasonable size.
Another suggestion is that for sufficiently small values of h and or
k, it may be quicker to compute the dedekind sum directly rather than
use the gcd formula.
A lookup table would also be useful for the tail end of the gcd/
dedekind sum computation.
I'd be shocked if the computation for n = 10^9 couldn't be done in
under 10s on sage.math.
Bill.
On 27 Jul, 00:39, Bill Hart <[EMAIL PROTECTED]> wrote:
> Here we go. Apparently g(h,q) = q*s(h,q) where s(h,q) is a dedekind
> sum.
>
> Then for h < q if r_0, r_1, ..., r_{n+1} are the remainders occurring
> in the Euclidean algorithm when computing gcd(h,q) (of which there
> should be about log q) then:
>
> s(h,q) = 1/12 * sum_{i=1}^{n+1}((-1)^(i+1) (r_i^2 + r_{i-1}^2 + 1) /
> (r_i * r_{i-1}) - ((-1)^n +1)/8)
>
> This, with the other ideas I gave above, will definitely let us beat
> Mathematica convincingly, as a simple back of the envelope calculation
> shows. It will have the added benefit of allowing us to beat Magma at
> something.
>
> Bill.
>
> On 27 Jul, 00:12, Bill Hart <[EMAIL PROTECTED]> wrote:
>
> > OK, apparently the Dedekind sums should be done using an algorithm due
> > to Apostol. I haven't tracked it down yet, but this is clearly what we
> > need.
>
> > Bill.
>
> > On 26 Jul, 23:37, Bill Hart <[EMAIL PROTECTED]> wrote:
>
> > > On 26 Jul, 23:18, Jonathan Bober <[EMAIL PROTECTED]> wrote:
>
> > > > s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
> > > > (1/k) sum_{j=1}^{k-1} j*floor(hj/k).
>
> > > > (To compute the floor in the inner sum, you can just use integer
> > > > division.)
>
> > > Yes, this will be faster. Of course integer division is not faster
> > > than floating point division, but since we are doing a sum from j = 1
> > > to k-1 we only need to know when floor(hj/k) increases by 1. This is
> > > simply a matter of adding h each iteration and seeing if we have gone
> > > over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
> > > and so on.
>
> > > This gets the inner computation down to about 8 cycles or so on
> > > average. Not enough to beat Mathematica though, unless I made a
> > > mistake in my back of the envelope computation somewhere.
>
> > > Bill.
--~--~---------~--~----~------------~-------~--~----~
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/
-~----------~----~----~----~------~----~------~--~---
