LevelS SHT library  3.50
sharp_legendre_roots.c
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  }
#define UTIL_ASSERT(cond, msg)
void sharp_legendre_roots(int n, double *x, double *w)

Generated on Mon Dec 10 2018 10:24:20 for LevelS SHT library