18 #define EPS 100*DBL_EPSILON
22 #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a))
32 #define POLYNOMIAL_DEFLATION 0
36 #ifdef POLYNOMIAL_DEFLATION
40 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
46 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
58 static void Jacobz (
const int n,
double *z,
const double alpha,
68 static void TriQL(
const int,
double *,
double *,
double **);
70 double gammaF (
const double);
72 static void RecCoeff(
const int,
double *,
double *,
const double,
76 void JKMatrix(
int,
double *,
double *);
78 void chri1(
int,
double*,
double*,
double*,
double*,
double);
102 void zwgj (
double *z,
double *w,
const int np,
const double alpha,
110 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
116 jacobd (np,z,w,np,alpha,beta);
120 fac = pow(two,apb + one)*
gammaF(alpha + np + one)*
gammaF(beta + np + one);
126 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
158 void zwgrjm(
double *z,
double *w,
const int np,
const double alpha,
178 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
184 jacobz (np-1,z+1,alpha,beta+1);
186 jacobfd (np,z,w,NULL,np-1,alpha,beta);
190 fac = pow(two,apb)*
gammaF(alpha + np)*
gammaF(beta + np);
196 for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
198 w[0] *= (beta + one);
234 void zwgrjp(
double *z,
double *w,
const int np,
const double alpha,
254 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
258 jacobz (np-1,z,alpha+1,beta);
262 jacobfd (np,z,w,NULL,np-1,alpha,beta);
266 fac = pow(two,apb)*
gammaF(alpha + np)*
gammaF(beta + np);
272 for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
274 w[np-1] *= (alpha + one);
306 void zwglj(
double *z,
double *w,
const int np,
const double alpha,
326 double fac, one = 1.0, apb = alpha + beta, two = 2.0;
334 jacobz (np-2,z + 1,alpha + one,beta + one);
336 jacobfd (np,z,w,NULL,np-1,alpha,beta);
340 fac = pow(two,apb + 1)*
gammaF(alpha + np)*
gammaF(beta + np);
342 fac /= (np-1)*
gammaF(np)*
gammaF(alpha + beta + np + one);
346 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
348 w[0] *= (beta + one);
350 w[np-1] *= (alpha + one);
380 void zwgk(
double *z,
double *w,
const int npt ,
const double alpha,
398 int kpoints = 2*np + 1;
404 int ncoeffs = (int)floor(3.0*(np+1)/2);
414 double *a =
new double[kpoints];
416 double *b =
new double[kpoints];
422 for(i = 0; i < kpoints; i++)
448 double** zmatrix =
new double*[kpoints];
450 for(i = 0; i < kpoints; i++)
454 zmatrix[i] =
new double[kpoints];
456 for(j = 0; j < kpoints; j++)
466 for(i = 0; i < kpoints; i++)
478 TriQL(kpoints, a, b, zmatrix);
482 for(i = 0; i < kpoints; i++)
514 void zwrk(
double* z,
double* w,
const int npt ,
const double alpha,
530 fprintf(stderr,
"too few points in formula\n");
554 int ncoeffs = (int)ceil(3.0*np/2);
560 double *a =
new double[ncoeffs+1];
562 double *b =
new double[ncoeffs+1];
568 for(i = 0; i < ncoeffs+1; i++)
586 double* a0 =
new double[ncoeffs];
588 double* b0 =
new double[ncoeffs];
592 chri1(ncoeffs,a,b,a0,b0,end0);
596 double s = b0[0]/fabs(b0[0]);
604 double* z1 =
new double[2*np-1];
606 double* w1 =
new double[2*np-1];
608 for(i = 0; i < ncoeffs; i++)
622 double** zmatrix =
new double*[2*np-1];
624 for(i = 0; i < 2*np-1; i++)
628 zmatrix[i] =
new double[2*np-1];
630 for(j = 0; j < 2*np-1; j++)
640 for(i = 0; i < 2*np-1; i++)
652 TriQL(2*np-1, z1, w1, zmatrix);
658 for(i = 0; i < 2*np-1; i++)
662 w1[i] = s*w1[i]/(z1[i]-end0);
674 for(i = 1; i < kpoints; i++)
706 void zwlk(
double* z,
double* w,
const int npt,
708 const double alpha,
const double beta)
722 fprintf (stderr,
"too few points in formula\n");
740 int kpoints = 2*np-1;
746 int ncoeffs = (int)ceil(3.0*np/2)-1;
752 double *a =
new double[ncoeffs+1];
754 double *b =
new double[ncoeffs+1];
760 for(i = 0; i < ncoeffs+1; i++)
780 double* a0 =
new double[ncoeffs];
782 double* b0 =
new double[ncoeffs];
786 chri1(ncoeffs,a,b,a0,b0,endl);
790 double* a1 =
new double[ncoeffs-1];
792 double* b1 =
new double[ncoeffs-1];
796 chri1(ncoeffs-1,a0,b0,a1,b1,endr);
802 double s = b1[0]/fabs(b1[0]);
810 double* z1 =
new double[2*np-3];
812 double* w1 =
new double[2*np-3];
814 for(i = 0; i < ncoeffs; i++)
828 double** zmatrix =
new double*[2*np-3];
830 for(i = 0; i < 2*np-3; i++)
834 zmatrix[i] =
new double[2*np-3];
836 for(j = 0; j < 2*np-3; j++)
846 for(i = 0; i < 2*np-3; i++)
858 TriQL(2*np-3, z1, w1, zmatrix);
864 double sumW1Z1 = 0.0;
866 for(i = 0; i < 2*np-3; i++)
870 w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
874 sumW1Z1 += z1[i]*w1[i];
880 double c0 = b[0]-sumW1;
882 double c1 = a[0]*b[0]-sumW1Z1;
890 w[0] = (c0*endr-c1)/(endr-endl);
892 w[2*np-2] = (c1-c0*endl)/(endr-endl);
896 for(i = 1; i < kpoints-1; i++)
932 void Dgj(
double *D,
const double *z,
const int np,
const double alpha,
940 double one = 1.0, two = 2.0;
958 pd = (
double *)malloc(np*
sizeof(
double));
960 jacobd(np,z,pd,np,alpha,beta);
964 for (i = 0; i < np; i++){
966 for (j = 0; j < np; j++){
972 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
976 D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
978 (two*(one - z[j]*z[j]));
1018 void Dgrjm(
double *D,
const double *z,
const int np,
const double alpha,
1036 double one = 1.0, two = 2.0;
1042 pd = (
double *)malloc(np*
sizeof(
double));
1046 pd[0] = pow(-one,np-1)*
gammaF(np+beta+one);
1050 jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
1052 for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
1056 for (i = 0; i < np; i++)
1058 for (j = 0; j < np; j++){
1062 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1068 D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
1074 D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
1076 (two*(one - z[j]*z[j]));
1116 void Dgrjp(
double *D,
const double *z,
const int np,
const double alpha,
1134 double one = 1.0, two = 2.0;
1140 pd = (
double *)malloc(np*
sizeof(
double));
1146 jacobd(np-1,z,pd,np-1,alpha+1,beta);
1148 for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
1150 pd[np-1] = -
gammaF(np+alpha+one);
1156 for (i = 0; i < np; i++)
1158 for (j = 0; j < np; j++){
1162 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1168 D[i*np+j] = (np + alpha + beta + one)*(np - one)/
1170 (two*(alpha + two));
1174 D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
1176 (two*(one - z[j]*z[j]));
1216 void Dglj(
double *D,
const double *z,
const int np,
const double alpha,
1234 double one = 1.0, two = 2.0;
1240 pd = (
double *)malloc(np*
sizeof(
double));
1244 pd[0] = two*pow(-one,np)*
gammaF(np + beta);
1248 jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
1250 for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
1252 pd[np-1] = -two*
gammaF(np + alpha);
1258 for (i = 0; i < np; i++)
1260 for (j = 0; j < np; j++){
1264 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1270 D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
1274 D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
1278 D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
1280 (two*(one - z[j]*z[j]));
1340 double hgj (
const int i,
const double z,
const double *zgj,
1342 const int np,
const double alpha,
const double beta)
1348 double zi, dz, p, pd, h;
1356 if (fabs(dz) <
EPS)
return 1.0;
1360 jacobd (1, &zi, &pd , np, alpha, beta);
1362 jacobfd(1, &z , &p, NULL , np, alpha, beta);
1416 double hgrjm (
const int i,
const double z,
const double *zgrj,
const int np,
1418 const double alpha,
const double beta)
1424 double zi, dz, p, pd, h;
1432 if (fabs(dz) <
EPS)
return 1.0;
1436 jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
1440 jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
1442 h = (1.0 + zi)*pd + p;
1444 jacobfd (1, &z, &p, NULL, np-1, alpha, beta + 1);
1446 h = (1.0 + z )*p/(h*dz);
1500 double hgrjp (
const int i,
const double z,
const double *zgrj,
const int np,
1502 const double alpha,
const double beta)
1508 double zi, dz, p, pd, h;
1516 if (fabs(dz) <
EPS)
return 1.0;
1520 jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
1524 jacobd (1, &zi, &pd, np-1, alpha+1, beta );
1526 h = (1.0 - zi)*pd - p;
1528 jacobfd (1, &z, &p, NULL, np-1, alpha+1, beta);
1530 h = (1.0 - z )*p/(h*dz);
1584 double hglj (
const int i,
const double z,
const double *zglj,
const int np,
1586 const double alpha,
const double beta)
1590 double one = 1., two = 2.;
1592 double zi, dz, p, pd, h;
1600 if (fabs(dz) <
EPS)
return 1.0;
1604 jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
1608 jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
1610 h = (one - zi*zi)*pd - two*zi*p;
1612 jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
1614 h = (one - z*z)*p/(h*dz);
1652 void Imgj(
double *im,
const double *zgj,
const double *zm,
const int nz,
1654 const int mz,
const double alpha,
const double beta){
1662 for (i = 0; i < nz; ++i) {
1664 for (j = 0; j < mz; ++j)
1670 im [i*mz+j] =
hgj(i, zp, zgj, nz, alpha, beta);
1710 void Imgrjm(
double *im,
const double *zgrj,
const double *zm,
const int nz,
1712 const int mz,
const double alpha,
const double beta)
1722 for (i = 0; i < nz; i++) {
1724 for (j = 0; j < mz; j++)
1730 im [i*mz+j] =
hgrjm(i, zp, zgrj, nz, alpha, beta);
1770 void Imgrjp(
double *im,
const double *zgrj,
const double *zm,
const int nz,
1772 const int mz,
const double alpha,
const double beta)
1782 for (i = 0; i < nz; i++) {
1784 for (j = 0; j < mz; j++)
1790 im [i*mz+j] =
hgrjp(i, zp, zgrj, nz, alpha, beta);
1832 void Imglj(
double *im,
const double *zglj,
const double *zm,
const int nz,
1834 const int mz,
const double alpha,
const double beta)
1844 for (i = 0; i < nz; i++) {
1846 for (j = 0; j < mz; j++)
1852 im[i*mz+j] =
hglj(i, zp, zglj, nz, alpha, beta);
1946 void jacobfd(
const int np,
const double *z,
double *poly_in,
double *polyd,
1948 const int n,
const double alpha,
const double beta){
1952 double zero = 0.0, one = 1.0, two = 2.0;
1966 for(i = 0; i < np; ++i)
1972 for(i = 0; i < np; ++i)
1982 for(i = 0; i < np; ++i)
1984 poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1988 for(i = 0; i < np; ++i)
1990 polyd[i] = 0.5*(alpha + beta + two);
2000 double two = 2.0, apb = alpha + beta;
2002 double *poly, *polyn1,*polyn2;
2008 polyn1 = (
double *)malloc(2*np*
sizeof(
double));
2018 polyn1 = (
double *)malloc(3*np*
sizeof(
double));
2028 for(i = 0; i < np; ++i){
2032 polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
2038 for(k = 2; k <= n; ++k){
2040 a1 = two*k*(k + apb)*(two*k + apb - two);
2042 a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
2044 a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
2046 a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
2058 for(i = 0; i < np; ++i){
2060 poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
2062 polyn2[i] = polyn1[i];
2064 polyn1[i] = poly [i];
2074 a1 = n*(alpha - beta);
2076 a2 = n*(two*n + alpha + beta);
2078 a3 = two*(n + alpha)*(n + beta);
2080 a4 = (two*n + alpha + beta);
2082 a1 /= a4; a2 /= a4; a3 /= a4;
2088 for(i = 0; i < np; ++i){
2090 polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
2092 polyd[i] /= (one - z[i]*z[i]);
2146 void jacobd(
const int np,
const double *z,
double *polyd,
const int n,
2148 const double alpha,
const double beta)
2158 for(i = 0; i < np; ++i) polyd[i] = 0.0;
2164 jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2166 for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (
double)n + one);
2206 if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2208 else if (!x)
return gamma;
2210 else if ((x-(
int)x) == 0.5){
2230 else if ((x-(
int)x) == 0.0){
2250 fprintf(stderr,
"%lf is not of integer or half order\n",x);
2274 static void Jacobz(
const int n,
double *z,
const double alpha,
2280 double dth = M_PI/(2.0*(double)n);
2282 double poly,pder,rlast=0.0;
2286 double one = 1.0, two = 2.0;
2296 for(k = 0; k < n; ++k){
2298 r = -cos((two*(
double)k + one) * dth);
2300 if(k) r = 0.5*(r + rlast);
2304 for(j = 1; j <
STOP; ++j){
2306 jacobfd(1,&r,&poly, &pder, n, alpha, beta);
2310 for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
2314 delr = -poly / (pder - sum * poly);
2318 if( fabs(delr) <
EPS )
break;
2384 void JacZeros(
const int n,
double *a,
double*b,
const double alpha,
2396 double **z =
new double*[n];
2398 for(i = 0; i < n; i++)
2402 z[i] =
new double[n];
2404 for(j = 0; j < n; j++)
2414 for(i = 0; i < n; i++)
2444 static void RecCoeff(
const int n,
double *a,
double *b,
const double alpha,
2452 double apb, apbi,a2b2;
2472 a[0] = (beta-alpha)/apbi;
2474 b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2478 a2b2 = beta*beta-alpha*alpha;
2482 for(i = 1; i < n-1; i++){
2484 apbi = 2.0*(i+1) + apb;
2486 a[i] = a2b2/((apbi-2.0)*apbi);
2488 b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2490 ((apbi*apbi-1)*apbi*apbi));
2498 a[n-1] = a2b2/((apbi-2.0)*apbi);
2560 static void TriQL(
const int n,
double *d,
double *e,
double **z){
2564 double s,r,p,g,f,dd,c,b;
2568 double MuZero = e[0];
2574 for(i = 0; i < n-1; i++)
2578 e[i] = sqrt(e[i+1]);
2594 for (m=l;m<n-1;m++) {
2596 dd=fabs(d[m])+fabs(d[m+1]);
2598 if (fabs(e[m])+dd == dd)
break;
2604 if (iter++ ==
STOP){
2606 fprintf(stderr,
"triQL: Too many iterations in TQLI");
2612 g=(d[l+1]-d[l])/(2.0*e[l]);
2616 g=d[m]-d[l]+e[l]/(g+
sign(r,g));
2622 for (i=m-1;i>=l;i--) {
2628 if (fabs(f) >= fabs(g)) {
2652 r=(d[i]-g)*s+2.0*c*b;
2664 for(k = 0; k < n; k++)
2670 z[k][i+1] = s*z[k][i] + c*f;
2672 z[k][i] = c*z[k][i] - s*f;
2700 for(i = 0; i < n-1; ++i){
2706 for(l = i+1; l < n; ++l)
2722 double temp = z[0][k];
2734 for(i =0 ; i < n; i++)
2738 e[i] = MuZero*z[0][i]*z[0][i];
2772 int size = (int)floor(n/2.0)+2;
2774 double *s =
new double[size];
2776 double *t =
new double[size];
2782 for(i = 0; i < size; i++)
2796 for(m = 0; m <= n-2; m++)
2802 for(k = (
int)floor((m+1)/2.0); k >= 0; k--)
2808 u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
2830 for(j = (
int)floor(n/2.0); j >= 0; j--)
2840 for(m = n-1; m <= 2*n-3; m++)
2846 for(k = m+1-n; k <= floor((m-1)/2.0); k++)
2854 u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
2868 a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
2878 b[k+n+1] = s[j+1]/s[j+2];
2898 a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
2926 void chri1(
int n,
double* a,
double* b,
double* a0,
2928 double* b0,
double z)
2934 double q = ceil(3.0*n/2);
2936 int size = (int)q+1;
2942 fprintf(stderr,
"input arrays a and b are too short\n");
2946 double* r =
new double[n+1];
2950 r[1] = z - a[1] - b[1]/r[0];
2952 a0[0] = a[1] + r[1] - r[0];
2968 for(k = 1; k < n; k++)
2972 r[k+1] = z - a[k+1] - b[k+1]/r[k];
2974 a0[k] = a[k+1] + r[k+1] - r[k];
2976 b0[k] = b[k] * r[k]/r[k-1];
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
#define sign(a, b)
return the sign(b)*a
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through.
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
void zwlk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
void JacZeros(const int n, double *a, double *b, const double alpha, const double beta)
Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal.
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points.
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
The namespace associated with the the Polylib library (Polylib introduction)
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
static void Jacobz(const int n, double *z, const double alpha, const double beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e.
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Radau-Kronrod-Jacobi zeros and weights.
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
#define EPS
Precision tolerance for two points to be similar.
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an.
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
double gammaF(const double)
Calculate the Gamma function , , for integer.
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.