zetaderiv: numerically unstable

[behackl]

The implementation of the derivative of the zeta function zetaderiv is numerically unstable and very slow (for large negative values):

sage: zetaderiv(1, CIF(-600))
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation

sage: timeit('zetaderiv(1, CIF(-500))')
5 loops, best of 3: 3.08 s per loop

Could either the current implementation be improved or an alternative numerical implementation be used?

---

This also causes errors when testing relations like

sage: bool(gamma(k*pi) * zetaderiv(1, k*pi*I)/log(2)^2 == 0)
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation

CC: @cheuberg @dkrenn @rwst

Component: numerical

Author: Benjamin Hackl

Branch/Commit: u/behackl/symbolic/test_relation/noconvergence @ 4ddf10c

Issue created by migration from https://trac.sagemath.org/ticket/20127

[behackl]

Branch: u/behackl/symbolic/test_relation/noconvergence

[behackl]

Commit: 4ddf10c

[behackl]

New commits:

4ddf10c

catch NoConvergence error in test_relation

[rwst]

comment:2

Looks fine, can you please add a doctest?

[behackl]

comment:3

In principle I could, yes. However, the simplest (reliable) doctest I can construct takes about 25 sec., and that's not really ideal. I'd really like to understand why this takes so much time.

[rwst]

comment:4

Which code would that be? Have you tried callgrind? I tried to confirm the random behaviour of the above but the first 16 tries all gave False. Do you have a reliable example?

[behackl]

comment:5

Replying to @rwst:

Which code would that be? Have you tried callgrind? I tried to confirm the random behaviour of the above but the first 16 tries all gave False. Do you have a reliable example?

Could you try it with the following expression? This is the original one which brought me to this error. If it works (which is rarely; maybe one out of 20 times on my laptop here), it takes several seconds.

sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)

[rwst]

comment:6

I get 10 times immediately False with develop. I can even do

sage: %timeit assert not bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1,1+k*pi*I/log(2))/log(2)^10 == 0)
10 loops, best of 3: 45.3 ms per loop

[rwst]

comment:7

I think the difference to your result is that I had a command k = var('k', domain='complex') instead of k = var('k'). With the latter I can confirm the issue. I thought I needed complex because I used k as integer before that. Now isn't the default complex anyway? It gets more mysterious.

[dkrenn]

comment:8

Replying to @behackl:

Could you try it with the following expression? This is the original one which brought me to this error. If it works (which is rarely; maybe one out of 20 times on my laptop here), it takes several seconds.

sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)

from datetime import datetime
var('k')

for _ in range(10):
    tic = datetime.now()
    try:
        result = bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)
    except Exception as e:
        print datetime.now()-tic, 'exception:', e
    else:
        print datetime.now()-tic, 'successful:', result

returns

0:00:11.445669 exception: zeta: too much cancellation
0:00:08.250215 exception: zeta: too much cancellation
0:00:03.706793 exception: zeta: too much cancellation
0:00:02.428823 successful: False
0:00:00.757367 successful: False
0:00:16.853168 exception: zeta: too much cancellation
0:00:00.228982 successful: False
0:00:03.824248 successful: False
0:00:37.852729 exception: zeta: too much cancellation
0:00:15.264686 exception: zeta: too much cancellation

Seems like I am lucky to get that much exceptions ;)

[rwst]

comment:9

See #14305. Variables may by default be complex in Pynac, but they are not in Maxima, which is (still) used by __nonzero__ (see #19040).