legendre transforms: Polishing and bugfixes

vbroadcast didn't work properly for some reason, vload does..

python/libsharp/libsharp.pyx

@@ -14,7 +14,9 @@ cdef extern from "sharp.h":
 def legendre_transform(x, bl, out=None):
     if out is None:
         out = np.empty_like(x)
-    if x.dtype == np.float64:
+    if out.shape[0] == 0:
+        return out
+    elif x.dtype == np.float64:
         if bl.dtype != np.float64:
             bl = bl.astype(np.float64)
         return _legendre_transform(x, bl, out=out)

python/libsharp/tests/test_legendre.py

@@ -2,14 +2,17 @@ import numpy as np
 from scipy.special import legendre
 import libsharp
 
+from numpy.testing import assert_allclose
 
-def test_legendre_transform():
-    lmax = 20
-    ntheta = 19
 
+def check_legendre_transform(lmax, ntheta):
     l = np.arange(lmax + 1)
-    bl = np.exp(-l*(l+1))
-    bl *= (2 * l + 1)
+    if lmax >= 1:
+        sigma = -np.log(1e-3) / lmax / (lmax + 1)
+        bl = np.exp(-sigma*l*(l+1))
+        bl *= (2 * l + 1)
+    else:
+        bl = np.asarray([1], dtype=np.double)
 
     theta = np.linspace(0, np.pi, ntheta, endpoint=True)
     x = np.cos(theta)
@@ -20,10 +23,19 @@ def test_legendre_transform():
         P[:, l] = legendre(l)(x)
     y0 = np.dot(P, bl)
 
+
     # double-precision
     y = libsharp.legendre_transform(x, bl)
-    assert np.max(np.abs(y - y) / np.abs(y)) < 1e-12
+
+    assert_allclose(y, y0, rtol=1e-12, atol=1e-12)
 
     # single-precision
     y32 = libsharp.legendre_transform(x.astype(np.float32), bl)
-    assert np.max(np.abs(y32 - y) / np.abs(y32)) < 1e-4
+    assert_allclose(y, y0, rtol=1e-5, atol=1e-5)
+
+
+def test_legendre_transform():
+    nthetas_to_try = [0, 9, 17, 19] + list(np.random.randint(500, size=20))
+    for ntheta in nthetas_to_try:
+        for lmax in [0, 1, 2, 3, 20] + list(np.random.randint(50, size=4)):
+            yield check_legendre_transform, lmax, ntheta