fixed bug in clpmn violating array bounds #482

Merged
merged 2 commits into from Mar 26, 2013

Projects

None yet

2 participants

@gertingold
Contributor

Hi,

it seems that there is a bug in the Fortran code of CLPMN where data are written Into the array CPM beyond its bounds provided that M=N. A corresponding bug originally present in the real version LPMN was fixed in 3d5294d3 about eight years ago. This bug concerns not only the associated Legendre polynomials but also the spherical harmonics for complex angles. In Python, the bug becomes apparent in the following minimal code:

'''
from scipy.special.specfun import clpmn
clpmn(1, 1, 5, 3)
'''

The third and fourth argument should not be very important as long as they do not correspond to an argument 1 or -1.

This bug was found by Michael Hartmann and analyzed together with him.

The change in line 275 in b73f9f9 should fix the problem.

Best regards,
Gert

@pv
Member
pv commented Mar 26, 2013

Thanks. Adding a test to test_basic.py should also be done.

@gertingold
Contributor

ok, I will try to come up with a test. It seems that so far there are no tests for the complex Legendre polynomials but only for the real ones. Should I add a test to the existing class TestLegendreFunctions (so far only real case or rather add a new class for the complex case?

@pv
Member
pv commented Mar 26, 2013

Adding it to the existing class should be OK. Here just checking e.g. the result from clpmn(1, 1, 5, 3) is probably enough to verify this fix --- writing a comprehensive test suite would take some more effort.

@gertingold
Contributor

I have added a test to test_basic.py but I must admit that I am not completely happy with it. First of all, the test itself runs ok even with the bug in CLPMN but segfaults at the end if the bug is not fixed. This clearly indicates the presence of the bug but it is not the way how the problem should appear in a test. Secondly, it would be nicer to calculate the expected array from the known formulae but this requires importing sqrt from cmath. I am not sure whether it would be acceptable to do such an import in test_basic.py.

In the test, I have chosen z=0.5+0.3j instead of 5+3j as an argument because then |z|<1 and everything is fine. For |z|>1, there seems to be an extra factor of 1j but I would have to look more deeply into the associated Legendre polynomials for complex arguments to be sure about it. This relates to your remark about writing more comprehensive tests for CLPMN.

@pv pv merged commit 5b75d9c into scipy:master Mar 26, 2013
@pv
Member
pv commented Mar 26, 2013

Thanks, looks good.

And you're absolutely correct, factors of i and -1 (in the real-valued code) apparently go wrong for |z| > 1, need to investigate -> http://projects.scipy.org/scipy/ticket/1877

@gertingold
Contributor

The following picture displaying the behavior of the result of CLPNM for N=M=1 on the imaginary axis gives further evidence that the behavior at |z|=1 is strange. After all, the associated Legendre polynomials should be analytic at z=i.

clpnm_on_imagaxis

The problem arises from LS in the code of CLPNM.

The above picture was generated with the following code:

import sys, os
from cmath import sqrt
import numpy as np
from scipy import special
from pyx import canvas, color, graph, style, text

"""This script shows that the evaluation of the associated Legendre polynomials
   in scipy by means of the function CLPNM has a problem when |z|>1. The graph
   shows the behavior of the real and imaginary part for n=m=1 on the imaginary
   axis. The result should be analytic at z=i, but a jump is observed. This
   behavior is due to the variable LS in the code of CLPNM. For comparison, the
   analytical continuation of the known result for real arguments is shown.
"""

xvals = np.linspace(0, 2, 500)
yvals_real = [special.specfun.clpmn(1, 1, 0, x)[0][1, 1].real for x in xvals]
yvals_imag = [special.specfun.clpmn(1, 1, 0, x)[0][1, 1].imag for x in xvals]
correct_real = [-sqrt(1+x*x).real for x in xvals]

c = canvas.canvas()
text.set(mode="latex")

g = graph.graphxy(width=8,
                  key=graph.key.key(pos="tl", dist=0.05),
                  x=graph.axis.linear(title="$y$"))
g.plot([graph.data.points(zip(xvals, yvals_real), x=1, y=2, title="CLPMN, real part"),
        graph.data.points(zip(xvals, yvals_imag), x=1, y=2, title="CLPMN, imaginary part"),
        graph.data.points(zip(xvals, correct_real), x=1, y=2, title="expected (is real)")],
        styles=[graph.style.line([style.linestyle.solid, color.gradient.Rainbow])])

c.insert(g)
c.text(g.xpos+0.5*g.width, g.ypos+g.height+1,
       u"CLPMN test on imaginary axis, purely imaginary argument $z=\mathrm{i}y$",
       [text.halign.center])
c.writeGSfile("%s.png" % os.path.splitext(sys.argv[0])[0],
              device="png256", resolution=100)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment