Nektar++
Macros | Functions
Polylib_test.cpp File Reference
#include "Polylib.h"
#include <cstdio>
#include <cmath>

Go to the source code of this file.

Macros

#define NPLOWER   5
 
#define NPUPPER   15
 
#define EPS   1e-12
 
#define GAUSS_INT   1
 
#define GAUSS_RADAUM_INT   1
 
#define GAUSS_RADAUP_INT   1
 
#define GAUSS_LOBATTO_INT   1
 
#define GAUSS_DIFF   1
 
#define GAUSS_RADAUM_DIFF   1
 
#define GAUSS_RADAUP_DIFF   1
 
#define GAUSS_LOBATTO_DIFF   1
 
#define GAUSS_INTERP   1
 
#define GAUSS_RADAUM_INTERP   1
 
#define GAUSS_RADAUP_INTERP   1
 
#define GAUSS_LOBATTO_INTERP   1
 

Functions

double ddot (int, double *, int, double *, int)
 
double * dvector (int, int)
 
int test_gamma_fraction ()
 
int main ()
 

Macro Definition Documentation

◆ EPS

#define EPS   1e-12

Definition at line 96 of file Polylib_test.cpp.

◆ GAUSS_DIFF

#define GAUSS_DIFF   1

Definition at line 102 of file Polylib_test.cpp.

◆ GAUSS_INT

#define GAUSS_INT   1

Definition at line 98 of file Polylib_test.cpp.

◆ GAUSS_INTERP

#define GAUSS_INTERP   1

Definition at line 106 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_DIFF

#define GAUSS_LOBATTO_DIFF   1

Definition at line 105 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_INT

#define GAUSS_LOBATTO_INT   1

Definition at line 101 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_INTERP

#define GAUSS_LOBATTO_INTERP   1

Definition at line 109 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_DIFF

#define GAUSS_RADAUM_DIFF   1

Definition at line 103 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_INT

#define GAUSS_RADAUM_INT   1

Definition at line 99 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_INTERP

#define GAUSS_RADAUM_INTERP   1

Definition at line 107 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_DIFF

#define GAUSS_RADAUP_DIFF   1

Definition at line 104 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_INT

#define GAUSS_RADAUP_INT   1

Definition at line 100 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_INTERP

#define GAUSS_RADAUP_INTERP   1

Definition at line 108 of file Polylib_test.cpp.

◆ NPLOWER

#define NPLOWER   5

Definition at line 94 of file Polylib_test.cpp.

◆ NPUPPER

#define NPUPPER   15

Definition at line 95 of file Polylib_test.cpp.

Function Documentation

◆ ddot()

double ddot ( int  n,
double *  x,
int  incx,
double *  y,
int  incy 
)

Definition at line 485 of file Polylib_test.cpp.

486 {
487  register double sum = 0.;
488 
489  while (n--) {
490  sum += (*x) * (*y);
491  x += incx;
492  y += incy;
493  }
494  return sum;
495 }

Referenced by main().

◆ dvector()

double * dvector ( int  nl,
int  nh 
)

Definition at line 498 of file Polylib_test.cpp.

499 {
500  double *v;
501 
502  v = (double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
503  return v-nl;
504 }

Referenced by main().

◆ main()

int main ( )

Definition at line 116 of file Polylib_test.cpp.

116  {
117  int np,n,i;
118  double *z,*w,*p,sum=0,alpha,beta,*d,*dt;
119 
120  z = dvector(0,NPUPPER-1);
121  w = dvector(0,NPUPPER-1);
122  p = dvector(0,NPUPPER-1);
123 
124  d = dvector(0,NPUPPER*NPUPPER-1);
125  dt = dvector(0,NPUPPER*NPUPPER-1);
126 
127 #if GAUSS_INT
128  /* Gauss Integration */
129  printf("Begin checking Gauss integration\n");
130  alpha = -0.5;
131  while(alpha <= 5.0){
132  beta = -0.5;
133  while(beta <= 5.0){
134 
135  for(np = NPLOWER; np <= NPUPPER; ++np){
136  zwgj(z,w,np,alpha,beta);
137  for(n = 2; n < 2*np-1; ++n){
138  jacobfd(np,z,p,NULL,n,alpha,beta);
139  sum = ddot(np,w,1,p,1);
140  if(fabs(sum)>EPS)
141  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
142  ,alpha,beta,np,n,sum);
143  }
144  }
145 
146  beta += 0.5;
147  }
148  printf("finished checking all beta values for alpha = %lf\n",alpha);
149  alpha += 0.5;
150  }
151  printf("Finished checking Gauss Integration\n");
152 #endif
153 
154 #if GAUSS_RADAUM_INT
155  /* Gauss Radau Integration */
156  printf("Begin checking Gauss Radau Integration\n");
157  alpha = -0.5;
158  while(alpha <= 5.0){
159  beta = -0.5;
160  while(beta <= 5.0){
161  for(np = NPLOWER; np <= NPUPPER; ++np){
162  zwgrjm(z,w,np,alpha,beta);
163  for(n = 2; n < 2*np-2; ++n){
164  jacobfd(np,z,p,NULL,n,alpha,beta);
165  sum = ddot(np,w,1,p,1);
166  if(fabs(sum)>EPS)
167  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
168  ,alpha,beta,np,n,sum);
169  }
170  }
171 
172  beta += 0.5;
173  }
174  printf("finished checking all beta values for alpha = %lf\n",alpha);
175  alpha += 0.5;
176  }
177  printf("Finished checking Gauss Radau (z=-1) Integration\n");
178 #endif
179 
180 
181 #if GAUSS_RADAUP_INT
182  /* Gauss Radau Integration */
183  printf("Begin checking Gauss Radau Integration\n");
184  alpha = -0.5;
185  while(alpha <= 5.0){
186  beta = -0.5;
187  while(beta <= 5.0){
188  for(np = NPLOWER; np <= NPUPPER; ++np){
189  zwgrjp(z,w,np,alpha,beta);
190  for(n = 2; n < 2*np-2; ++n){
191  jacobfd(np,z,p,NULL,n,alpha,beta);
192  sum = ddot(np,w,1,p,1);
193  if(fabs(sum)>EPS)
194  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
195  ,alpha,beta,np,n,sum);
196  }
197  }
198 
199  beta += 0.5;
200  }
201  printf("finished checking all beta values for alpha = %lf\n",alpha);
202  alpha += 0.5;
203  }
204  printf("Finished checking Gauss Radau (z=1) Integration\n");
205 #endif
206 
207 #if GAUSS_LOBATTO_INT
208  /* Gauss Lobatto Integration */
209  printf("Begin checking Gauss Lobatto integration\n");
210  alpha = -0.5;
211  while(alpha <= 5.0){
212  beta = -0.5;
213  while(beta <= 5.0){
214 
215  for(np = NPLOWER; np <= NPUPPER; ++np){
216  zwglj(z,w,np,alpha,beta);
217  for(n = 2; n < 2*np-3; ++n){
218  jacobfd(np,z,p,NULL,n,alpha,beta);
219  sum = ddot(np,w,1,p,1);
220  if(fabs(sum)>EPS)
221  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
222  ,alpha,beta,np,n,sum);
223  }
224  }
225 
226  beta += 0.5;
227  }
228  printf("finished checking all beta values for alpha = %lf\n",alpha);
229  alpha += 0.5;
230  }
231  printf("Finished checking Gauss Lobatto Integration\n");
232 #endif
233 
234 #if GAUSS_DIFF
235  printf("Begin checking differentiation through Gauss points\n");
236  alpha = -0.5;
237  while(alpha <= 5.0){
238  beta = -0.5;
239  while(beta <= 5.0){
240 
241  for(np = NPLOWER; np <= NPUPPER; ++np){
242  zwgj(z,w,np,alpha,beta);
243  for(n = 2; n < np-1; ++n){
244  Dgj(d,z,np,alpha,beta);
245 
246  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
247  sum = 0;
248  for(i = 0; i < np; ++i)
249  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
250  sum /= np;
251  if(fabs(sum)>EPS)
252  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
253  ,alpha,beta,np,n,sum);
254  }
255  }
256  beta += 0.5;
257  }
258  printf("finished checking all beta values for alpha = %lf\n",alpha);
259  alpha += 0.5;
260  }
261  printf("Finished checking Gauss Jacobi differentiation\n");
262 #endif
263 
264 #if GAUSS_RADAUM_DIFF
265  printf("Begin checking differentiation through Gauss Radau points\n");
266  alpha = -0.5;
267  while(alpha <= 5.0){
268  beta = -0.5;
269  while(beta <= 5.0){
270 
271  for(np = NPLOWER; np <= NPUPPER; ++np){
272  zwgrjm(z,w,np,alpha,beta);
273  for(n = 2; n < np-1; ++n){
274  Dgrjm(d,z,np,alpha,beta);
275  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
276  sum = 0;
277  for(i = 0; i < np; ++i)
278  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
279  sum /= np;
280  if(fabs(sum)>EPS)
281  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
282  ,alpha,beta,np,n,sum);
283  }
284  }
285  beta += 0.5;
286  }
287  printf("finished checking all beta values for alpha = %lf\n",alpha);
288  alpha += 0.5;
289  }
290  printf("Finished checking Gauss Radau (z=-1) differentiation\n");
291 #endif
292 
293 #if GAUSS_RADAUP_DIFF
294  printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
295  alpha = -0.5;
296  while(alpha <= 5.0){
297  beta = -0.5;
298  while(beta <= 5.0){
299 
300  for(np = NPLOWER; np <= NPUPPER; ++np){
301  zwgrjp(z,w,np,alpha,beta);
302  for(n = 2; n < np-1; ++n){
303  Dgrjp(d,z,np,alpha,beta);
304  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
305  sum = 0;
306  for(i = 0; i < np; ++i)
307  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
308  sum /= np;
309  if(fabs(sum)>EPS)
310  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
311  ,alpha,beta,np,n,sum);
312  }
313  }
314  beta += 0.5;
315  }
316  printf("finished checking all beta values for alpha = %lf\n",alpha);
317  alpha += 0.5;
318  }
319  printf("Finished checking Gauss Radau (z=1) differentiation\n");
320 #endif
321 
322 #if GAUSS_LOBATTO_DIFF
323  printf("Begin checking differentiation through Gauss Lobatto points\n");
324  alpha = -0.5;
325  while(alpha <= 5.0){
326  beta = -0.5;
327  while(beta <= 5.0){
328 
329  for(np = NPLOWER; np <= NPUPPER; ++np){
330  zwglj(z,w,np,alpha,beta);
331  for(n = 2; n < np-1; ++n){
332  Dglj(d,z,np,alpha,beta);
333  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
334  sum = 0;
335  for(i = 0; i < np; ++i)
336  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
337  sum /= np;
338  if(fabs(sum)>EPS)
339  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
340  ,alpha,beta,np,n,sum);
341  }
342  }
343  beta += 0.5;
344  }
345  printf("finished checking all beta values for alpha = %lf\n",alpha);
346  alpha += 0.5;
347  }
348  printf("Finished checking Gauss Lobatto differentiation\n");
349 #endif
350 
351  /* check interpolation routines */
352 #if GAUSS_INTERP
353  printf("Begin checking interpolation through Gauss points\n");
354  alpha = -0.5;
355  while(alpha <= 5.0){
356  beta = -0.5;
357  while(beta <= 5.0){
358 
359  for(np = NPLOWER; np <= NPUPPER; ++np){
360  zwgj(z,w,np,alpha,beta);
361  for(n = 2; n < np-1; ++n){
362  for(i = 0; i < np; ++i) {
363  w[i] = 2.0*i/(double)(np-1)-1.0;
364  p[i] = pow(z[i],n);
365  }
366  Imgj(d,z,w,np,np,alpha,beta);
367  sum = 0;
368  for(i = 0; i < np; ++i)
369  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
370  sum /= np;
371  if(fabs(sum)>EPS)
372  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
373  ,alpha,beta,np,n,sum);
374  }
375  }
376  beta += 0.5;
377  }
378  printf("finished checking all beta values for alpha = %lf\n",alpha);
379  alpha += 0.5;
380  }
381  printf("Finished checking Gauss Jacobi interpolation\n");
382 #endif
383 
384 #if GAUSS_RADAUM_INTERP
385  printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
386  alpha = -0.5;
387  while(alpha <= 5.0){
388  beta = -0.5;
389  while(beta <= 5.0){
390 
391  for(np = NPLOWER; np <= NPUPPER; ++np){
392  zwgrjm(z,w,np,alpha,beta);
393  for(n = 2; n < np-1; ++n){
394  for(i = 0; i < np; ++i) {
395  w[i] = 2.0*i/(double)(np-1)-1.0;
396  p[i] = pow(z[i],n);
397  }
398  Imgrjm(d,z,w,np,np,alpha,beta);
399  sum = 0;
400  for(i = 0; i < np; ++i)
401  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
402  sum /= np;
403  if(fabs(sum)>EPS)
404  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
405  ,alpha,beta,np,n,sum);
406  }
407  }
408  beta += 0.5;
409  }
410  printf("finished checking all beta values for alpha = %lf\n",alpha);
411  alpha += 0.5;
412  }
413  printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
414 #endif
415 #if GAUSS_RADAUP_INTERP
416  printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
417  alpha = -0.5;
418  while(alpha <= 5.0){
419  beta = -0.5;
420  while(beta <= 5.0){
421 
422  for(np = NPLOWER; np <= NPUPPER; ++np){
423  zwgrjp(z,w,np,alpha,beta);
424  for(n = 2; n < np-1; ++n){
425  for(i = 0; i < np; ++i) {
426  w[i] = 2.0*i/(double)(np-1)-1.0;
427  p[i] = pow(z[i],n);
428  }
429  Imgrjp(d,z,w,np,np,alpha,beta);
430  sum = 0;
431  for(i = 0; i < np; ++i)
432  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
433  sum /= np;
434  if(fabs(sum)>EPS)
435  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
436  ,alpha,beta,np,n,sum);
437  }
438  }
439  beta += 0.5;
440  }
441  printf("finished checking all beta values for alpha = %lf\n",alpha);
442  alpha += 0.5;
443  }
444  printf("Finished checking Gauss Radau (z=1) interpolation\n");
445 #endif
446 
447 #if GAUSS_LOBATTO_INTERP
448  printf("Begin checking Interpolation through Gauss Lobatto points\n");
449  alpha = -0.5;
450  while(alpha <= 5.0){
451  beta = -0.5;
452  while(beta <= 5.0){
453 
454  for(np = NPLOWER; np <= NPUPPER; ++np){
455  zwglj(z,w,np,alpha,beta);
456  for(n = 2; n < np-1; ++n){
457  for(i = 0; i < np; ++i) {
458  w[i] = 2.0*i/(double)(np-1)-1.0;
459  p[i] = pow(z[i],n);
460  }
461  Imglj(d,z,w,np,np,alpha,beta);
462  sum = 0;
463  for(i = 0; i < np; ++i)
464  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
465  sum /= np;
466  if(fabs(sum)>EPS)
467  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
468  ,alpha,beta,np,n,sum);
469  }
470  }
471  beta += 0.5;
472  }
473  printf("finished checking all beta values for alpha = %lf\n",alpha);
474  alpha += 0.5;
475  }
476  printf("Finished checking Gauss Lobatto interploation\n");
477 #endif
478  test_gamma_fraction();
479 
480 
481  free(z); free(w); free(p); free(d); free(dt);
482  return 0;
483 }
dvector
double * dvector(int, int)
Definition: Polylib_test.cpp:498
NPUPPER
#define NPUPPER
Definition: Polylib_test.cpp:95
test_gamma_fraction
int test_gamma_fraction()
Definition: Polylib_test.cpp:506
EPS
#define EPS
Definition: Polylib_test.cpp:96
ddot
double ddot(int, double *, int, double *, int)
Definition: Polylib_test.cpp:485
NPLOWER
#define NPLOWER
Definition: Polylib_test.cpp:94
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Polylib::Dgj
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
Definition: Polylib.cpp:546
Polylib::zwgrjm
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.
Definition: Polylib.cpp:119
Polylib::zwgrjp
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.
Definition: Polylib.cpp:158
Polylib::Dglj
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
Definition: Polylib.cpp:690
Polylib::Imgrjp
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 (including z=1) to an arbitrary distrubtion at ...
Definition: Polylib.cpp:946
Polylib::zwglj
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.
Definition: Polylib.cpp:195
Polylib::zwgj
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:91
Polylib::Imgj
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 arbitrary distribution at points zm.
Definition: Polylib.cpp:887
Polylib::Imgrjm
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 (including z=-1) to an arbitrary distrubtion at...
Definition: Polylib.cpp:916
Polylib::Dgrjm
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
Definition: Polylib.cpp:589
Polylib::Dgrjp
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.
Definition: Polylib.cpp:639
Polylib::Imglj
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 to an arbitrary distrubtion at points zm.
Definition: Polylib.cpp:977
Polylib::jacobfd
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, .
Definition: Polylib.cpp:1034

References ddot(), Polylib::Dgj(), Polylib::Dglj(), Polylib::Dgrjm(), Polylib::Dgrjp(), dvector(), EPS, Polylib::Imgj(), Polylib::Imglj(), Polylib::Imgrjm(), Polylib::Imgrjp(), Polylib::jacobfd(), NPLOWER, NPUPPER, CellMLToNektar.cellml_metadata::p, test_gamma_fraction(), Polylib::zwgj(), Polylib::zwglj(), Polylib::zwgrjm(), and Polylib::zwgrjp().

◆ test_gamma_fraction()

int test_gamma_fraction ( )

Definition at line 506 of file Polylib_test.cpp.

506  {
507  double a = 362880.;
508  double c = Polylib::gammaFracGammaF(10, 0., 0, 1.);
509  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.,1., c/a-1.);
510  c = Polylib::gammaFracGammaF(1, 0., 12, -2.);
511  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 1.,10., c*a-1.);
512 
513  a = 30.;
514  c = Polylib::gammaFracGammaF(4, 3., 5, 0.);
515  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.,7., c/a-1.);
516  c = Polylib::gammaFracGammaF(5, 0., 7, 0.);
517  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 7.,5., c*a-1.);
518 
519  a = 21651.09375;
520  c = Polylib::gammaFracGammaF(7, 3.5, 5, 0.5);
521  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,5.5, c/a-1.);
522  c = Polylib::gammaFracGammaF(5, 0.5, 10, 0.5);
523  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5,10.5, c*a-1.);
524 
525  a = 97429.921875;
526  c = Polylib::gammaFracGammaF(10, 0.5, 5, -0.5);
527  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,4.5, c/a-1.);
528  c = Polylib::gammaFracGammaF(5, -0.5, 11, -0.5);
529  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 4.5,10.5, c*a-1.);
530 
531  a = 2279.0625;
532  c = Polylib::gammaFracGammaF(10, -0.5, 5, 0.5);
533  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 9.5,5.5, c/a-1.);
534  c = Polylib::gammaFracGammaF(5, 0.5, 10, -0.5);
535  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5,9.5, c*a-1.);
536 
537  a = 639383.8623046875;
538  c = Polylib::gammaFracGammaF(10, 0.5, 0, 0.5);
539  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,0.5, c/a-1.);
540  c = Polylib::gammaFracGammaF(0, 0.5, 10, 0.5);
541  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 0.5,10.5, c*a-1.);
542 
543  a = 5967498288235982848.;
544  c = Polylib::gammaFracGammaF(100, 0., 90, 0.5);
545  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 100.,90.5, c/a-1.);
546  c = Polylib::gammaFracGammaF(90, 0.5, 100, 0.);
547  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 90.5,100., c*a-1.);
548 
549  a = 5618641603777298169856.;
550  c = Polylib::gammaFracGammaF(200, 0., 191, -0.5);
551  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200.,190.5, c/a-1.);
552  c = Polylib::gammaFracGammaF(190, 0.5, 200, 0.);
553  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.5,200., c*a-1.);
554 
555  a = 77396694214720029196288.;
556  c = Polylib::gammaFracGammaF(200, 0., 190, 0.);
557  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200.,190., c/a-1.);
558  c = Polylib::gammaFracGammaF(190, 0., 200, 0.);
559  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.,200., c*a-1.);
560 
561  int Q=2;
562  for(double alpha=-0.5; alpha<=150.; alpha+=0.5)
563  {
564  for(double beta=-0.5; beta<=150.; beta+=0.5)
565  {
566  a = Polylib::gammaF(Q+alpha)/Polylib::gammaF(Q+beta);
567  c = Polylib::gammaFracGammaF(Q, alpha, Q, beta) / a - 1.;
568  if(fabs(c)>5.e-15)
569  {
570  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", Q+alpha,Q+beta, c);
571  }
572  }
573  }
574 
575  return 0;
576 }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1161
Polylib::gammaFracGammaF
double gammaFracGammaF(const int, const double, const int, const double)
Calculate fraction of two Gamma functions, , for integer values and halves.
Definition: Polylib.cpp:1203

References Polylib::gammaF(), and Polylib::gammaFracGammaF().

Referenced by main().

Generated on Wed Nov 24 2021 10:33:30 for Nektar++ by doxygen 1.9.1