13 static inline double one_minus_x2 (
double x)
14 {
return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
18 const double pi = 3.141592653589793238462643383279502884197;
19 const double eps = 3e-14;
22 double t0 = 1 - (1-1./n) / (8.*n*n);
23 double t1 = 1./(4.*n+2.);
28 #pragma omp for schedule(dynamic,100) 31 double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
42 for (
int k=2; k<=n; k++)
47 P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
50 dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
58 if (fabs(dx)<=eps) dobreak=1;
64 w[i-1] = w[n-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
#define UTIL_ASSERT(cond, msg)
void sharp_legendre_roots(int n, double *x, double *w)