Expose gauss_legendre_tbl publicly as gauss_legendre_roots

libsharp/planck.make

@@ -8,7 +8,7 @@ FULL_INCLUDE+= -I$(SD)
 HDR_$(PKG):=$(SD)/*.h
 LIB_$(PKG):=$(LIBDIR)/libsharp.a
 BIN:=sharp_testsuite
-LIBOBJ:=sharp_ylmgen_c.o sharp.o sharp_announce.o sharp_geomhelpers.o sharp_almhelpers.o sharp_core.o sharp_legendre.o
+LIBOBJ:=sharp_ylmgen_c.o sharp.o sharp_announce.o sharp_geomhelpers.o sharp_almhelpers.o sharp_core.o sharp_legendre.o sharp_legendre_roots.o
 ALLOBJ:=$(LIBOBJ) sharp_testsuite.o
 LIBOBJ:=$(LIBOBJ:%=$(OD)/%)
 ALLOBJ:=$(ALLOBJ:%=$(OD)/%)

libsharp/sharp.h

@@ -40,5 +40,6 @@
 
 #include "sharp_lowlevel.h"
 #include "sharp_legendre.h"
+#include "sharp_legendre_roots.h"
 
 #endif

libsharp/sharp_geomhelpers.c

@@ -32,6 +32,7 @@
 
 #include <math.h>
 #include "sharp_geomhelpers.h"
+#include "sharp_legendre_roots.h"
 #include "c_utils.h"
 #include "ls_fft.h"
 #include <stdio.h>
@@ -106,69 +107,6 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
   sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, NULL, weight, geom_info);
   }
 
-static inline double one_minus_x2 (double x)
-  { return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
-
-/* Function adapted from GNU GSL file glfixed.c
-   Original author: Pavel Holoborodko (http://www.holoborodko.com)
-
-   Adjustments by M. Reinecke
-    - adjusted interface (keep epsilon internal, return full number of points)
-    - removed precomputed tables
-    - tweaked Newton iteration to obtain higher accuracy */
-static void gauss_legendre_tbl(int n, double *x, double *w)
-  {
-  const double pi = 3.141592653589793238462643383279502884197;
-  const double eps = 3e-14;
-  int m = (n+1)>>1;
-
-  double t0 = 1 - (1-1./n) / (8.*n*n);
-  double t1 = 1./(4.*n+2.);
-
-#pragma omp parallel
-{
-  int i;
-#pragma omp for schedule(dynamic,100) 
-  for (i=1; i<=m; ++i)
-    {
-    double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
-
-    int dobreak=0;
-    int j=0;
-    double dpdx;
-    while(1)
-      {
-      double P_1 = 1.0;
-      double P0 = x0;
-      double dx, x1;
-
-      for (int k=2; k<=n; k++)
-        {
-        double P_2 = P_1;
-        P_1 = P0;
-//        P0 = ((2*k-1)*x0*P_1-(k-1)*P_2)/k;
-        P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
-        }
-
-      dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
-
-      /* Newton step */
-      x1 = x0 - P0/dpdx;
-      dx = x0-x1;
-      x0 = x1;
-      if (dobreak) break;
-
-      if (fabs(dx)<=eps) dobreak=1;
-      UTIL_ASSERT(++j<100,"convergence problem");
-      }
-
-    x[i-1] = -x0;
-    x[n-i] = x0;
-    w[i-1] = w[n-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
-    }
-} // end of parallel region
-  }
-
 void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
   int stride_lon, int stride_lat, sharp_geom_info **geom_info)
   {
@@ -181,7 +119,7 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
   ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
   int *stride_=RALLOC(int,nrings);
 
-  gauss_legendre_tbl(nrings,theta,weight);
+  sharp_legendre_roots(nrings,theta,weight);
   for (int m=0; m<nrings; ++m)
     {
     theta[m] = acos(-theta[m]);

libsharp/sharp_legendre_roots.c

@@ -0,0 +1,67 @@
+/* Function adapted from GNU GSL file glfixed.c
+   Original author: Pavel Holoborodko (http://www.holoborodko.com)
+
+   Adjustments by M. Reinecke
+    - adjusted interface (keep epsilon internal, return full number of points)
+    - removed precomputed tables
+    - tweaked Newton iteration to obtain higher accuracy */
+
+#include <math.h>
+#include "sharp_legendre_roots.h"
+#include "c_utils.h"
+
+static inline double one_minus_x2 (double x)
+  { return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
+
+void sharp_legendre_roots(int n, double *x, double *w)
+  {
+  const double pi = 3.141592653589793238462643383279502884197;
+  const double eps = 3e-14;
+  int m = (n+1)>>1;
+
+  double t0 = 1 - (1-1./n) / (8.*n*n);
+  double t1 = 1./(4.*n+2.);
+
+#pragma omp parallel
+{
+  int i;
+#pragma omp for schedule(dynamic,100)
+  for (i=1; i<=m; ++i)
+    {
+    double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
+
+    int dobreak=0;
+    int j=0;
+    double dpdx;
+    while(1)
+      {
+      double P_1 = 1.0;
+      double P0 = x0;
+      double dx, x1;
+
+      for (int k=2; k<=n; k++)
+        {
+        double P_2 = P_1;
+        P_1 = P0;
+//        P0 = ((2*k-1)*x0*P_1-(k-1)*P_2)/k;
+        P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
+        }
+
+      dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
+
+      /* Newton step */
+      x1 = x0 - P0/dpdx;
+      dx = x0-x1;
+      x0 = x1;
+      if (dobreak) break;
+
+      if (fabs(dx)<=eps) dobreak=1;
+      UTIL_ASSERT(++j<100,"convergence problem");
+      }
+
+    x[i-1] = -x0;
+    x[n-i] = x0;
+    w[i-1] = w[n-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
+    }
+} // end of parallel region
+  }

libsharp/sharp_legendre_roots.h

@@ -0,0 +1,50 @@
+/*
+ *  This file is part of libsharp.
+ *
+ *  libsharp is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libsharp is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libsharp; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file sharp_legendre_roots.h
+ *
+ *  Copyright (C) 2006-2012 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#ifndef SHARP_LEGENDRE_ROOTS_H
+#define SHARP_LEGENDRE_ROOTS_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*! Computes roots and Gaussian quadrature weights for Legendre polynomial
+    of degree \a n.
+    \param n Order of Legendre polynomial
+    \param x Array of length \a n for output (root position)
+    \param w Array of length \a w for output (weight for Gaussian quadrature)
+ */
+void sharp_legendre_roots(int n, double *x, double *w);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

python/libsharp/libsharp.pyx

@@ -1,5 +1,7 @@
 import numpy as np
 
+__all__ = ['legendre_transform', 'legendre_roots']
+
 cdef extern from "sharp.h":
     ctypedef long ptrdiff_t
 
@@ -9,6 +11,7 @@ cdef extern from "sharp.h":
                                   double *out, ptrdiff_t nx)
     void sharp_legendre_transform_recfac(double *r, ptrdiff_t lmax)
     void sharp_legendre_transform_recfac_s(float *r, ptrdiff_t lmax)
+    void sharp_legendre_roots(int n, double *x, double *w)
 
 
 def legendre_transform(x, bl, out=None):
@@ -41,3 +44,13 @@ def _legendre_transform_s(float[::1] x, float[::1] bl, float[::1] out):
     sharp_legendre_transform_s(&bl[0], NULL, bl.shape[0] - 1, &x[0], &out[0], x.shape[0])
     return np.asarray(out)
 
+
+def legendre_roots(n):
+    x = np.empty(n, np.double)
+    w = np.empty(n, np.double)
+    cdef double[::1] x_buf = x, w_buf = w
+    if not (x_buf.shape[0] == w_buf.shape[0] == n):
+        raise AssertionError()
+    if n > 0:
+        sharp_legendre_roots(n, &x_buf[0], &w_buf[0])
+    return x, w

python/libsharp/tests/test_legendre.py

@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.special import legendre
+from scipy.special import p_roots
 import libsharp
 
 from numpy.testing import assert_allclose
@@ -39,3 +40,20 @@ def test_legendre_transform():
     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
+
+def check_legendre_roots(n):
+    xs, ws = ([], []) if n == 0 else p_roots(n) # from SciPy
+    xl, wl = libsharp.legendre_roots(n)
+    assert_allclose(xs, xl)
+    assert_allclose(ws, wl)
+
+def test_legendre_roots():
+    """
+    Test the Legendre root-finding algorithm from libsharp by comparing it with
+    the SciPy version.
+    """
+    yield check_legendre_roots, 0
+    yield check_legendre_roots, 1
+    yield check_legendre_roots, 32
+    yield check_legendre_roots, 33
+    yield check_legendre_roots, 128