diff --git a/libsharp/planck.make b/libsharp/planck.make index fcf51df..8bbd7e1 100644 --- a/libsharp/planck.make +++ b/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)/%) diff --git a/libsharp/sharp.h b/libsharp/sharp.h index da6fc4a..4497dd4 100644 --- a/libsharp/sharp.h +++ b/libsharp/sharp.h @@ -40,5 +40,6 @@ #include "sharp_lowlevel.h" #include "sharp_legendre.h" +#include "sharp_legendre_roots.h" #endif diff --git a/libsharp/sharp_geomhelpers.c b/libsharp/sharp_geomhelpers.c index cd088b3..dbb44e0 100644 --- a/libsharp/sharp_geomhelpers.c +++ b/libsharp/sharp_geomhelpers.c @@ -32,6 +32,7 @@ #include #include "sharp_geomhelpers.h" +#include "sharp_legendre_roots.h" #include "c_utils.h" #include "ls_fft.h" #include @@ -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 +#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 + } diff --git a/libsharp/sharp_legendre_roots.h b/libsharp/sharp_legendre_roots.h new file mode 100644 index 0000000..2a056b2 --- /dev/null +++ b/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 diff --git a/python/libsharp/libsharp.pyx b/python/libsharp/libsharp.pyx index e622bb3..9e83d3d 100644 --- a/python/libsharp/libsharp.pyx +++ b/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 diff --git a/python/libsharp/tests/test_legendre.py b/python/libsharp/tests/test_legendre.py index 66b4c18..7860d2a 100644 --- a/python/libsharp/tests/test_legendre.py +++ b/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