Skip to content

Commit 7ecd1dd

Browse files
author
Dag Sverre Seljebotn
committed
Expose gauss_legendre_tbl publicly as gauss_legendre_roots
1 parent 765831e commit 7ecd1dd

7 files changed

Lines changed: 152 additions & 65 deletions

File tree

‎libsharp/planck.make‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ FULL_INCLUDE+= -I$(SD)
88
HDR_$(PKG):=$(SD)/*.h
99
LIB_$(PKG):=$(LIBDIR)/libsharp.a
1010
BIN:=sharp_testsuite
11-
LIBOBJ:=sharp_ylmgen_c.o sharp.o sharp_announce.o sharp_geomhelpers.o sharp_almhelpers.o sharp_core.o sharp_legendre.o
11+
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
1212
ALLOBJ:=$(LIBOBJ) sharp_testsuite.o
1313
LIBOBJ:=$(LIBOBJ:%=$(OD)/%)
1414
ALLOBJ:=$(ALLOBJ:%=$(OD)/%)

‎libsharp/sharp.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,5 +40,6 @@
4040

4141
#include "sharp_lowlevel.h"
4242
#include "sharp_legendre.h"
43+
#include "sharp_legendre_roots.h"
4344

4445
#endif

‎libsharp/sharp_geomhelpers.c‎

Lines changed: 2 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@
3232

3333
#include <math.h>
3434
#include "sharp_geomhelpers.h"
35+
#include "sharp_legendre_roots.h"
3536
#include "c_utils.h"
3637
#include "ls_fft.h"
3738
#include <stdio.h>
@@ -106,69 +107,6 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
106107
sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, NULL, weight, geom_info);
107108
}
108109

109-
static inline double one_minus_x2 (double x)
110-
{ return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
111-
112-
/* Function adapted from GNU GSL file glfixed.c
113-
Original author: Pavel Holoborodko (http://www.holoborodko.com)
114-
115-
Adjustments by M. Reinecke
116-
- adjusted interface (keep epsilon internal, return full number of points)
117-
- removed precomputed tables
118-
- tweaked Newton iteration to obtain higher accuracy */
119-
static void gauss_legendre_tbl(int n, double *x, double *w)
120-
{
121-
const double pi = 3.141592653589793238462643383279502884197;
122-
const double eps = 3e-14;
123-
int m = (n+1)>>1;
124-
125-
double t0 = 1 - (1-1./n) / (8.*n*n);
126-
double t1 = 1./(4.*n+2.);
127-
128-
#pragma omp parallel
129-
{
130-
int i;
131-
#pragma omp for schedule(dynamic,100)
132-
for (i=1; i<=m; ++i)
133-
{
134-
double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
135-
136-
int dobreak=0;
137-
int j=0;
138-
double dpdx;
139-
while(1)
140-
{
141-
double P_1 = 1.0;
142-
double P0 = x0;
143-
double dx, x1;
144-
145-
for (int k=2; k<=n; k++)
146-
{
147-
double P_2 = P_1;
148-
P_1 = P0;
149-
// P0 = ((2*k-1)*x0*P_1-(k-1)*P_2)/k;
150-
P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
151-
}
152-
153-
dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
154-
155-
/* Newton step */
156-
x1 = x0 - P0/dpdx;
157-
dx = x0-x1;
158-
x0 = x1;
159-
if (dobreak) break;
160-
161-
if (fabs(dx)<=eps) dobreak=1;
162-
UTIL_ASSERT(++j<100,"convergence problem");
163-
}
164-
165-
x[i-1] = -x0;
166-
x[n-i] = x0;
167-
w[i-1] = w[n-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
168-
}
169-
} // end of parallel region
170-
}
171-
172110
void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
173111
int stride_lon, int stride_lat, sharp_geom_info **geom_info)
174112
{
@@ -181,7 +119,7 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
181119
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
182120
int *stride_=RALLOC(int,nrings);
183121

184-
gauss_legendre_tbl(nrings,theta,weight);
122+
sharp_legendre_roots(nrings,theta,weight);
185123
for (int m=0; m<nrings; ++m)
186124
{
187125
theta[m] = acos(-theta[m]);

‎libsharp/sharp_legendre_roots.c‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
/* Function adapted from GNU GSL file glfixed.c
2+
Original author: Pavel Holoborodko (http://www.holoborodko.com)
3+
4+
Adjustments by M. Reinecke
5+
- adjusted interface (keep epsilon internal, return full number of points)
6+
- removed precomputed tables
7+
- tweaked Newton iteration to obtain higher accuracy */
8+
9+
#include <math.h>
10+
#include "sharp_legendre_roots.h"
11+
#include "c_utils.h"
12+
13+
static inline double one_minus_x2 (double x)
14+
{ return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
15+
16+
void sharp_legendre_roots(int n, double *x, double *w)
17+
{
18+
const double pi = 3.141592653589793238462643383279502884197;
19+
const double eps = 3e-14;
20+
int m = (n+1)>>1;
21+
22+
double t0 = 1 - (1-1./n) / (8.*n*n);
23+
double t1 = 1./(4.*n+2.);
24+
25+
#pragma omp parallel
26+
{
27+
int i;
28+
#pragma omp for schedule(dynamic,100)
29+
for (i=1; i<=m; ++i)
30+
{
31+
double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
32+
33+
int dobreak=0;
34+
int j=0;
35+
double dpdx;
36+
while(1)
37+
{
38+
double P_1 = 1.0;
39+
double P0 = x0;
40+
double dx, x1;
41+
42+
for (int k=2; k<=n; k++)
43+
{
44+
double P_2 = P_1;
45+
P_1 = P0;
46+
// P0 = ((2*k-1)*x0*P_1-(k-1)*P_2)/k;
47+
P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
48+
}
49+
50+
dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
51+
52+
/* Newton step */
53+
x1 = x0 - P0/dpdx;
54+
dx = x0-x1;
55+
x0 = x1;
56+
if (dobreak) break;
57+
58+
if (fabs(dx)<=eps) dobreak=1;
59+
UTIL_ASSERT(++j<100,"convergence problem");
60+
}
61+
62+
x[i-1] = -x0;
63+
x[n-i] = x0;
64+
w[i-1] = w[n-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
65+
}
66+
} // end of parallel region
67+
}

‎libsharp/sharp_legendre_roots.h‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
/*
2+
* This file is part of libsharp.
3+
*
4+
* libsharp is free software; you can redistribute it and/or modify
5+
* it under the terms of the GNU General Public License as published by
6+
* the Free Software Foundation; either version 2 of the License, or
7+
* (at your option) any later version.
8+
*
9+
* libsharp is distributed in the hope that it will be useful,
10+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12+
* GNU General Public License for more details.
13+
*
14+
* You should have received a copy of the GNU General Public License
15+
* along with libsharp; if not, write to the Free Software
16+
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17+
*/
18+
19+
/*
20+
* libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
21+
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22+
* (DLR).
23+
*/
24+
25+
/*! \file sharp_legendre_roots.h
26+
*
27+
* Copyright (C) 2006-2012 Max-Planck-Society
28+
* \author Martin Reinecke
29+
*/
30+
31+
#ifndef SHARP_LEGENDRE_ROOTS_H
32+
#define SHARP_LEGENDRE_ROOTS_H
33+
34+
#ifdef __cplusplus
35+
extern "C" {
36+
#endif
37+
38+
/*! Computes roots and Gaussian quadrature weights for Legendre polynomial
39+
of degree \a n.
40+
\param n Order of Legendre polynomial
41+
\param x Array of length \a n for output (root position)
42+
\param w Array of length \a w for output (weight for Gaussian quadrature)
43+
*/
44+
void sharp_legendre_roots(int n, double *x, double *w);
45+
46+
#ifdef __cplusplus
47+
}
48+
#endif
49+
50+
#endif

‎python/libsharp/libsharp.pyx‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
import numpy as np
22

3+
__all__ = ['legendre_transform', 'legendre_roots']
4+
35
cdef extern from "sharp.h":
46
ctypedef long ptrdiff_t
57

@@ -9,6 +11,7 @@ cdef extern from "sharp.h":
911
double *out, ptrdiff_t nx)
1012
void sharp_legendre_transform_recfac(double *r, ptrdiff_t lmax)
1113
void sharp_legendre_transform_recfac_s(float *r, ptrdiff_t lmax)
14+
void sharp_legendre_roots(int n, double *x, double *w)
1215

1316

1417
def legendre_transform(x, bl, out=None):
@@ -41,3 +44,13 @@ def _legendre_transform_s(float[::1] x, float[::1] bl, float[::1] out):
4144
sharp_legendre_transform_s(&bl[0], NULL, bl.shape[0] - 1, &x[0], &out[0], x.shape[0])
4245
return np.asarray(out)
4346

47+
48+
def legendre_roots(n):
49+
x = np.empty(n, np.double)
50+
w = np.empty(n, np.double)
51+
cdef double[::1] x_buf = x, w_buf = w
52+
if not (x_buf.shape[0] == w_buf.shape[0] == n):
53+
raise AssertionError()
54+
if n > 0:
55+
sharp_legendre_roots(n, &x_buf[0], &w_buf[0])
56+
return x, w

‎python/libsharp/tests/test_legendre.py‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
import numpy as np
22
from scipy.special import legendre
3+
from scipy.special import p_roots
34
import libsharp
45

56
from numpy.testing import assert_allclose
@@ -39,3 +40,20 @@ def test_legendre_transform():
3940
for ntheta in nthetas_to_try:
4041
for lmax in [0, 1, 2, 3, 20] + list(np.random.randint(50, size=4)):
4142
yield check_legendre_transform, lmax, ntheta
43+
44+
def check_legendre_roots(n):
45+
xs, ws = ([], []) if n == 0 else p_roots(n) # from SciPy
46+
xl, wl = libsharp.legendre_roots(n)
47+
assert_allclose(xs, xl)
48+
assert_allclose(ws, wl)
49+
50+
def test_legendre_roots():
51+
"""
52+
Test the Legendre root-finding algorithm from libsharp by comparing it with
53+
the SciPy version.
54+
"""
55+
yield check_legendre_roots, 0
56+
yield check_legendre_roots, 1
57+
yield check_legendre_roots, 32
58+
yield check_legendre_roots, 33
59+
yield check_legendre_roots, 128

0 commit comments

Comments
 (0)