diff options
| author | Joel Sherrill <joel.sherrill@OARcorp.com> | 2011-04-08 17:33:11 +0000 |
|---|---|---|
| committer | Joel Sherrill <joel.sherrill@OARcorp.com> | 2011-04-08 17:33:11 +0000 |
| commit | 73f643f3f4a55310b2c8c1a9858906b2dd676e72 (patch) | |
| tree | b1df97f18dace4a5702b0bc7aafdfee8a2f25ada /gsl-1.9/specfunc/gamma.c | |
| parent | 2011-04-08 Joel Sherrill <joel.sherrill@oarcorp.com> (diff) | |
| download | rtems-addon-packages-73f643f3f4a55310b2c8c1a9858906b2dd676e72.tar.bz2 | |
2011-04-08 Joel Sherrill <joel.sherrill@oarcorp.com>
* AUTHORS, BUGS, COPYING, ChangeLog, INSTALL, Makefile.am, Makefile.in,
NEWS, README, SUPPORT, THANKS, TODO, acconfig.h, aclocal.m4,
autogen.sh, config.guess, config.h.in, config.sub, configure,
configure.ac, gsl-config.in, gsl-histogram.c, gsl-randist.c, gsl.m4,
gsl.pc.in, gsl.spec.in, gsl_machine.h, gsl_math.h, gsl_mode.h,
gsl_nan.h, gsl_pow_int.h, gsl_precision.h, gsl_types.h,
gsl_version.h.in, install-sh, ltmain.sh, mdate-sh, missing,
mkinstalldirs, templates_off.h, templates_on.h,
test_gsl_histogram.sh, version.c, blas/ChangeLog, blas/Makefile.am,
blas/Makefile.in, blas/TODO, blas/blas.c, blas/gsl_blas.h,
blas/gsl_blas_types.h, block/ChangeLog, block/Makefile.am,
block/Makefile.in, block/block.c, block/block_source.c, block/file.c,
block/fprintf_source.c, block/fwrite_source.c, block/gsl_block.h,
block/gsl_block_char.h, block/gsl_block_complex_double.h,
block/gsl_block_complex_float.h,
block/gsl_block_complex_long_double.h, block/gsl_block_double.h,
block/gsl_block_float.h, block/gsl_block_int.h,
block/gsl_block_long.h, block/gsl_block_long_double.h,
block/gsl_block_short.h, block/gsl_block_uchar.h,
block/gsl_block_uint.h, block/gsl_block_ulong.h,
block/gsl_block_ushort.h, block/gsl_check_range.h, block/init.c,
block/init_source.c, block/test.c, block/test_complex_io.c,
block/test_complex_source.c, block/test_io.c, block/test_source.c,
bspline/ChangeLog, bspline/Makefile.am, bspline/Makefile.in,
bspline/TODO, bspline/bspline.c, bspline/gsl_bspline.h,
bspline/test.c, cblas/ChangeLog, cblas/Makefile.am,
cblas/Makefile.in, cblas/TODO, cblas/caxpy.c, cblas/cblas.h,
cblas/ccopy.c, cblas/cdotc_sub.c, cblas/cdotu_sub.c, cblas/cgbmv.c,
cblas/cgemm.c, cblas/cgemv.c, cblas/cgerc.c, cblas/cgeru.c,
cblas/chbmv.c, cblas/chemm.c, cblas/chemv.c, cblas/cher.c,
cblas/cher2.c, cblas/cher2k.c, cblas/cherk.c, cblas/chpmv.c,
cblas/chpr.c, cblas/chpr2.c, cblas/cscal.c, cblas/csscal.c,
cblas/cswap.c, cblas/csymm.c, cblas/csyr2k.c, cblas/csyrk.c,
cblas/ctbmv.c, cblas/ctbsv.c, cblas/ctpmv.c, cblas/ctpsv.c,
cblas/ctrmm.c, cblas/ctrmv.c, cblas/ctrsm.c, cblas/ctrsv.c,
cblas/dasum.c, cblas/daxpy.c, cblas/dcopy.c, cblas/ddot.c,
cblas/dgbmv.c, cblas/dgemm.c, cblas/dgemv.c, cblas/dger.c,
cblas/dnrm2.c, cblas/drot.c, cblas/drotg.c, cblas/drotm.c,
cblas/drotmg.c, cblas/dsbmv.c, cblas/dscal.c, cblas/dsdot.c,
cblas/dspmv.c, cblas/dspr.c, cblas/dspr2.c, cblas/dswap.c,
cblas/dsymm.c, cblas/dsymv.c, cblas/dsyr.c, cblas/dsyr2.c,
cblas/dsyr2k.c, cblas/dsyrk.c, cblas/dtbmv.c, cblas/dtbsv.c,
cblas/dtpmv.c, cblas/dtpsv.c, cblas/dtrmm.c, cblas/dtrmv.c,
cblas/dtrsm.c, cblas/dtrsv.c, cblas/dzasum.c, cblas/dznrm2.c,
cblas/gsl_cblas.h, cblas/hypot.c, cblas/icamax.c, cblas/idamax.c,
cblas/isamax.c, cblas/izamax.c, cblas/sasum.c, cblas/saxpy.c,
cblas/scasum.c, cblas/scnrm2.c, cblas/scopy.c, cblas/sdot.c,
cblas/sdsdot.c, cblas/sgbmv.c, cblas/sgemm.c, cblas/sgemv.c,
cblas/sger.c, cblas/snrm2.c, cblas/source_asum_c.h,
cblas/source_asum_r.h, cblas/source_axpy_c.h, cblas/source_axpy_r.h,
cblas/source_copy_c.h, cblas/source_copy_r.h, cblas/source_dot_c.h,
cblas/source_dot_r.h, cblas/source_gbmv_c.h, cblas/source_gbmv_r.h,
cblas/source_gemm_c.h, cblas/source_gemm_r.h, cblas/source_gemv_c.h,
cblas/source_gemv_r.h, cblas/source_ger.h, cblas/source_gerc.h,
cblas/source_geru.h, cblas/source_hbmv.h, cblas/source_hemm.h,
cblas/source_hemv.h, cblas/source_her.h, cblas/source_her2.h,
cblas/source_her2k.h, cblas/source_herk.h, cblas/source_hpmv.h,
cblas/source_hpr.h, cblas/source_hpr2.h, cblas/source_iamax_c.h,
cblas/source_iamax_r.h, cblas/source_nrm2_c.h, cblas/source_nrm2_r.h,
cblas/source_rot.h, cblas/source_rotg.h, cblas/source_rotm.h,
cblas/source_rotmg.h, cblas/source_sbmv.h, cblas/source_scal_c.h,
cblas/source_scal_c_s.h, cblas/source_scal_r.h, cblas/source_spmv.h,
cblas/source_spr.h, cblas/source_spr2.h, cblas/source_swap_c.h,
cblas/source_swap_r.h, cblas/source_symm_c.h, cblas/source_symm_r.h,
cblas/source_symv.h, cblas/source_syr.h, cblas/source_syr2.h,
cblas/source_syr2k_c.h, cblas/source_syr2k_r.h,
cblas/source_syrk_c.h, cblas/source_syrk_r.h, cblas/source_tbmv_c.h,
cblas/source_tbmv_r.h, cblas/source_tbsv_c.h, cblas/source_tbsv_r.h,
cblas/source_tpmv_c.h, cblas/source_tpmv_r.h, cblas/source_tpsv_c.h,
cblas/source_tpsv_r.h, cblas/source_trmm_c.h, cblas/source_trmm_r.h,
cblas/source_trmv_c.h, cblas/source_trmv_r.h, cblas/source_trsm_c.h,
cblas/source_trsm_r.h, cblas/source_trsv_c.h, cblas/source_trsv_r.h,
cblas/srot.c, cblas/srotg.c, cblas/srotm.c, cblas/srotmg.c,
cblas/ssbmv.c, cblas/sscal.c, cblas/sspmv.c, cblas/sspr.c,
cblas/sspr2.c, cblas/sswap.c, cblas/ssymm.c, cblas/ssymv.c,
cblas/ssyr.c, cblas/ssyr2.c, cblas/ssyr2k.c, cblas/ssyrk.c,
cblas/stbmv.c, cblas/stbsv.c, cblas/stpmv.c, cblas/stpsv.c,
cblas/strmm.c, cblas/strmv.c, cblas/strsm.c, cblas/strsv.c,
cblas/test.c, cblas/test_amax.c, cblas/test_asum.c,
cblas/test_axpy.c, cblas/test_copy.c, cblas/test_dot.c,
cblas/test_gbmv.c, cblas/test_gemm.c, cblas/test_gemv.c,
cblas/test_ger.c, cblas/test_hbmv.c, cblas/test_hemm.c,
cblas/test_hemv.c, cblas/test_her.c, cblas/test_her2.c,
cblas/test_her2k.c, cblas/test_herk.c, cblas/test_hpmv.c,
cblas/test_hpr.c, cblas/test_hpr2.c, cblas/test_nrm2.c,
cblas/test_rot.c, cblas/test_rotg.c, cblas/test_rotm.c,
cblas/test_rotmg.c, cblas/test_sbmv.c, cblas/test_scal.c,
cblas/test_spmv.c, cblas/test_spr.c, cblas/test_spr2.c,
cblas/test_swap.c, cblas/test_symm.c, cblas/test_symv.c,
cblas/test_syr.c, cblas/test_syr2.c, cblas/test_syr2k.c,
cblas/test_syrk.c, cblas/test_tbmv.c, cblas/test_tbsv.c,
cblas/test_tpmv.c, cblas/test_tpsv.c, cblas/test_trmm.c,
cblas/test_trmv.c, cblas/test_trsm.c, cblas/test_trsv.c,
cblas/tests.c, cblas/tests.h, cblas/xerbla.c, cblas/zaxpy.c,
cblas/zcopy.c, cblas/zdotc_sub.c, cblas/zdotu_sub.c, cblas/zdscal.c,
cblas/zgbmv.c, cblas/zgemm.c, cblas/zgemv.c, cblas/zgerc.c,
cblas/zgeru.c, cblas/zhbmv.c, cblas/zhemm.c, cblas/zhemv.c,
cblas/zher.c, cblas/zher2.c, cblas/zher2k.c, cblas/zherk.c,
cblas/zhpmv.c, cblas/zhpr.c, cblas/zhpr2.c, cblas/zscal.c,
cblas/zswap.c, cblas/zsymm.c, cblas/zsyr2k.c, cblas/zsyrk.c,
cblas/ztbmv.c, cblas/ztbsv.c, cblas/ztpmv.c, cblas/ztpsv.c,
cblas/ztrmm.c, cblas/ztrmv.c, cblas/ztrsm.c, cblas/ztrsv.c,
cdf/ChangeLog, cdf/Makefile.am, cdf/Makefile.in, cdf/beta.c,
cdf/beta_inc.c, cdf/betainv.c, cdf/binomial.c, cdf/cauchy.c,
cdf/cauchyinv.c, cdf/chisq.c, cdf/chisqinv.c, cdf/error.h,
cdf/exponential.c, cdf/exponentialinv.c, cdf/exppow.c, cdf/fdist.c,
cdf/fdistinv.c, cdf/flat.c, cdf/flatinv.c, cdf/gamma.c,
cdf/gammainv.c, cdf/gauss.c, cdf/gaussinv.c, cdf/geometric.c,
cdf/gsl_cdf.h, cdf/gumbel1.c, cdf/gumbel1inv.c, cdf/gumbel2.c,
cdf/gumbel2inv.c, cdf/hypergeometric.c, cdf/laplace.c,
cdf/laplaceinv.c, cdf/logistic.c, cdf/logisticinv.c, cdf/lognormal.c,
cdf/lognormalinv.c, cdf/nbinomial.c, cdf/pareto.c, cdf/paretoinv.c,
cdf/pascal.c, cdf/poisson.c, cdf/rat_eval.h, cdf/rayleigh.c,
cdf/rayleighinv.c, cdf/tdist.c, cdf/tdistinv.c, cdf/test.c,
cdf/test_auto.c, cdf/weibull.c, cdf/weibullinv.c, cheb/ChangeLog,
cheb/Makefile.am, cheb/Makefile.in, cheb/deriv.c, cheb/eval.c,
cheb/gsl_chebyshev.h, cheb/init.c, cheb/integ.c, cheb/test.c,
combination/ChangeLog, combination/Makefile.am,
combination/Makefile.in, combination/combination.c,
combination/file.c, combination/gsl_combination.h,
combination/init.c, combination/test.c, complex/ChangeLog,
complex/Makefile.am, complex/Makefile.in, complex/TODO,
complex/gsl_complex.h, complex/gsl_complex_math.h, complex/math.c,
complex/results.h, complex/results1.h, complex/results_real.h,
complex/test.c, const/ChangeLog, const/Makefile.am,
const/Makefile.in, const/TODO, const/gsl_const.h,
const/gsl_const_cgs.h, const/gsl_const_cgsm.h, const/gsl_const_mks.h,
const/gsl_const_mksa.h, const/gsl_const_num.h, const/test.c,
deriv/ChangeLog, deriv/Makefile.am, deriv/Makefile.in, deriv/deriv.c,
deriv/gsl_deriv.h, deriv/test.c, dht/ChangeLog, dht/Makefile.am,
dht/Makefile.in, dht/dht.c, dht/gsl_dht.h, dht/test.c,
diff/ChangeLog, diff/Makefile.am, diff/Makefile.in, diff/diff.c,
diff/gsl_diff.h, diff/test.c, doc/12-cities.eps, doc/ChangeLog,
doc/Makefile.am, doc/Makefile.in, doc/algorithm.sty,
doc/algorithmic.sty, doc/autoconf.texi, doc/blas.texi,
doc/bspline.eps, doc/bspline.texi, doc/calc.sty, doc/cblas.texi,
doc/cheb.eps, doc/cheb.texi, doc/combination.texi, doc/complex.texi,
doc/const.texi, doc/debug.texi, doc/dht.texi, doc/diff.texi,
doc/dwt-orig.eps, doc/dwt-samp.eps, doc/dwt.texi, doc/eigen.texi,
doc/err.texi, doc/fdl.texi, doc/fft-complex-radix2-f.eps,
doc/fft-complex-radix2-t.eps, doc/fft-complex-radix2.eps,
doc/fft-real-mixedradix.eps, doc/fft.texi, doc/fftalgorithms.bib,
doc/fftalgorithms.tex, doc/final-route.eps, doc/fit-exp.eps,
doc/fit-wlinear.eps, doc/fit-wlinear2.eps, doc/fitting.texi,
doc/freemanuals.texi, doc/gpl.texi, doc/gsl-config.1,
doc/gsl-design.texi, doc/gsl-histogram.1, doc/gsl-randist.1,
doc/gsl-ref.info, doc/gsl-ref.info-1, doc/gsl-ref.info-2,
doc/gsl-ref.info-3, doc/gsl-ref.info-4, doc/gsl-ref.info-5,
doc/gsl-ref.info-6, doc/gsl-ref.texi, doc/gsl.3, doc/histogram.eps,
doc/histogram.texi, doc/histogram2d.eps, doc/ieee754.texi,
doc/initial-route.eps, doc/integration.texi, doc/interp.texi,
doc/interp2.eps, doc/interpp2.eps, doc/intro.texi, doc/landau.dat,
doc/linalg.texi, doc/math.texi, doc/mdate-sh, doc/min-interval.eps,
doc/min.texi, doc/montecarlo.texi, doc/multifit.texi,
doc/multimin.eps, doc/multimin.texi, doc/multiroots.texi,
doc/ntuple.eps, doc/ntuple.texi, doc/ode-initval.texi,
doc/permutation.texi, doc/poly.texi, doc/qrng.eps, doc/qrng.texi,
doc/rand-bernoulli.tex, doc/rand-beta.tex, doc/rand-binomial.tex,
doc/rand-bivariate-gaussian.tex, doc/rand-cauchy.tex,
doc/rand-chisq.tex, doc/rand-erlang.tex, doc/rand-exponential.tex,
doc/rand-exppow.tex, doc/rand-fdist.tex, doc/rand-flat.tex,
doc/rand-gamma.tex, doc/rand-gaussian-tail.tex,
doc/rand-gaussian.tex, doc/rand-geometric.tex, doc/rand-gumbel.tex,
doc/rand-gumbel1.tex, doc/rand-gumbel2.tex,
doc/rand-hypergeometric.tex, doc/rand-landau.tex,
doc/rand-laplace.tex, doc/rand-levy.tex, doc/rand-levyskew.tex,
doc/rand-logarithmic.tex, doc/rand-logistic.tex,
doc/rand-lognormal.tex, doc/rand-nbinomial.tex, doc/rand-pareto.tex,
doc/rand-pascal.tex, doc/rand-poisson.tex,
doc/rand-rayleigh-tail.tex, doc/rand-rayleigh.tex,
doc/rand-tdist.tex, doc/rand-weibull.tex, doc/randist.texi,
doc/random-walk.tex, doc/randplots.gnp, doc/rng.texi,
doc/roots-bisection.eps, doc/roots-false-position.eps,
doc/roots-newtons-method.eps, doc/roots-secant-method.eps,
doc/roots.texi, doc/siman-energy.eps, doc/siman-test.eps,
doc/siman.texi, doc/sort.texi, doc/specfunc-airy.texi,
doc/specfunc-bessel.texi, doc/specfunc-clausen.texi,
doc/specfunc-coulomb.texi, doc/specfunc-coupling.texi,
doc/specfunc-dawson.texi, doc/specfunc-debye.texi,
doc/specfunc-dilog.texi, doc/specfunc-elementary.texi,
doc/specfunc-ellint.texi, doc/specfunc-elljac.texi,
doc/specfunc-erf.texi, doc/specfunc-exp.texi,
doc/specfunc-expint.texi, doc/specfunc-fermi-dirac.texi,
doc/specfunc-gamma.texi, doc/specfunc-gegenbauer.texi,
doc/specfunc-hyperg.texi, doc/specfunc-laguerre.texi,
doc/specfunc-lambert.texi, doc/specfunc-legendre.texi,
doc/specfunc-log.texi, doc/specfunc-mathieu.texi,
doc/specfunc-pow-int.texi, doc/specfunc-psi.texi,
doc/specfunc-synchrotron.texi, doc/specfunc-transport.texi,
doc/specfunc-trig.texi, doc/specfunc-zeta.texi, doc/specfunc.texi,
doc/stamp-vti, doc/statistics.texi, doc/sum.texi, doc/texinfo.tex,
doc/usage.texi, doc/vdp.eps, doc/vectors.texi, doc/version-ref.texi,
doc/examples/blas.c, doc/examples/blas.out, doc/examples/block.c,
doc/examples/block.out, doc/examples/bspline.c, doc/examples/cblas.c,
doc/examples/cblas.out, doc/examples/cdf.c, doc/examples/cdf.out,
doc/examples/cheb.c, doc/examples/combination.c,
doc/examples/combination.out, doc/examples/const.c,
doc/examples/const.out, doc/examples/demo_fn.c,
doc/examples/demo_fn.h, doc/examples/diff.c, doc/examples/diff.out,
doc/examples/dwt.c, doc/examples/dwt.dat, doc/examples/ecg.dat,
doc/examples/eigen.c, doc/examples/eigen_nonsymm.c,
doc/examples/expfit.c, doc/examples/fft.c, doc/examples/fftmr.c,
doc/examples/fftreal.c, doc/examples/fitting.c,
doc/examples/fitting2.c, doc/examples/fitting3.c,
doc/examples/histogram.c, doc/examples/histogram2d.c,
doc/examples/ieee.c, doc/examples/ieeeround.c,
doc/examples/integration.c, doc/examples/integration.out,
doc/examples/interp.c, doc/examples/interpp.c, doc/examples/intro.c,
doc/examples/intro.out, doc/examples/linalglu.c,
doc/examples/linalglu.out, doc/examples/matrix.c,
doc/examples/matrixw.c, doc/examples/min.c, doc/examples/min.out,
doc/examples/monte.c, doc/examples/nlfit.c, doc/examples/ntupler.c,
doc/examples/ntuplew.c, doc/examples/ode-initval.c,
doc/examples/odefixed.c, doc/examples/permseq.c,
doc/examples/permshuffle.c, doc/examples/polyroots.c,
doc/examples/polyroots.out, doc/examples/qrng.c,
doc/examples/randpoisson.2.out, doc/examples/randpoisson.c,
doc/examples/randpoisson.out, doc/examples/randwalk.c,
doc/examples/rng.c, doc/examples/rng.out, doc/examples/rngunif.2.out,
doc/examples/rngunif.c, doc/examples/rngunif.out,
doc/examples/rootnewt.c, doc/examples/roots.c, doc/examples/siman.c,
doc/examples/sortsmall.c, doc/examples/sortsmall.out,
doc/examples/specfun.c, doc/examples/specfun.out,
doc/examples/specfun_e.c, doc/examples/specfun_e.out,
doc/examples/stat.c, doc/examples/stat.out, doc/examples/statsort.c,
doc/examples/statsort.out, doc/examples/sum.c, doc/examples/sum.out,
doc/examples/vector.c, doc/examples/vectorr.c,
doc/examples/vectorview.c, doc/examples/vectorview.out,
doc/examples/vectorw.c, eigen/ChangeLog, eigen/Makefile.am,
eigen/Makefile.in, eigen/TODO, eigen/francis.c, eigen/gsl_eigen.h,
eigen/herm.c, eigen/hermv.c, eigen/jacobi.c, eigen/nonsymm.c,
eigen/nonsymmv.c, eigen/qrstep.c, eigen/schur.c, eigen/schur.h,
eigen/sort.c, eigen/symm.c, eigen/symmv.c, eigen/test.c,
err/ChangeLog, err/Makefile.am, err/Makefile.in, err/TODO,
err/error.c, err/gsl_errno.h, err/gsl_message.h, err/message.c,
err/stream.c, err/strerror.c, err/test.c, fft/ChangeLog,
fft/Makefile.am, fft/Makefile.in, fft/TODO, fft/bitreverse.c,
fft/bitreverse.h, fft/c_init.c, fft/c_main.c, fft/c_pass.h,
fft/c_pass_2.c, fft/c_pass_3.c, fft/c_pass_4.c, fft/c_pass_5.c,
fft/c_pass_6.c, fft/c_pass_7.c, fft/c_pass_n.c, fft/c_radix2.c,
fft/compare.h, fft/compare_source.c, fft/complex_internal.h,
fft/dft.c, fft/dft_source.c, fft/factorize.c, fft/factorize.h,
fft/fft.c, fft/gsl_dft_complex.h, fft/gsl_dft_complex_float.h,
fft/gsl_fft.h, fft/gsl_fft_complex.h, fft/gsl_fft_complex_float.h,
fft/gsl_fft_halfcomplex.h, fft/gsl_fft_halfcomplex_float.h,
fft/gsl_fft_real.h, fft/gsl_fft_real_float.h, fft/hc_init.c,
fft/hc_main.c, fft/hc_pass.h, fft/hc_pass_2.c, fft/hc_pass_3.c,
fft/hc_pass_4.c, fft/hc_pass_5.c, fft/hc_pass_n.c, fft/hc_radix2.c,
fft/hc_unpack.c, fft/real_init.c, fft/real_main.c, fft/real_pass.h,
fft/real_pass_2.c, fft/real_pass_3.c, fft/real_pass_4.c,
fft/real_pass_5.c, fft/real_pass_n.c, fft/real_radix2.c,
fft/real_unpack.c, fft/signals.c, fft/signals.h,
fft/signals_source.c, fft/test.c, fft/test_complex_source.c,
fft/test_real_source.c, fft/test_trap_source.c, fft/urand.c,
fit/ChangeLog, fit/Makefile.am, fit/Makefile.in, fit/gsl_fit.h,
fit/linear.c, fit/test.c, gsl/Makefile.am, gsl/Makefile.in,
histogram/ChangeLog, histogram/Makefile.am, histogram/Makefile.in,
histogram/TODO, histogram/add.c, histogram/add2d.c,
histogram/calloc_range.c, histogram/calloc_range2d.c,
histogram/copy.c, histogram/copy2d.c, histogram/file.c,
histogram/file2d.c, histogram/find.c, histogram/find2d.c,
histogram/get.c, histogram/get2d.c, histogram/gsl_histogram.h,
histogram/gsl_histogram2d.h, histogram/init.c, histogram/init2d.c,
histogram/maxval.c, histogram/maxval2d.c, histogram/oper.c,
histogram/oper2d.c, histogram/params.c, histogram/params2d.c,
histogram/pdf.c, histogram/pdf2d.c, histogram/reset.c,
histogram/reset2d.c, histogram/stat.c, histogram/stat2d.c,
histogram/test.c, histogram/test1d.c, histogram/test1d_resample.c,
histogram/test1d_trap.c, histogram/test2d.c,
histogram/test2d_resample.c, histogram/test2d_trap.c,
histogram/urand.c, ieee-utils/ChangeLog, ieee-utils/Makefile.am,
ieee-utils/Makefile.in, ieee-utils/TODO, ieee-utils/endian.c,
ieee-utils/env.c, ieee-utils/fp-aix.c, ieee-utils/fp-darwin.c,
ieee-utils/fp-darwin86.c, ieee-utils/fp-freebsd.c,
ieee-utils/fp-gnuc99.c, ieee-utils/fp-gnum68k.c,
ieee-utils/fp-gnuppc.c, ieee-utils/fp-gnusparc.c,
ieee-utils/fp-gnux86.c, ieee-utils/fp-hpux.c, ieee-utils/fp-hpux11.c,
ieee-utils/fp-irix.c, ieee-utils/fp-netbsd.c,
ieee-utils/fp-openbsd.c, ieee-utils/fp-os2emx.c,
ieee-utils/fp-solaris.c, ieee-utils/fp-sunos4.c,
ieee-utils/fp-tru64.c, ieee-utils/fp-unknown.c, ieee-utils/fp.c,
ieee-utils/gsl_ieee_utils.h, ieee-utils/make_rep.c,
ieee-utils/print.c, ieee-utils/read.c, ieee-utils/standardize.c,
ieee-utils/test.c, integration/ChangeLog, integration/Makefile.am,
integration/Makefile.in, integration/TODO, integration/append.c,
integration/err.c, integration/gsl_integration.h,
integration/initialise.c, integration/positivity.c,
integration/ptsort.c, integration/qag.c, integration/qagp.c,
integration/qags.c, integration/qawc.c, integration/qawf.c,
integration/qawo.c, integration/qaws.c, integration/qc25c.c,
integration/qc25f.c, integration/qc25s.c, integration/qcheb.c,
integration/qelg.c, integration/qk.c, integration/qk15.c,
integration/qk21.c, integration/qk31.c, integration/qk41.c,
integration/qk51.c, integration/qk61.c, integration/qmomo.c,
integration/qmomof.c, integration/qng.c, integration/qng.h,
integration/qpsrt.c, integration/qpsrt2.c, integration/reset.c,
integration/set_initial.c, integration/test.c, integration/tests.c,
integration/tests.h, integration/util.c, integration/workspace.c,
interpolation/ChangeLog, interpolation/Makefile.am,
interpolation/Makefile.in, interpolation/TODO, interpolation/accel.c,
interpolation/akima.c, interpolation/bsearch.c,
interpolation/bsearch.h, interpolation/cspline.c,
interpolation/gsl_interp.h, interpolation/gsl_spline.h,
interpolation/integ_eval.h, interpolation/interp.c,
interpolation/linear.c, interpolation/poly.c, interpolation/spline.c,
interpolation/test.c, linalg/ChangeLog, linalg/Makefile.am,
linalg/Makefile.in, linalg/TODO, linalg/apply_givens.c,
linalg/balance.c, linalg/balancemat.c, linalg/bidiag.c,
linalg/cholesky.c, linalg/exponential.c, linalg/givens.c,
linalg/gsl_linalg.h, linalg/hermtd.c, linalg/hessenberg.c,
linalg/hh.c, linalg/householder.c, linalg/householdercomplex.c,
linalg/lq.c, linalg/lu.c, linalg/luc.c, linalg/multiply.c,
linalg/ptlq.c, linalg/qr.c, linalg/qrpt.c, linalg/svd.c,
linalg/svdstep.c, linalg/symmtd.c, linalg/test.c, linalg/tridiag.c,
linalg/tridiag.h, matrix/ChangeLog, matrix/Makefile.am,
matrix/Makefile.in, matrix/TODO, matrix/copy.c, matrix/copy_source.c,
matrix/file.c, matrix/file_source.c, matrix/getset.c,
matrix/getset_source.c, matrix/gsl_matrix.h,
matrix/gsl_matrix_char.h, matrix/gsl_matrix_complex_double.h,
matrix/gsl_matrix_complex_float.h,
matrix/gsl_matrix_complex_long_double.h, matrix/gsl_matrix_double.h,
matrix/gsl_matrix_float.h, matrix/gsl_matrix_int.h,
matrix/gsl_matrix_long.h, matrix/gsl_matrix_long_double.h,
matrix/gsl_matrix_short.h, matrix/gsl_matrix_uchar.h,
matrix/gsl_matrix_uint.h, matrix/gsl_matrix_ulong.h,
matrix/gsl_matrix_ushort.h, matrix/init.c, matrix/init_source.c,
matrix/matrix.c, matrix/matrix_source.c, matrix/minmax.c,
matrix/minmax_source.c, matrix/oper.c, matrix/oper_complex_source.c,
matrix/oper_source.c, matrix/prop.c, matrix/prop_source.c,
matrix/rowcol.c, matrix/rowcol_source.c, matrix/submatrix.c,
matrix/submatrix_source.c, matrix/swap.c, matrix/swap_source.c,
matrix/test.c, matrix/test_complex_source.c, matrix/test_source.c,
matrix/test_static.c, matrix/view.c, matrix/view.h,
matrix/view_source.c, min/ChangeLog, min/Makefile.am,
min/Makefile.in, min/bracketing.c, min/brent.c, min/convergence.c,
min/fsolver.c, min/golden.c, min/gsl_min.h, min/min.h, min/test.c,
min/test.h, min/test_funcs.c, monte/ChangeLog, monte/Makefile.am,
monte/Makefile.in, monte/README, monte/TODO, monte/gsl_monte.h,
monte/gsl_monte_miser.h, monte/gsl_monte_plain.h,
monte/gsl_monte_vegas.h, monte/miser.c, monte/plain.c, monte/test.c,
monte/test_main.c, monte/vegas.c, multifit/ChangeLog,
multifit/Makefile.am, multifit/Makefile.in, multifit/TODO,
multifit/convergence.c, multifit/covar.c, multifit/fdfsolver.c,
multifit/fsolver.c, multifit/gradient.c, multifit/gsl_multifit.h,
multifit/gsl_multifit_nlin.h, multifit/lmder.c, multifit/lmiterate.c,
multifit/lmpar.c, multifit/lmset.c, multifit/lmutil.c,
multifit/multilinear.c, multifit/qrsolv.c, multifit/test.c,
multifit/test_brown.c, multifit/test_enso.c,
multifit/test_estimator.c, multifit/test_filip.c, multifit/test_fn.c,
multifit/test_hahn1.c, multifit/test_kirby2.c,
multifit/test_longley.c, multifit/test_nelson.c,
multifit/test_pontius.c, multifit/work.c, multimin/ChangeLog,
multimin/Makefile.am, multimin/Makefile.in, multimin/TODO,
multimin/conjugate_fr.c, multimin/conjugate_pr.c,
multimin/convergence.c, multimin/diff.c,
multimin/directional_minimize.c, multimin/fdfminimizer.c,
multimin/fminimizer.c, multimin/gsl_multimin.h,
multimin/linear_minimize.c, multimin/linear_wrapper.c,
multimin/simplex.c, multimin/steepest_descent.c, multimin/test.c,
multimin/test_funcs.c, multimin/test_funcs.h, multimin/vector_bfgs.c,
multimin/vector_bfgs2.c, multiroots/ChangeLog,
multiroots/Makefile.am, multiroots/Makefile.in, multiroots/broyden.c,
multiroots/convergence.c, multiroots/dnewton.c, multiroots/dogleg.c,
multiroots/enorm.c, multiroots/fdfsolver.c, multiroots/fdjac.c,
multiroots/fsolver.c, multiroots/gnewton.c,
multiroots/gsl_multiroots.h, multiroots/hybrid.c,
multiroots/hybridj.c, multiroots/newton.c, multiroots/test.c,
multiroots/test_funcs.c, multiroots/test_funcs.h, ntuple/ChangeLog,
ntuple/Makefile.am, ntuple/Makefile.in, ntuple/gsl_ntuple.h,
ntuple/ntuple.c, ntuple/test.c, ode-initval/ChangeLog,
ode-initval/Makefile.am, ode-initval/Makefile.in, ode-initval/TODO,
ode-initval/bsimp.c, ode-initval/control.c, ode-initval/cscal.c,
ode-initval/cstd.c, ode-initval/evolve.c, ode-initval/gear1.c,
ode-initval/gear2.c, ode-initval/gsl_odeiv.h,
ode-initval/odeiv_util.h, ode-initval/rk2.c, ode-initval/rk2imp.c,
ode-initval/rk2simp.c, ode-initval/rk4.c, ode-initval/rk4imp.c,
ode-initval/rk8pd.c, ode-initval/rkck.c, ode-initval/rkf45.c,
ode-initval/step.c, ode-initval/test.c, permutation/ChangeLog,
permutation/Makefile.am, permutation/Makefile.in,
permutation/canonical.c, permutation/file.c,
permutation/gsl_permutation.h, permutation/gsl_permute.h,
permutation/gsl_permute_char.h,
permutation/gsl_permute_complex_double.h,
permutation/gsl_permute_complex_float.h,
permutation/gsl_permute_complex_long_double.h,
permutation/gsl_permute_double.h, permutation/gsl_permute_float.h,
permutation/gsl_permute_int.h, permutation/gsl_permute_long.h,
permutation/gsl_permute_long_double.h,
permutation/gsl_permute_short.h, permutation/gsl_permute_uchar.h,
permutation/gsl_permute_uint.h, permutation/gsl_permute_ulong.h,
permutation/gsl_permute_ushort.h, permutation/gsl_permute_vector.h,
permutation/gsl_permute_vector_char.h,
permutation/gsl_permute_vector_complex_double.h,
permutation/gsl_permute_vector_complex_float.h,
permutation/gsl_permute_vector_complex_long_double.h,
permutation/gsl_permute_vector_double.h,
permutation/gsl_permute_vector_float.h,
permutation/gsl_permute_vector_int.h,
permutation/gsl_permute_vector_long.h,
permutation/gsl_permute_vector_long_double.h,
permutation/gsl_permute_vector_short.h,
permutation/gsl_permute_vector_uchar.h,
permutation/gsl_permute_vector_uint.h,
permutation/gsl_permute_vector_ulong.h,
permutation/gsl_permute_vector_ushort.h, permutation/init.c,
permutation/permutation.c, permutation/permute.c,
permutation/permute_source.c, permutation/test.c, poly/ChangeLog,
poly/Makefile.am, poly/Makefile.in, poly/TODO, poly/balance.c,
poly/companion.c, poly/dd.c, poly/eval.c, poly/gsl_poly.h, poly/qr.c,
poly/solve_cubic.c, poly/solve_quadratic.c, poly/test.c,
poly/zsolve.c, poly/zsolve_cubic.c, poly/zsolve_init.c,
poly/zsolve_quadratic.c, qrng/ChangeLog, qrng/Makefile.am,
qrng/Makefile.in, qrng/TODO, qrng/gsl_qrng.h, qrng/niederreiter-2.c,
qrng/qrng.c, qrng/sobol.c, qrng/test.c, randist/ChangeLog,
randist/Makefile.am, randist/Makefile.in, randist/TODO,
randist/bernoulli.c, randist/beta.c, randist/bigauss.c,
randist/binomial.c, randist/binomial_tpe.c, randist/cauchy.c,
randist/chisq.c, randist/dirichlet.c, randist/discrete.c,
randist/erlang.c, randist/exponential.c, randist/exppow.c,
randist/fdist.c, randist/flat.c, randist/gamma.c, randist/gauss.c,
randist/gausstail.c, randist/gausszig.c, randist/geometric.c,
randist/gsl_randist.h, randist/gumbel.c, randist/hyperg.c,
randist/landau.c, randist/laplace.c, randist/levy.c,
randist/logarithmic.c, randist/logistic.c, randist/lognormal.c,
randist/multinomial.c, randist/nbinomial.c, randist/pareto.c,
randist/pascal.c, randist/poisson.c, randist/rayleigh.c,
randist/shuffle.c, randist/sphere.c, randist/tdist.c, randist/test.c,
randist/weibull.c, rng/ChangeLog, rng/Makefile.am, rng/Makefile.in,
rng/TODO, rng/borosh13.c, rng/cmrg.c, rng/coveyou.c, rng/default.c,
rng/file.c, rng/fishman18.c, rng/fishman20.c, rng/fishman2x.c,
rng/gfsr4.c, rng/gsl_rng.h, rng/knuthran.c, rng/knuthran2.c,
rng/knuthran2002.c, rng/lecuyer21.c, rng/minstd.c, rng/mrg.c,
rng/mt.c, rng/r250.c, rng/ran0.c, rng/ran1.c, rng/ran2.c, rng/ran3.c,
rng/rand.c, rng/rand48.c, rng/random.c, rng/randu.c, rng/ranf.c,
rng/ranlux.c, rng/ranlxd.c, rng/ranlxs.c, rng/ranmar.c, rng/rng.c,
rng/schrage.c, rng/slatec.c, rng/taus.c, rng/taus113.c, rng/test.c,
rng/transputer.c, rng/tt.c, rng/types.c, rng/uni.c, rng/uni32.c,
rng/vax.c, rng/waterman14.c, rng/zuf.c, roots/ChangeLog,
roots/Makefile.am, roots/Makefile.in, roots/TODO, roots/bisection.c,
roots/brent.c, roots/convergence.c, roots/falsepos.c,
roots/fdfsolver.c, roots/fsolver.c, roots/gsl_roots.h,
roots/newton.c, roots/roots.h, roots/secant.c, roots/steffenson.c,
roots/test.c, roots/test.h, roots/test_funcs.c, siman/ChangeLog,
siman/Makefile.am, siman/Makefile.in, siman/TODO, siman/gsl_siman.h,
siman/siman.c, siman/siman_test_driver.sh, siman/siman_tsp.c,
siman/test.c, sort/ChangeLog, sort/Makefile.am, sort/Makefile.in,
sort/TODO, sort/gsl_heapsort.h, sort/gsl_sort.h,
sort/gsl_sort_char.h, sort/gsl_sort_double.h, sort/gsl_sort_float.h,
sort/gsl_sort_int.h, sort/gsl_sort_long.h,
sort/gsl_sort_long_double.h, sort/gsl_sort_short.h,
sort/gsl_sort_uchar.h, sort/gsl_sort_uint.h, sort/gsl_sort_ulong.h,
sort/gsl_sort_ushort.h, sort/gsl_sort_vector.h,
sort/gsl_sort_vector_char.h, sort/gsl_sort_vector_double.h,
sort/gsl_sort_vector_float.h, sort/gsl_sort_vector_int.h,
sort/gsl_sort_vector_long.h, sort/gsl_sort_vector_long_double.h,
sort/gsl_sort_vector_short.h, sort/gsl_sort_vector_uchar.h,
sort/gsl_sort_vector_uint.h, sort/gsl_sort_vector_ulong.h,
sort/gsl_sort_vector_ushort.h, sort/sort.c, sort/sortind.c,
sort/sortvec.c, sort/sortvec_source.c, sort/sortvecind.c,
sort/sortvecind_source.c, sort/subset.c, sort/subset_source.c,
sort/subsetind.c, sort/subsetind_source.c, sort/test.c,
sort/test_heapsort.c, sort/test_source.c, specfunc/ChangeLog,
specfunc/Makefile.am, specfunc/Makefile.in, specfunc/TODO,
specfunc/airy.c, specfunc/airy_der.c, specfunc/airy_zero.c,
specfunc/atanint.c, specfunc/bessel.c, specfunc/bessel.h,
specfunc/bessel_I0.c, specfunc/bessel_I1.c, specfunc/bessel_In.c,
specfunc/bessel_Inu.c, specfunc/bessel_J0.c, specfunc/bessel_J1.c,
specfunc/bessel_Jn.c, specfunc/bessel_Jnu.c, specfunc/bessel_K0.c,
specfunc/bessel_K1.c, specfunc/bessel_Kn.c, specfunc/bessel_Knu.c,
specfunc/bessel_Y0.c, specfunc/bessel_Y1.c, specfunc/bessel_Yn.c,
specfunc/bessel_Ynu.c, specfunc/bessel_amp_phase.c,
specfunc/bessel_amp_phase.h, specfunc/bessel_i.c,
specfunc/bessel_j.c, specfunc/bessel_k.c, specfunc/bessel_olver.c,
specfunc/bessel_olver.h, specfunc/bessel_sequence.c,
specfunc/bessel_temme.c, specfunc/bessel_temme.h,
specfunc/bessel_y.c, specfunc/bessel_zero.c, specfunc/beta.c,
specfunc/beta_inc.c, specfunc/cheb_eval.c, specfunc/cheb_eval_mode.c,
specfunc/chebyshev.h, specfunc/check.h, specfunc/clausen.c,
specfunc/coulomb.c, specfunc/coulomb_bound.c, specfunc/coupling.c,
specfunc/dawson.c, specfunc/debye.c, specfunc/dilog.c,
specfunc/elementary.c, specfunc/ellint.c, specfunc/elljac.c,
specfunc/erfc.c, specfunc/error.h, specfunc/eval.h, specfunc/exp.c,
specfunc/expint.c, specfunc/expint3.c, specfunc/fermi_dirac.c,
specfunc/gamma.c, specfunc/gamma_inc.c, specfunc/gegenbauer.c,
specfunc/gsl_sf.h, specfunc/gsl_sf_airy.h, specfunc/gsl_sf_bessel.h,
specfunc/gsl_sf_clausen.h, specfunc/gsl_sf_coulomb.h,
specfunc/gsl_sf_coupling.h, specfunc/gsl_sf_dawson.h,
specfunc/gsl_sf_debye.h, specfunc/gsl_sf_dilog.h,
specfunc/gsl_sf_elementary.h, specfunc/gsl_sf_ellint.h,
specfunc/gsl_sf_elljac.h, specfunc/gsl_sf_erf.h,
specfunc/gsl_sf_exp.h, specfunc/gsl_sf_expint.h,
specfunc/gsl_sf_fermi_dirac.h, specfunc/gsl_sf_gamma.h,
specfunc/gsl_sf_gegenbauer.h, specfunc/gsl_sf_hyperg.h,
specfunc/gsl_sf_laguerre.h, specfunc/gsl_sf_lambert.h,
specfunc/gsl_sf_legendre.h, specfunc/gsl_sf_log.h,
specfunc/gsl_sf_mathieu.h, specfunc/gsl_sf_pow_int.h,
specfunc/gsl_sf_psi.h, specfunc/gsl_sf_result.h,
specfunc/gsl_sf_synchrotron.h, specfunc/gsl_sf_transport.h,
specfunc/gsl_sf_trig.h, specfunc/gsl_sf_zeta.h,
specfunc/gsl_specfunc.h, specfunc/hyperg.c, specfunc/hyperg.h,
specfunc/hyperg_0F1.c, specfunc/hyperg_1F1.c, specfunc/hyperg_2F0.c,
specfunc/hyperg_2F1.c, specfunc/hyperg_U.c, specfunc/laguerre.c,
specfunc/lambert.c, specfunc/legendre.h, specfunc/legendre_H3d.c,
specfunc/legendre_Qn.c, specfunc/legendre_con.c,
specfunc/legendre_poly.c, specfunc/log.c, specfunc/mathieu_angfunc.c,
specfunc/mathieu_charv.c, specfunc/mathieu_coeff.c,
specfunc/mathieu_radfunc.c, specfunc/mathieu_workspace.c,
specfunc/poch.c, specfunc/pow_int.c, specfunc/psi.c,
specfunc/recurse.h, specfunc/result.c, specfunc/shint.c,
specfunc/sinint.c, specfunc/synchrotron.c, specfunc/test_airy.c,
specfunc/test_bessel.c, specfunc/test_coulomb.c,
specfunc/test_dilog.c, specfunc/test_gamma.c, specfunc/test_hyperg.c,
specfunc/test_legendre.c, specfunc/test_mathieu.c,
specfunc/test_sf.c, specfunc/test_sf.h, specfunc/transport.c,
specfunc/trig.c, specfunc/zeta.c, statistics/ChangeLog,
statistics/Makefile.am, statistics/Makefile.in, statistics/TODO,
statistics/absdev.c, statistics/absdev_source.c,
statistics/covariance.c, statistics/covariance_source.c,
statistics/gsl_statistics.h, statistics/gsl_statistics_char.h,
statistics/gsl_statistics_double.h,
statistics/gsl_statistics_float.h, statistics/gsl_statistics_int.h,
statistics/gsl_statistics_long.h,
statistics/gsl_statistics_long_double.h,
statistics/gsl_statistics_short.h, statistics/gsl_statistics_uchar.h,
statistics/gsl_statistics_uint.h, statistics/gsl_statistics_ulong.h,
statistics/gsl_statistics_ushort.h, statistics/kurtosis.c,
statistics/kurtosis_source.c, statistics/lag1.c,
statistics/lag1_source.c, statistics/mean.c,
statistics/mean_source.c, statistics/median.c,
statistics/median_source.c, statistics/minmax.c,
statistics/minmax_source.c, statistics/p_variance.c,
statistics/p_variance_source.c, statistics/quantiles.c,
statistics/quantiles_source.c, statistics/skew.c,
statistics/skew_source.c, statistics/test.c,
statistics/test_float_source.c, statistics/test_int_source.c,
statistics/test_nist.c, statistics/ttest.c,
statistics/ttest_source.c, statistics/variance.c,
statistics/variance_source.c, statistics/wabsdev.c,
statistics/wabsdev_source.c, statistics/wkurtosis.c,
statistics/wkurtosis_source.c, statistics/wmean.c,
statistics/wmean_source.c, statistics/wskew.c,
statistics/wskew_source.c, statistics/wvariance.c,
statistics/wvariance_source.c, sum/ChangeLog, sum/Makefile.am,
sum/Makefile.in, sum/gsl_sum.h, sum/levin_u.c, sum/levin_utrunc.c,
sum/test.c, sum/work_u.c, sum/work_utrunc.c, sys/ChangeLog,
sys/Makefile.am, sys/Makefile.in, sys/coerce.c, sys/expm1.c,
sys/fcmp.c, sys/fdiv.c, sys/gsl_sys.h, sys/hypot.c, sys/infnan.c,
sys/invhyp.c, sys/ldfrexp.c, sys/log1p.c, sys/minmax.c,
sys/pow_int.c, sys/prec.c, sys/test.c, test/ChangeLog,
test/Makefile.am, test/Makefile.in, test/gsl_test.h, test/results.c,
utils/Makefile.am, utils/Makefile.in, utils/README, utils/memcpy.c,
utils/memmove.c, utils/placeholder.c, utils/strdup.c, utils/strtol.c,
utils/strtoul.c, utils/system.h, vector/ChangeLog,
vector/Makefile.am, vector/Makefile.in, vector/TODO, vector/copy.c,
vector/copy_source.c, vector/file.c, vector/file_source.c,
vector/gsl_vector.h, vector/gsl_vector_char.h,
vector/gsl_vector_complex.h, vector/gsl_vector_complex_double.h,
vector/gsl_vector_complex_float.h,
vector/gsl_vector_complex_long_double.h, vector/gsl_vector_double.h,
vector/gsl_vector_float.h, vector/gsl_vector_int.h,
vector/gsl_vector_long.h, vector/gsl_vector_long_double.h,
vector/gsl_vector_short.h, vector/gsl_vector_uchar.h,
vector/gsl_vector_uint.h, vector/gsl_vector_ulong.h,
vector/gsl_vector_ushort.h, vector/init.c, vector/init_source.c,
vector/minmax.c, vector/minmax_source.c, vector/oper.c,
vector/oper_source.c, vector/prop.c, vector/prop_source.c,
vector/reim.c, vector/reim_source.c, vector/subvector.c,
vector/subvector_source.c, vector/swap.c, vector/swap_source.c,
vector/test.c, vector/test_complex_source.c, vector/test_source.c,
vector/test_static.c, vector/vector.c, vector/vector_source.c,
vector/view.c, vector/view.h, vector/view_source.c,
wavelet/ChangeLog, wavelet/Makefile.am, wavelet/Makefile.in,
wavelet/TODO, wavelet/bspline.c, wavelet/daubechies.c, wavelet/dwt.c,
wavelet/gsl_wavelet.h, wavelet/gsl_wavelet2d.h, wavelet/haar.c,
wavelet/test.c, wavelet/wavelet.c: New files.
Diffstat (limited to 'gsl-1.9/specfunc/gamma.c')
| -rw-r--r-- | gsl-1.9/specfunc/gamma.c | 1685 |
1 files changed, 1685 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/gamma.c b/gsl-1.9/specfunc/gamma.c new file mode 100644 index 0000000..4703cee --- /dev/null +++ b/gsl-1.9/specfunc/gamma.c @@ -0,0 +1,1685 @@ +/* specfunc/gamma.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_exp.h> +#include <gsl/gsl_sf_log.h> +#include <gsl/gsl_sf_psi.h> +#include <gsl/gsl_sf_trig.h> +#include <gsl/gsl_sf_gamma.h> + +#include "error.h" +#include "check.h" + +#include "chebyshev.h" +#include "cheb_eval.c" + +#define LogRootTwoPi_ 0.9189385332046727418 + + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +static struct {int n; double f; long i; } fact_table[GSL_SF_FACT_NMAX + 1] = { + { 0, 1.0, 1L }, + { 1, 1.0, 1L }, + { 2, 2.0, 2L }, + { 3, 6.0, 6L }, + { 4, 24.0, 24L }, + { 5, 120.0, 120L }, + { 6, 720.0, 720L }, + { 7, 5040.0, 5040L }, + { 8, 40320.0, 40320L }, + + { 9, 362880.0, 362880L }, + { 10, 3628800.0, 3628800L }, + { 11, 39916800.0, 39916800L }, + { 12, 479001600.0, 479001600L }, + + { 13, 6227020800.0, 0 }, + { 14, 87178291200.0, 0 }, + { 15, 1307674368000.0, 0 }, + { 16, 20922789888000.0, 0 }, + { 17, 355687428096000.0, 0 }, + { 18, 6402373705728000.0, 0 }, + { 19, 121645100408832000.0, 0 }, + { 20, 2432902008176640000.0, 0 }, + { 21, 51090942171709440000.0, 0 }, + { 22, 1124000727777607680000.0, 0 }, + { 23, 25852016738884976640000.0, 0 }, + { 24, 620448401733239439360000.0, 0 }, + { 25, 15511210043330985984000000.0, 0 }, + { 26, 403291461126605635584000000.0, 0 }, + { 27, 10888869450418352160768000000.0, 0 }, + { 28, 304888344611713860501504000000.0, 0 }, + { 29, 8841761993739701954543616000000.0, 0 }, + { 30, 265252859812191058636308480000000.0, 0 }, + { 31, 8222838654177922817725562880000000.0, 0 }, + { 32, 263130836933693530167218012160000000.0, 0 }, + { 33, 8683317618811886495518194401280000000.0, 0 }, + { 34, 2.95232799039604140847618609644e38, 0 }, + { 35, 1.03331479663861449296666513375e40, 0 }, + { 36, 3.71993326789901217467999448151e41, 0 }, + { 37, 1.37637530912263450463159795816e43, 0 }, + { 38, 5.23022617466601111760007224100e44, 0 }, + { 39, 2.03978820811974433586402817399e46, 0 }, + { 40, 8.15915283247897734345611269600e47, 0 }, + { 41, 3.34525266131638071081700620534e49, 0 }, + { 42, 1.40500611775287989854314260624e51, 0 }, + { 43, 6.04152630633738356373551320685e52, 0 }, + { 44, 2.65827157478844876804362581101e54, 0 }, + { 45, 1.19622220865480194561963161496e56, 0 }, + { 46, 5.50262215981208894985030542880e57, 0 }, + { 47, 2.58623241511168180642964355154e59, 0 }, + { 48, 1.24139155925360726708622890474e61, 0 }, + { 49, 6.08281864034267560872252163321e62, 0 }, + { 50, 3.04140932017133780436126081661e64, 0 }, + { 51, 1.55111875328738228022424301647e66, 0 }, + { 52, 8.06581751709438785716606368564e67, 0 }, + { 53, 4.27488328406002556429801375339e69, 0 }, + { 54, 2.30843697339241380472092742683e71, 0 }, + { 55, 1.26964033536582759259651008476e73, 0 }, + { 56, 7.10998587804863451854045647464e74, 0 }, + { 57, 4.05269195048772167556806019054e76, 0 }, + { 58, 2.35056133128287857182947491052e78, 0 }, + { 59, 1.38683118545689835737939019720e80, 0 }, + { 60, 8.32098711274139014427634118320e81, 0 }, + { 61, 5.07580213877224798800856812177e83, 0 }, + { 62, 3.14699732603879375256531223550e85, 0 }, + { 63, 1.982608315404440064116146708360e87, 0 }, + { 64, 1.268869321858841641034333893350e89, 0 }, + { 65, 8.247650592082470666723170306800e90, 0 }, + { 66, 5.443449390774430640037292402480e92, 0 }, + { 67, 3.647111091818868528824985909660e94, 0 }, + { 68, 2.480035542436830599600990418570e96, 0 }, + { 69, 1.711224524281413113724683388810e98, 0 }, + { 70, 1.197857166996989179607278372170e100, 0 }, + { 71, 8.504785885678623175211676442400e101, 0 }, + { 72, 6.123445837688608686152407038530e103, 0 }, + { 73, 4.470115461512684340891257138130e105, 0 }, + { 74, 3.307885441519386412259530282210e107, 0 }, + { 75, 2.480914081139539809194647711660e109, 0 }, + { 76, 1.885494701666050254987932260860e111, 0 }, + { 77, 1.451830920282858696340707840860e113, 0 }, + { 78, 1.132428117820629783145752115870e115, 0 }, + { 79, 8.946182130782975286851441715400e116, 0 }, + { 80, 7.156945704626380229481153372320e118, 0 }, + { 81, 5.797126020747367985879734231580e120, 0 }, + { 82, 4.753643337012841748421382069890e122, 0 }, + { 83, 3.945523969720658651189747118010e124, 0 }, + { 84, 3.314240134565353266999387579130e126, 0 }, + { 85, 2.817104114380550276949479442260e128, 0 }, + { 86, 2.422709538367273238176552320340e130, 0 }, + { 87, 2.107757298379527717213600518700e132, 0 }, + { 88, 1.854826422573984391147968456460e134, 0 }, + { 89, 1.650795516090846108121691926250e136, 0 }, + { 90, 1.485715964481761497309522733620e138, 0 }, + { 91, 1.352001527678402962551665687590e140, 0 }, + { 92, 1.243841405464130725547532432590e142, 0 }, + { 93, 1.156772507081641574759205162310e144, 0 }, + { 94, 1.087366156656743080273652852570e146, 0 }, + { 95, 1.032997848823905926259970209940e148, 0 }, + { 96, 9.916779348709496892095714015400e149, 0 }, + { 97, 9.619275968248211985332842594960e151, 0 }, + { 98, 9.426890448883247745626185743100e153, 0 }, + { 99, 9.332621544394415268169923885600e155, 0 }, + { 100, 9.33262154439441526816992388563e157, 0 }, + { 101, 9.42594775983835942085162312450e159, 0 }, + { 102, 9.61446671503512660926865558700e161, 0 }, + { 103, 9.90290071648618040754671525458e163, 0 }, + { 104, 1.02990167451456276238485838648e166, 0 }, + { 105, 1.08139675824029090050410130580e168, 0 }, + { 106, 1.146280563734708354534347384148e170, 0 }, + { 107, 1.226520203196137939351751701040e172, 0 }, + { 108, 1.324641819451828974499891837120e174, 0 }, + { 109, 1.443859583202493582204882102460e176, 0 }, + { 110, 1.588245541522742940425370312710e178, 0 }, + { 111, 1.762952551090244663872161047110e180, 0 }, + { 112, 1.974506857221074023536820372760e182, 0 }, + { 113, 2.231192748659813646596607021220e184, 0 }, + { 114, 2.543559733472187557120132004190e186, 0 }, + { 115, 2.925093693493015690688151804820e188, 0 }, + { 116, 3.393108684451898201198256093590e190, 0 }, + { 117, 3.96993716080872089540195962950e192, 0 }, + { 118, 4.68452584975429065657431236281e194, 0 }, + { 119, 5.57458576120760588132343171174e196, 0 }, + { 120, 6.68950291344912705758811805409e198, 0 }, + { 121, 8.09429852527344373968162284545e200, 0 }, + { 122, 9.87504420083360136241157987140e202, 0 }, + { 123, 1.21463043670253296757662432419e205, 0 }, + { 124, 1.50614174151114087979501416199e207, 0 }, + { 125, 1.88267717688892609974376770249e209, 0 }, + { 126, 2.37217324288004688567714730514e211, 0 }, + { 127, 3.01266001845765954480997707753e213, 0 }, + { 128, 3.85620482362580421735677065923e215, 0 }, + { 129, 4.97450422247728744039023415041e217, 0 }, + { 130, 6.46685548922047367250730439554e219, 0 }, + { 131, 8.47158069087882051098456875820e221, 0 }, + { 132, 1.11824865119600430744996307608e224, 0 }, + { 133, 1.48727070609068572890845089118e226, 0 }, + { 134, 1.99294274616151887673732419418e228, 0 }, + { 135, 2.69047270731805048359538766215e230, 0 }, + { 136, 3.65904288195254865768972722052e232, 0 }, + { 137, 5.01288874827499166103492629211e234, 0 }, + { 138, 6.91778647261948849222819828311e236, 0 }, + { 139, 9.61572319694108900419719561353e238, 0 }, + { 140, 1.34620124757175246058760738589e241, 0 }, + { 141, 1.89814375907617096942852641411e243, 0 }, + { 142, 2.69536413788816277658850750804e245, 0 }, + { 143, 3.85437071718007277052156573649e247, 0 }, + { 144, 5.55029383273930478955105466055e249, 0 }, + { 145, 8.04792605747199194484902925780e251, 0 }, + { 146, 1.17499720439091082394795827164e254, 0 }, + { 147, 1.72724589045463891120349865931e256, 0 }, + { 148, 2.55632391787286558858117801578e258, 0 }, + { 149, 3.80892263763056972698595524351e260, 0 }, + { 150, 5.71338395644585459047893286526e262, 0 }, + { 151, 8.62720977423324043162318862650e264, 0 }, + { 152, 1.31133588568345254560672467123e267, 0 }, + { 153, 2.00634390509568239477828874699e269, 0 }, + { 154, 3.08976961384735088795856467036e271, 0 }, + { 155, 4.78914290146339387633577523906e273, 0 }, + { 156, 7.47106292628289444708380937294e275, 0 }, + { 157, 1.17295687942641442819215807155e278, 0 }, + { 158, 1.85327186949373479654360975305e280, 0 }, + { 159, 2.94670227249503832650433950735e282, 0 }, + { 160, 4.71472363599206132240694321176e284, 0 }, + { 161, 7.59070505394721872907517857094e286, 0 }, + { 162, 1.22969421873944943411017892849e289, 0 }, + { 163, 2.00440157654530257759959165344e291, 0 }, + { 164, 3.28721858553429622726333031164e293, 0 }, + { 165, 5.42391066613158877498449501421e295, 0 }, + { 166, 9.00369170577843736647426172359e297, 0 }, + { 167, 1.50361651486499904020120170784e300, 0 }, + { 168, 2.52607574497319838753801886917e302, 0 }, + { 169, 4.26906800900470527493925188890e304, 0 }, + { 170, 7.25741561530799896739672821113e306, 0 }, + + /* + { 171, 1.24101807021766782342484052410e309, 0 }, + { 172, 2.13455108077438865629072570146e311, 0 }, + { 173, 3.69277336973969237538295546352e313, 0 }, + { 174, 6.42542566334706473316634250653e315, 0 }, + { 175, 1.12444949108573632830410993864e318, 0 }, + { 176, 1.97903110431089593781523349201e320, 0 }, + { 177, 3.50288505463028580993296328086e322, 0 }, + { 178, 6.23513539724190874168067463993e324, 0 }, + { 179, 1.11608923610630166476084076055e327, 0 }, + { 180, 2.00896062499134299656951336898e329, 0 }, + { 181, 3.63621873123433082379081919786e331, 0 }, + { 182, 6.61791809084648209929929094011e333, 0 }, + { 183, 1.21107901062490622417177024204e336, 0 }, + { 184, 2.22838537954982745247605724535e338, 0 }, + { 185, 4.12251295216718078708070590390e340, 0 }, + { 186, 7.66787409103095626397011298130e342, 0 }, + { 187, 1.43389245502278882136241112750e345, 0 }, + { 188, 2.69571781544284298416133291969e347, 0 }, + { 189, 5.09490667118697324006491921822e349, 0 }, + { 190, 9.68032267525524915612334651460e351, 0 }, + { 191, 1.84894163097375258881955918429e354, 0 }, + { 192, 3.54996793146960497053355363384e356, 0 }, + { 193, 6.85143810773633759312975851330e358, 0 }, + { 194, 1.32917899290084949306717315158e361, 0 }, + { 195, 2.59189903615665651148098764559e363, 0 }, + { 196, 5.08012211086704676250273578535e365, 0 }, + { 197, 1.00078405584080821221303894971e368, 0 }, + { 198, 1.98155243056480026018181712043e370, 0 }, + { 199, 3.94328933682395251776181606966e372, 0 }, + { 200, 7.88657867364790503552363213932e374, 0 } + */ +}; + +static struct {int n; double f; long i; } doub_fact_table[GSL_SF_DOUBLEFACT_NMAX + 1] = { + { 0, 1.000000000000000000000000000, 1L }, + { 1, 1.000000000000000000000000000, 1L }, + { 2, 2.000000000000000000000000000, 2L }, + { 3, 3.000000000000000000000000000, 3L }, + { 4, 8.000000000000000000000000000, 8L }, + { 5, 15.00000000000000000000000000, 15L }, + { 6, 48.00000000000000000000000000, 48L }, + { 7, 105.0000000000000000000000000, 105L }, + { 8, 384.0000000000000000000000000, 384L }, + { 9, 945.0000000000000000000000000, 945L }, + { 10, 3840.000000000000000000000000, 3840L }, + { 11, 10395.00000000000000000000000, 10395L }, + { 12, 46080.00000000000000000000000, 46080L }, + { 13, 135135.0000000000000000000000, 135135L }, + { 14, 645120.00000000000000000000000, 645120L }, + { 15, 2.02702500000000000000000000000e6, 2027025L }, + { 16, 1.03219200000000000000000000000e7, 10321920L }, + { 17, 3.4459425000000000000000000000e7, 34459425L }, + { 18, 1.85794560000000000000000000000e8, 185794560L }, + { 19, 6.5472907500000000000000000000e8, 0 }, + { 20, 3.7158912000000000000000000000e9, 0 }, + { 21, 1.37493105750000000000000000000e10, 0 }, + { 22, 8.1749606400000000000000000000e10, 0 }, + { 23, 3.1623414322500000000000000000e11, 0 }, + { 24, 1.96199055360000000000000000000e12, 0 }, + { 25, 7.9058535806250000000000000000e12, 0 }, + { 26, 5.1011754393600000000000000000e13, 0 }, + { 27, 2.13458046676875000000000000000e14, 0 }, + { 28, 1.42832912302080000000000000000e15, 0 }, + { 29, 6.1902833536293750000000000000e15, 0 }, + { 30, 4.2849873690624000000000000000e16, 0 }, + { 31, 1.91898783962510625000000000000e17, 0 }, + { 32, 1.37119595809996800000000000000e18, 0 }, + { 33, 6.3326598707628506250000000000e18, 0 }, + { 34, 4.6620662575398912000000000000e19, 0 }, + { 35, 2.21643095476699771875000000000e20, 0 }, + { 36, 1.67834385271436083200000000000e21, 0 }, + { 37, 8.2007945326378915593750000000e21, 0 }, + { 38, 6.3777066403145711616000000000e22, 0 }, + { 39, 3.1983098677287777081562500000e23, 0 }, + { 40, 2.55108265612582846464000000000e24, 0 }, + { 41, 1.31130704576879886034406250000e25, 0 }, + { 42, 1.07145471557284795514880000000e26, 0 }, + { 43, 5.6386202968058350994794687500e26, 0 }, + { 44, 4.7144007485205310026547200000e27, 0 }, + { 45, 2.53737913356262579476576093750e28, 0 }, + { 46, 2.16862434431944426122117120000e29, 0 }, + { 47, 1.19256819277443412353990764062e30, 0 }, + { 48, 1.04093968527333324538616217600e31, 0 }, + { 49, 5.8435841445947272053455474391e31, 0 }, + { 50, 5.2046984263666662269308108800e32, 0 }, + { 51, 2.98022791374331087472622919392e33, 0 }, + { 52, 2.70644318171066643800402165760e34, 0 }, + { 53, 1.57952079428395476360490147278e35, 0 }, + { 54, 1.46147931812375987652217169510e36, 0 }, + { 55, 8.6873643685617511998269581003e36, 0 }, + { 56, 8.1842841814930553085241614926e37, 0 }, + { 57, 4.9517976900801981839013661172e38, 0 }, + { 58, 4.7468848252659720789440136657e39, 0 }, + { 59, 2.92156063714731692850180600912e40, 0 }, + { 60, 2.84813089515958324736640819942e41, 0 }, + { 61, 1.78215198865986332638610166557e42, 0 }, + { 62, 1.76584115499894161336717308364e43, 0 }, + { 63, 1.12275575285571389562324404931e44, 0 }, + { 64, 1.13013833919932263255499077353e45, 0 }, + { 65, 7.2979123935621403215510863205e45, 0 }, + { 66, 7.4589130387155293748629391053e46, 0 }, + { 67, 4.8896013036866340154392278347e47, 0 }, + { 68, 5.0720608663265599749067985916e48, 0 }, + { 69, 3.3738248995437774706530672060e49, 0 }, + { 70, 3.5504426064285919824347590141e50, 0 }, + { 71, 2.39541567867608200416367771623e51, 0 }, + { 72, 2.55631867662858622735302649017e52, 0 }, + { 73, 1.74865344543353986303948473285e53, 0 }, + { 74, 1.89167582070515380824123960272e54, 0 }, + { 75, 1.31149008407515489727961354964e55, 0 }, + { 76, 1.43767362373591689426334209807e56, 0 }, + { 77, 1.00984736473786927090530243322e57, 0 }, + { 78, 1.12138542651401517752540683649e58, 0 }, + { 79, 7.9777941814291672401518892225e58, 0 }, + { 80, 8.9710834121121214202032546920e59, 0 }, + { 81, 6.4620132869576254645230302702e60, 0 }, + { 82, 7.3562883979319395645666688474e61, 0 }, + { 83, 5.3634710281748291355541151243e62, 0 }, + { 84, 6.1792822542628292342360018318e63, 0 }, + { 85, 4.5589503739486047652209978556e64, 0 }, + { 86, 5.3141827386660331414429615754e65, 0 }, + { 87, 3.9662868253352861457422681344e66, 0 }, + { 88, 4.6764808100261091644698061863e67, 0 }, + { 89, 3.5299952745484046697106186396e68, 0 }, + { 90, 4.2088327290234982480228255677e69, 0 }, + { 91, 3.2122956998390482494366629620e70, 0 }, + { 92, 3.8721261107016183881809995223e71, 0 }, + { 93, 2.98743500085031487197609655470e72, 0 }, + { 94, 3.6397985440595212848901395509e73, 0 }, + { 95, 2.83806325080779912837729172696e74, 0 }, + { 96, 3.4942066022971404334945339689e75, 0 }, + { 97, 2.75292135328356515452597297515e76, 0 }, + { 98, 3.4243224702511976248246432895e77, 0 }, + { 99, 2.72539213975072950298071324540e78, 0 }, + { 100, 3.4243224702511976248246432895e79, 0 }, + { 101, 2.75264606114823679801052037785e80, 0 }, + { 102, 3.4928089196562215773211361553e81, 0 }, + { 103, 2.83522544298268390195083598919e82, 0 }, + { 104, 3.6325212764424704404139816015e83, 0 }, + { 105, 2.97698671513181809704837778865e84, 0 }, + { 106, 3.8504725530290186668388204976e85, 0 }, + { 107, 3.1853757851910453638417642339e86, 0 }, + { 108, 4.1585103572713401601859261374e87, 0 }, + { 109, 3.4720596058582394465875230149e88, 0 }, + { 110, 4.5743613929984741762045187512e89, 0 }, + { 111, 3.8539861625026457857121505465e90, 0 }, + { 112, 5.1232847601582910773490610013e91, 0 }, + { 113, 4.3550043636279897378547301176e92, 0 }, + { 114, 5.8405446265804518281779295415e93, 0 }, + { 115, 5.0082550181721881985329396352e94, 0 }, + { 116, 6.7750317668333241206863982681e95, 0 }, + { 117, 5.8596583712614601922835393732e96, 0 }, + { 118, 7.9945374848633224624099499564e97, 0 }, + { 119, 6.9729934618011376288174118541e98, 0 }, + { 120, 9.5934449818359869548919399477e99, 0 }, + { 121, 8.4373220887793765308690683435e100, 0 }, + { 122, 1.17040028778399040849681667362e102, 0 }, + { 123, 1.03779061691986331329689540625e103, 0 }, + { 124, 1.45129635685214810653605267528e104, 0 }, + { 125, 1.29723827114982914162111925781e105, 0 }, + { 126, 1.82863340963370661423542637086e106, 0 }, + { 127, 1.64749260436028300985882145742e107, 0 }, + { 128, 2.34065076433114446622134575470e108, 0 }, + { 129, 2.12526545962476508271787968008e109, 0 }, + { 130, 3.04284599363048780608774948111e110, 0 }, + { 131, 2.78409775210844225836042238090e111, 0 }, + { 132, 4.0165567115922439040358293151e112, 0 }, + { 133, 3.7028500103042282036193617666e113, 0 }, + { 134, 5.3821859935336068314080112822e114, 0 }, + { 135, 4.9988475139107080748861383849e115, 0 }, + { 136, 7.3197729512057052907148953438e116, 0 }, + { 137, 6.8484210940576700625940095873e117, 0 }, + { 138, 1.01012866726638733011865555744e119, 0 }, + { 139, 9.5193053207401613870056733264e119, 0 }, + { 140, 1.41418013417294226216611778042e121, 0 }, + { 141, 1.34222205022436275556779993902e122, 0 }, + { 142, 2.00813579052557801227588724819e123, 0 }, + { 143, 1.91937753182083874046195391280e124, 0 }, + { 144, 2.89171553835683233767727763739e125, 0 }, + { 145, 2.78309742114021617366983317355e126, 0 }, + { 146, 4.2219046860009752130088253506e127, 0 }, + { 147, 4.0911532090761177752946547651e128, 0 }, + { 148, 6.2484189352814433152530615189e129, 0 }, + { 149, 6.0958182815234154851890356000e130, 0 }, + { 150, 9.3726284029221649728795922783e131, 0 }, + { 151, 9.2046856051003573826354437561e132, 0 }, + { 152, 1.42463951724416907587769802630e134, 0 }, + { 153, 1.40831689758035467954322289468e135, 0 }, + { 154, 2.19394485655602037685165496051e136, 0 }, + { 155, 2.18289119124954975329199548675e137, 0 }, + { 156, 3.4225539762273917878885817384e138, 0 }, + { 157, 3.4271391702617931126684329142e139, 0 }, + { 158, 5.4076352824392790248639591467e140, 0 }, + { 159, 5.4491512807162510491428083336e141, 0 }, + { 160, 8.6522164519028464397823346347e142, 0 }, + { 161, 8.7731335619531641891199214170e143, 0 }, + { 162, 1.40165906520826112324473821082e145, 0 }, + { 163, 1.43002077059836576282654719098e146, 0 }, + { 164, 2.29872086694154824212137066574e147, 0 }, + { 165, 2.35953427148730350866380286512e148, 0 }, + { 166, 3.8158766391229700819214753051e149, 0 }, + { 167, 3.9404222333837968594685507847e150, 0 }, + { 168, 6.4106727537265897376280785126e151, 0 }, + { 169, 6.6593135744186166925018508262e152, 0 }, + { 170, 1.08981436813352025539677334714e154, 0 }, + { 171, 1.13874262122558345441781649128e155, 0 }, + { 172, 1.87448071318965483928245015709e156, 0 }, + { 173, 1.97002473472025937614282252992e157, 0 }, + { 174, 3.2615964409499994203514632733e158, 0 }, + { 175, 3.4475432857604539082499394274e159, 0 }, + { 176, 5.7404097360719989798185753611e160, 0 }, + { 177, 6.1021516157960034176023927864e161, 0 }, + { 178, 1.02179293302081581840770641427e163, 0 }, + { 179, 1.09228513922748461175082830877e164, 0 }, + { 180, 1.83922727943746847313387154568e165, 0 }, + { 181, 1.97703610200174714726899923887e166, 0 }, + { 182, 3.3473936485761926211036462131e167, 0 }, + { 183, 3.6179760666631972795022686071e168, 0 }, + { 184, 6.1592043133801944228307090322e169, 0 }, + { 185, 6.6932557233269149670791969232e170, 0 }, + { 186, 1.14561200228871616264651187999e172, 0 }, + { 187, 1.25163882026213309884380982464e173, 0 }, + { 188, 2.15375056430278638577544233437e174, 0 }, + { 189, 2.36559737029543155681480056857e175, 0 }, + { 190, 4.0921260721752941329733404353e176, 0 }, + { 191, 4.5182909772642742735162690860e177, 0 }, + { 192, 7.8568820585765647353088136358e178, 0 }, + { 193, 8.7203015861200493478863993359e179, 0 }, + { 194, 1.52423511936385355864990984535e181, 0 }, + { 195, 1.70045880929340962283784787050e182, 0 }, + { 196, 2.98750083395315297495382329688e183, 0 }, + { 197, 3.3499038543080169569905603049e184, 0 }, + { 198, 5.9152516512272428904085701278e185, 0 }, + { 199, 6.6663086700729537444112150067e186, 0 }, + { 200, 1.18305033024544857808171402556e188, 0 }, + { 201, 1.33992804268466370262665421635e189, 0 }, + { 202, 2.38976166709580612772506233164e190, 0 }, + { 203, 2.72005392664986731633210805920e191, 0 }, + { 204, 4.8751138008754445005591271565e192, 0 }, + { 205, 5.5761105496322279984808215214e193, 0 }, + { 206, 1.00427344298034156711518019425e195, 0 }, + { 207, 1.15425488377387119568553005492e196, 0 }, + { 208, 2.08888876139911045959957480403e197, 0 }, + { 209, 2.41239270708739079898275781478e198, 0 }, + { 210, 4.3866663989381319651591070885e199, 0 }, + { 211, 5.0901486119543945858536189892e200, 0 }, + { 212, 9.2997327657488397661373070276e201, 0 }, + { 213, 1.08420165434628604678682084470e203, 0 }, + { 214, 1.99014281187025170995338370390e204, 0 }, + { 215, 2.33103355684451500059166481610e205, 0 }, + { 216, 4.2987084736397436934993088004e206, 0 }, + { 217, 5.0583428183525975512839126509e207, 0 }, + { 218, 9.3711844725346412518284931849e208, 0 }, + { 219, 1.10777707721921886373117687056e210, 0 }, + { 220, 2.06166058395762107540226850068e211, 0 }, + { 221, 2.44818734065447368884590088393e212, 0 }, + { 222, 4.5768864963859187873930360715e213, 0 }, + { 223, 5.4594577696594763261263589712e214, 0 }, + { 224, 1.02522257519044580837604008002e216, 0 }, + { 225, 1.22837799817338217337843076851e217, 0 }, + { 226, 2.31700301993040752692985058084e218, 0 }, + { 227, 2.78841805585357753356903784452e219, 0 }, + { 228, 5.2827668854413291614000593243e220, 0 }, + { 229, 6.3854773479046925518730966640e221, 0 }, + { 230, 1.21503638365150570712201364459e223, 0 }, + { 231, 1.47504526736598397948268532937e224, 0 }, + { 232, 2.81888441007149324052307165546e225, 0 }, + { 233, 3.4368554729627426721946568174e226, 0 }, + { 234, 6.5961895195672941828239876738e227, 0 }, + { 235, 8.0766103614624452796574435210e228, 0 }, + { 236, 1.55670072661788142714646109101e230, 0 }, + { 237, 1.91415665566659953127881411447e231, 0 }, + { 238, 3.7049477293505577966085773966e232, 0 }, + { 239, 4.5748344070431728797563657336e233, 0 }, + { 240, 8.8918745504413387118605857518e234, 0 }, + { 241, 1.10253509209740466402128414180e236, 0 }, + { 242, 2.15183364120680396827026175195e237, 0 }, + { 243, 2.67916027379669333357172046456e238, 0 }, + { 244, 5.2504740845446016825794386748e239, 0 }, + { 245, 6.5639426708018986672507151382e240, 0 }, + { 246, 1.29161662479797201391454191399e242, 0 }, + { 247, 1.62129383968806897081092663913e243, 0 }, + { 248, 3.2032092294989705945080639467e244, 0 }, + { 249, 4.0370216608232917373192073314e245, 0 }, + { 250, 8.0080230737474264862701598667e246, 0 }, + { 251, 1.01329243686664622606712104019e248, 0 }, + { 252, 2.01802181458435147454008028642e249, 0 }, + { 253, 2.56362986527261495194981623168e250, 0 }, + { 254, 5.1257754090442527453318039275e251, 0 }, + { 255, 6.5372561564451681274720313908e252, 0 }, + { 256, 1.31219850471532870280494180544e254, 0 }, + { 257, 1.68007483220640820876031206743e255, 0 }, + { 258, 3.3854721421655480532367498580e256, 0 }, + { 259, 4.3513938154145972606892082546e257, 0 }, + { 260, 8.8022275696304249384155496309e258, 0 }, + { 261, 1.13571378582320988503988335446e260, 0 }, + { 262, 2.30618362324317133386487400329e261, 0 }, + { 263, 2.98692725671504199765489322224e262, 0 }, + { 264, 6.0883247653619723214032673687e263, 0 }, + { 265, 7.9153572302948612937854670389e264, 0 }, + { 266, 1.61949438758628463749326912007e266, 0 }, + { 267, 2.11340038048872796544071969939e267, 0 }, + { 268, 4.3402449587312428284819612418e268, 0 }, + { 269, 5.6850470235146782270355359914e269, 0 }, + { 270, 1.17186613885743556369012953528e271, 0 }, + { 271, 1.54064774337247779952663025366e272, 0 }, + { 272, 3.1874758976922247332371523360e273, 0 }, + { 273, 4.2059683394068643927077005925e274, 0 }, + { 274, 8.7336839596766957690697974006e275, 0 }, + { 275, 1.15664129333688770799461766294e277, 0 }, + { 276, 2.41049677287076803226326408256e278, 0 }, + { 277, 3.2038963825431789511450909263e279, 0 }, + { 278, 6.7011810285807351296918741495e280, 0 }, + { 279, 8.9388709072954692736948036845e281, 0 }, + { 280, 1.87633068800260583631372476186e283, 0 }, + { 281, 2.51182272495002686590823983534e284, 0 }, + { 282, 5.2912525401673484584047038284e285, 0 }, + { 283, 7.1084583116085760305203187340e286, 0 }, + { 284, 1.50271572140752696218693588728e288, 0 }, + { 285, 2.02591061880844416869829083919e289, 0 }, + { 286, 4.2977669632255271118546366376e290, 0 }, + { 287, 5.8143634759802347641640947085e291, 0 }, + { 288, 1.23775688540895180821413535163e293, 0 }, + { 289, 1.68035104455828784684342337075e294, 0 }, + { 290, 3.5894949676859602438209925197e295, 0 }, + { 291, 4.8898215396646176343143620089e296, 0 }, + { 292, 1.04813253056430039119572981576e298, 0 }, + { 293, 1.43271771112173296685410806860e299, 0 }, + { 294, 3.08150963985904315011544565835e300, 0 }, + { 295, 4.2265172478091122522196188024e301, 0 }, + { 296, 9.1212685339827677243417191487e302, 0 }, + { 297, 1.25527562259930633890922678431e304, 0 }, + /* + { 298, 2.71813802312686478185383230631e305, 0 }, + { 299, 3.7532741115719259533385880851e306, 0 }, + { 300, 8.1544140693805943455614969189e307, } + */ +}; + + +/* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1 + */ +static double gstar_a_data[30] = { + 2.16786447866463034423060819465, + -0.05533249018745584258035832802, + 0.01800392431460719960888319748, + -0.00580919269468937714480019814, + 0.00186523689488400339978881560, + -0.00059746524113955531852595159, + 0.00019125169907783353925426722, + -0.00006124996546944685735909697, + 0.00001963889633130842586440945, + -6.3067741254637180272515795142e-06, + 2.0288698405861392526872789863e-06, + -6.5384896660838465981983750582e-07, + 2.1108698058908865476480734911e-07, + -6.8260714912274941677892994580e-08, + 2.2108560875880560555583978510e-08, + -7.1710331930255456643627187187e-09, + 2.3290892983985406754602564745e-09, + -7.5740371598505586754890405359e-10, + 2.4658267222594334398525312084e-10, + -8.0362243171659883803428749516e-11, + 2.6215616826341594653521346229e-11, + -8.5596155025948750540420068109e-12, + 2.7970831499487963614315315444e-12, + -9.1471771211886202805502562414e-13, + 2.9934720198063397094916415927e-13, + -9.8026575909753445931073620469e-14, + 3.2116773667767153777571410671e-14, + -1.0518035333878147029650507254e-14, + 3.4144405720185253938994854173e-15, + -1.0115153943081187052322643819e-15 +}; +static cheb_series gstar_a_cs = { + gstar_a_data, + 29, + -1, 1, + 17 +}; + + +/* Chebyshev coefficients for + * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1 + */ +static double gstar_b_data[] = { + 0.0057502277273114339831606096782, + 0.0004496689534965685038254147807, + -0.0001672763153188717308905047405, + 0.0000615137014913154794776670946, + -0.0000223726551711525016380862195, + 8.0507405356647954540694800545e-06, + -2.8671077107583395569766746448e-06, + 1.0106727053742747568362254106e-06, + -3.5265558477595061262310873482e-07, + 1.2179216046419401193247254591e-07, + -4.1619640180795366971160162267e-08, + 1.4066283500795206892487241294e-08, + -4.6982570380537099016106141654e-09, + 1.5491248664620612686423108936e-09, + -5.0340936319394885789686867772e-10, + 1.6084448673736032249959475006e-10, + -5.0349733196835456497619787559e-11, + 1.5357154939762136997591808461e-11, + -4.5233809655775649997667176224e-12, + 1.2664429179254447281068538964e-12, + -3.2648287937449326771785041692e-13, + 7.1528272726086133795579071407e-14, + -9.4831735252566034505739531258e-15, + -2.3124001991413207293120906691e-15, + 2.8406613277170391482590129474e-15, + -1.7245370321618816421281770927e-15, + 8.6507923128671112154695006592e-16, + -3.9506563665427555895391869919e-16, + 1.6779342132074761078792361165e-16, + -6.0483153034414765129837716260e-17 +}; +static cheb_series gstar_b_cs = { + gstar_b_data, + 29, + -1, 1, + 18 +}; + + +/* coefficients for gamma=7, kmax=8 Lanczos method */ +static double lanczos_7_c[9] = { + 0.99999999999980993227684700473478, + 676.520368121885098567009190444019, + -1259.13921672240287047156078755283, + 771.3234287776530788486528258894, + -176.61502916214059906584551354, + 12.507343278686904814458936853, + -0.13857109526572011689554707, + 9.984369578019570859563e-6, + 1.50563273514931155834e-7 +}; + +/* complex version of Lanczos method; this is not safe for export + * since it becomes bad in the left half-plane + */ +static +int +lngamma_lanczos_complex(double zr, double zi, gsl_sf_result * yr, gsl_sf_result * yi) +{ + int k; + gsl_sf_result log1_r, log1_i; + gsl_sf_result logAg_r, logAg_i; + double Ag_r, Ag_i; + double yi_tmp_val, yi_tmp_err; + + zr -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ + + Ag_r = lanczos_7_c[0]; + Ag_i = 0.0; + for(k=1; k<=8; k++) { + double R = zr + k; + double I = zi; + double a = lanczos_7_c[k] / (R*R + I*I); + Ag_r += a * R; + Ag_i -= a * I; + } + + gsl_sf_complex_log_e(zr + 7.5, zi, &log1_r, &log1_i); + gsl_sf_complex_log_e(Ag_r, Ag_i, &logAg_r, &logAg_i); + + /* (z+0.5)*log(z+7.5) - (z+7.5) + LogRootTwoPi_ + log(Ag(z)) */ + yr->val = (zr+0.5)*log1_r.val - zi*log1_i.val - (zr+7.5) + LogRootTwoPi_ + logAg_r.val; + yi->val = zi*log1_r.val + (zr+0.5)*log1_i.val - zi + logAg_i.val; + yr->err = 4.0 * GSL_DBL_EPSILON * fabs(yr->val); + yi->err = 4.0 * GSL_DBL_EPSILON * fabs(yi->val); + yi_tmp_val = yi->val; + yi_tmp_err = yi->err; + gsl_sf_angle_restrict_symm_err_e(yi_tmp_val, yi); + yi->err += yi_tmp_err; + return GSL_SUCCESS; +} + + +/* Lanczos method for real x > 0; + * gamma=7, truncated at 1/(z+8) + * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86] + */ +static +int +lngamma_lanczos(double x, gsl_sf_result * result) +{ + int k; + double Ag; + double term1, term2; + + x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ + + Ag = lanczos_7_c[0]; + for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); } + + /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */ + term1 = (x+0.5)*log((x+7.5)/M_E); + term2 = LogRootTwoPi_ + log(Ag); + result->val = term1 + (term2 - 7.0); + result->err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + 7.0); + result->err += GSL_DBL_EPSILON * fabs(result->val); + + return GSL_SUCCESS; +} + +/* x = eps near zero + * gives double-precision for |eps| < 0.02 + */ +static +int +lngamma_sgn_0(double eps, gsl_sf_result * lng, double * sgn) +{ + /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */ + const double c1 = -0.07721566490153286061; + const double c2 = -0.01094400467202744461; + const double c3 = 0.09252092391911371098; + const double c4 = -0.01827191316559981266; + const double c5 = 0.01800493109685479790; + const double c6 = -0.00685088537872380685; + const double c7 = 0.00399823955756846603; + const double c8 = -0.00189430621687107802; + const double c9 = 0.00097473237804513221; + const double c10 = -0.00048434392722255893; + const double g6 = c6+eps*(c7+eps*(c8 + eps*(c9 + eps*c10))); + const double g = eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*g6))))); + + /* calculate Gamma(eps) eps, a positive quantity */ + const double gee = g + 1.0/(1.0+eps) + 0.5*eps; + + lng->val = log(gee/fabs(eps)); + lng->err = 4.0 * GSL_DBL_EPSILON * fabs(lng->val); + *sgn = GSL_SIGN(eps); + + return GSL_SUCCESS; +} + + +/* x near a negative integer + * Calculates sign as well as log(|gamma(x)|). + * x = -N + eps + * assumes N >= 1 + */ +static +int +lngamma_sgn_sing(int N, double eps, gsl_sf_result * lng, double * sgn) +{ + if(eps == 0.0) { + lng->val = 0.0; + lng->err = 0.0; + *sgn = 0.0; + GSL_ERROR ("error", GSL_EDOM); + } + else if(N == 1) { + /* calculate series for + * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2) + * double-precision for |eps| < 0.02 + */ + const double c0 = 0.07721566490153286061; + const double c1 = 0.08815966957356030521; + const double c2 = -0.00436125434555340577; + const double c3 = 0.01391065882004640689; + const double c4 = -0.00409427227680839100; + const double c5 = 0.00275661310191541584; + const double c6 = -0.00124162645565305019; + const double c7 = 0.00065267976121802783; + const double c8 = -0.00032205261682710437; + const double c9 = 0.00016229131039545456; + const double g5 = c5 + eps*(c6 + eps*(c7 + eps*(c8 + eps*c9))); + const double g = eps*(c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*g5))))); + + /* calculate eps gamma(-1+eps), a negative quantity */ + const double gam_e = g - 1.0 - 0.5*eps*(1.0+3.0*eps)/(1.0 - eps*eps); + + lng->val = log(fabs(gam_e)/fabs(eps)); + lng->err = 2.0 * GSL_DBL_EPSILON * fabs(lng->val); + *sgn = ( eps > 0.0 ? -1.0 : 1.0 ); + return GSL_SUCCESS; + } + else { + double g; + + /* series for sin(Pi(N+1-eps))/(Pi eps) modulo the sign + * double-precision for |eps| < 0.02 + */ + const double cs1 = -1.6449340668482264365; + const double cs2 = 0.8117424252833536436; + const double cs3 = -0.1907518241220842137; + const double cs4 = 0.0261478478176548005; + const double cs5 = -0.0023460810354558236; + const double e2 = eps*eps; + const double sin_ser = 1.0 + e2*(cs1+e2*(cs2+e2*(cs3+e2*(cs4+e2*cs5)))); + + /* calculate series for ln(gamma(1+N-eps)) + * double-precision for |eps| < 0.02 + */ + double aeps = fabs(eps); + double c1, c2, c3, c4, c5, c6, c7; + double lng_ser; + gsl_sf_result c0; + gsl_sf_result psi_0; + gsl_sf_result psi_1; + gsl_sf_result psi_2; + gsl_sf_result psi_3; + gsl_sf_result psi_4; + gsl_sf_result psi_5; + gsl_sf_result psi_6; + psi_2.val = 0.0; + psi_3.val = 0.0; + psi_4.val = 0.0; + psi_5.val = 0.0; + psi_6.val = 0.0; + gsl_sf_lnfact_e(N, &c0); + gsl_sf_psi_int_e(N+1, &psi_0); + gsl_sf_psi_1_int_e(N+1, &psi_1); + if(aeps > 0.00001) gsl_sf_psi_n_e(2, N+1.0, &psi_2); + if(aeps > 0.0002) gsl_sf_psi_n_e(3, N+1.0, &psi_3); + if(aeps > 0.001) gsl_sf_psi_n_e(4, N+1.0, &psi_4); + if(aeps > 0.005) gsl_sf_psi_n_e(5, N+1.0, &psi_5); + if(aeps > 0.01) gsl_sf_psi_n_e(6, N+1.0, &psi_6); + c1 = psi_0.val; + c2 = psi_1.val/2.0; + c3 = psi_2.val/6.0; + c4 = psi_3.val/24.0; + c5 = psi_4.val/120.0; + c6 = psi_5.val/720.0; + c7 = psi_6.val/5040.0; + lng_ser = c0.val-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7)))))); + + /* calculate + * g = ln(|eps gamma(-N+eps)|) + * = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|) + */ + g = -lng_ser - log(sin_ser); + + lng->val = g - log(fabs(eps)); + lng->err = c0.err + 2.0 * GSL_DBL_EPSILON * (fabs(g) + fabs(lng->val)); + + *sgn = ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 ); + + return GSL_SUCCESS; + } +} + + +/* This gets bad near the negative half axis. However, this + * region can be avoided by use of the reflection formula, as usual. + * Only the first two terms of the series are kept. + */ +#if 0 +static +int +lngamma_complex_stirling(const double zr, const double zi, double * lg_r, double * arg) +{ + double re_zinv, im_zinv; + double re_zinv2, im_zinv2; + double re_zinv3, im_zinv3; + double re_zhlnz, im_zhlnz; + double r, lnr, theta; + gsl_sf_complex_log_e(zr, zi, &lnr, &theta); /* z = r e^{i theta} */ + r = exp(lnr); + re_zinv = (zr/r)/r; + im_zinv = -(zi/r)/r; + re_zinv2 = re_zinv*re_zinv - im_zinv*im_zinv; + re_zinv2 = 2.0*re_zinv*im_zinv; + re_zinv3 = re_zinv2*re_zinv - im_zinv2*im_zinv; + re_zinv3 = re_zinv2*im_zinv + im_zinv2*re_zinv; + re_zhlnz = (zr - 0.5)*lnr - zi*theta; + im_zhlnz = zi*lnr + zr*theta; + *lg_r = re_zhlnz - zr + 0.5*(M_LN2+M_LNPI) + re_zinv/12.0 - re_zinv3/360.0; + *arg = im_zhlnz - zi + 1.0/12.0*im_zinv - im_zinv3/360.0; + return GSL_SUCCESS; +} +#endif /* 0 */ + + +inline +static +int +lngamma_1_pade(const double eps, gsl_sf_result * result) +{ + /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps + * plus a correction series. + */ + const double n1 = -1.0017419282349508699871138440; + const double n2 = 1.7364839209922879823280541733; + const double d1 = 1.2433006018858751556055436011; + const double d2 = 5.0456274100274010152489597514; + const double num = (eps + n1) * (eps + n2); + const double den = (eps + d1) * (eps + d2); + const double pade = 2.0816265188662692474880210318 * num / den; + const double c0 = 0.004785324257581753; + const double c1 = -0.01192457083645441; + const double c2 = 0.01931961413960498; + const double c3 = -0.02594027398725020; + const double c4 = 0.03141928755021455; + const double eps5 = eps*eps*eps*eps*eps; + const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); + result->val = eps * (pade + corr); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; +} + +inline +static +int +lngamma_2_pade(const double eps, gsl_sf_result * result) +{ + /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps + * plus a correction series. + */ + const double n1 = 1.000895834786669227164446568; + const double n2 = 4.209376735287755081642901277; + const double d1 = 2.618851904903217274682578255; + const double d2 = 10.85766559900983515322922936; + const double num = (eps + n1) * (eps + n2); + const double den = (eps + d1) * (eps + d2); + const double pade = 2.85337998765781918463568869 * num/den; + const double c0 = 0.0001139406357036744; + const double c1 = -0.0001365435269792533; + const double c2 = 0.0001067287169183665; + const double c3 = -0.0000693271800931282; + const double c4 = 0.0000407220927867950; + const double eps5 = eps*eps*eps*eps*eps; + const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); + result->val = eps * (pade + corr); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; +} + + +/* series for gammastar(x) + * double-precision for x > 10.0 + */ +static +int +gammastar_ser(const double x, gsl_sf_result * result) +{ + /* Use the Stirling series for the correction to Log(Gamma(x)), + * which is better behaved and easier to compute than the + * regular Stirling series for Gamma(x). + */ + const double y = 1.0/(x*x); + const double c0 = 1.0/12.0; + const double c1 = -1.0/360.0; + const double c2 = 1.0/1260.0; + const double c3 = -1.0/1680.0; + const double c4 = 1.0/1188.0; + const double c5 = -691.0/360360.0; + const double c6 = 1.0/156.0; + const double c7 = -3617.0/122400.0; + const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7)))))); + result->val = exp(ser/x); + result->err = 2.0 * GSL_DBL_EPSILON * result->val * GSL_MAX_DBL(1.0, ser/x); + return GSL_SUCCESS; +} + + +/* Chebyshev expansion for log(gamma(x)/gamma(8)) + * 5 < x < 10 + * -1 < t < 1 + */ +static double gamma_5_10_data[24] = { + -1.5285594096661578881275075214, + 4.8259152300595906319768555035, + 0.2277712320977614992970601978, + -0.0138867665685617873604917300, + 0.0012704876495201082588139723, + -0.0001393841240254993658962470, + 0.0000169709242992322702260663, + -2.2108528820210580075775889168e-06, + 3.0196602854202309805163918716e-07, + -4.2705675000079118380587357358e-08, + 6.2026423818051402794663551945e-09, + -9.1993973208880910416311405656e-10, + 1.3875551258028145778301211638e-10, + -2.1218861491906788718519522978e-11, + 3.2821736040381439555133562600e-12, + -5.1260001009953791220611135264e-13, + 8.0713532554874636696982146610e-14, + -1.2798522376569209083811628061e-14, + 2.0417711600852502310258808643e-15, + -3.2745239502992355776882614137e-16, + 5.2759418422036579482120897453e-17, + -8.5354147151695233960425725513e-18, + 1.3858639703888078291599886143e-18, + -2.2574398807738626571560124396e-19 +}; +static const cheb_series gamma_5_10_cs = { + gamma_5_10_data, + 23, + -1, 1, + 11 +}; + + +/* gamma(x) for x >= 1/2 + * assumes x >= 1/2 + */ +static +int +gamma_xgthalf(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x == 0.5) { + result->val = 1.77245385090551602729817; + result->err = GSL_DBL_EPSILON * result->val; + return GSL_SUCCESS; + } else if (x <= (GSL_SF_FACT_NMAX + 1.0) && x == floor(x)) { + int n = (int) floor (x); + result->val = fact_table[n - 1].f; + result->err = GSL_DBL_EPSILON * result->val; + return GSL_SUCCESS; + } + else if(fabs(x - 1.0) < 0.01) { + /* Use series for Gamma[1+eps] - 1/(1+eps). + */ + const double eps = x - 1.0; + const double c1 = 0.4227843350984671394; + const double c2 = -0.01094400467202744461; + const double c3 = 0.09252092391911371098; + const double c4 = -0.018271913165599812664; + const double c5 = 0.018004931096854797895; + const double c6 = -0.006850885378723806846; + const double c7 = 0.003998239557568466030; + result->val = 1.0/x + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*c7)))))); + result->err = GSL_DBL_EPSILON; + return GSL_SUCCESS; + } + else if(fabs(x - 2.0) < 0.01) { + /* Use series for Gamma[1 + eps]. + */ + const double eps = x - 2.0; + const double c1 = 0.4227843350984671394; + const double c2 = 0.4118403304264396948; + const double c3 = 0.08157691924708626638; + const double c4 = 0.07424901075351389832; + const double c5 = -0.00026698206874501476832; + const double c6 = 0.011154045718130991049; + const double c7 = -0.002852645821155340816; + const double c8 = 0.0021039333406973880085; + result->val = 1.0 + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8))))))); + result->err = GSL_DBL_EPSILON; + return GSL_SUCCESS; + } + else if(x < 5.0) { + /* Exponentiating the logarithm is fine, as + * long as the exponential is not so large + * that it greatly amplifies the error. + */ + gsl_sf_result lg; + lngamma_lanczos(x, &lg); + result->val = exp(lg.val); + result->err = result->val * (lg.err + 2.0 * GSL_DBL_EPSILON); + return GSL_SUCCESS; + } + else if(x < 10.0) { + /* This is a sticky area. The logarithm + * is too large and the gammastar series + * is not good. + */ + const double gamma_8 = 5040.0; + const double t = (2.0*x - 15.0)/5.0; + gsl_sf_result c; + cheb_eval_e(&gamma_5_10_cs, t, &c); + result->val = exp(c.val) * gamma_8; + result->err = result->val * c.err; + result->err += 2.0 * GSL_DBL_EPSILON * result->val; + return GSL_SUCCESS; + } + else if(x < GSL_SF_GAMMA_XMAX) { + /* We do not want to exponentiate the logarithm + * if x is large because of the inevitable + * inflation of the error. So we carefully + * use pow() and exp() with exact quantities. + */ + double p = pow(x, 0.5*x); + double e = exp(-x); + double q = (p * e) * p; + double pre = M_SQRT2 * M_SQRTPI * q/sqrt(x); + gsl_sf_result gstar; + int stat_gs = gammastar_ser(x, &gstar); + result->val = pre * gstar.val; + result->err = (x + 2.5) * GSL_DBL_EPSILON * result->val; + return stat_gs; + } + else { + OVERFLOW_ERROR(result); + } +} + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + + +int gsl_sf_lngamma_e(double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(fabs(x - 1.0) < 0.01) { + /* Note that we must amplify the errors + * from the Pade evaluations because of + * the way we must pass the argument, i.e. + * writing (1-x) is a loss of precision + * when x is near 1. + */ + int stat = lngamma_1_pade(x - 1.0, result); + result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0)); + return stat; + } + else if(fabs(x - 2.0) < 0.01) { + int stat = lngamma_2_pade(x - 2.0, result); + result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0)); + return stat; + } + else if(x >= 0.5) { + return lngamma_lanczos(x, result); + } + else if(x == 0.0) { + DOMAIN_ERROR(result); + } + else if(fabs(x) < 0.02) { + double sgn; + return lngamma_sgn_0(x, result, &sgn); + } + else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) { + /* Try to extract a fractional + * part from x. + */ + double z = 1.0 - x; + double s = sin(M_PI*z); + double as = fabs(s); + if(s == 0.0) { + DOMAIN_ERROR(result); + } + else if(as < M_PI*0.015) { + /* x is near a negative integer, -N */ + if(x < INT_MIN + 2.0) { + result->val = 0.0; + result->err = 0.0; + GSL_ERROR ("error", GSL_EROUND); + } + else { + int N = -(int)(x - 0.5); + double eps = x + N; + double sgn; + return lngamma_sgn_sing(N, eps, result, &sgn); + } + } + else { + gsl_sf_result lg_z; + lngamma_lanczos(z, &lg_z); + result->val = M_LNPI - (log(as) + lg_z.val); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg_z.err; + return GSL_SUCCESS; + } + } + else { + /* |x| was too large to extract any fractional part */ + result->val = 0.0; + result->err = 0.0; + GSL_ERROR ("error", GSL_EROUND); + } +} + + +int gsl_sf_lngamma_sgn_e(double x, gsl_sf_result * result_lg, double * sgn) +{ + if(fabs(x - 1.0) < 0.01) { + int stat = lngamma_1_pade(x - 1.0, result_lg); + result_lg->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0)); + *sgn = 1.0; + return stat; + } + else if(fabs(x - 2.0) < 0.01) { + int stat = lngamma_2_pade(x - 2.0, result_lg); + result_lg->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0)); + *sgn = 1.0; + return stat; + } + else if(x >= 0.5) { + *sgn = 1.0; + return lngamma_lanczos(x, result_lg); + } + else if(x == 0.0) { + *sgn = 0.0; + DOMAIN_ERROR(result_lg); + } + else if(fabs(x) < 0.02) { + return lngamma_sgn_0(x, result_lg, sgn); + } + else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) { + /* Try to extract a fractional + * part from x. + */ + double z = 1.0 - x; + double s = sin(M_PI*x); + double as = fabs(s); + if(s == 0.0) { + *sgn = 0.0; + DOMAIN_ERROR(result_lg); + } + else if(as < M_PI*0.015) { + /* x is near a negative integer, -N */ + if(x < INT_MIN + 2.0) { + result_lg->val = 0.0; + result_lg->err = 0.0; + *sgn = 0.0; + GSL_ERROR ("error", GSL_EROUND); + } + else { + int N = -(int)(x - 0.5); + double eps = x + N; + return lngamma_sgn_sing(N, eps, result_lg, sgn); + } + } + else { + gsl_sf_result lg_z; + lngamma_lanczos(z, &lg_z); + *sgn = (s > 0.0 ? 1.0 : -1.0); + result_lg->val = M_LNPI - (log(as) + lg_z.val); + result_lg->err = 2.0 * GSL_DBL_EPSILON * fabs(result_lg->val) + lg_z.err; + return GSL_SUCCESS; + } + } + else { + /* |x| was too large to extract any fractional part */ + result_lg->val = 0.0; + result_lg->err = 0.0; + *sgn = 0.0; + GSL_ERROR ("error", GSL_EROUND); + } +} + + +int +gsl_sf_gamma_e(const double x, gsl_sf_result * result) +{ + if(x < 0.5) { + int rint_x = (int)floor(x+0.5); + double f_x = x - rint_x; + double sgn_gamma = ( GSL_IS_EVEN(rint_x) ? 1.0 : -1.0 ); + double sin_term = sgn_gamma * sin(M_PI * f_x) / M_PI; + + if(sin_term == 0.0) { + DOMAIN_ERROR(result); + } + else if(x > -169.0) { + gsl_sf_result g; + gamma_xgthalf(1.0-x, &g); + if(fabs(sin_term) * g.val * GSL_DBL_MIN < 1.0) { + result->val = 1.0/(sin_term * g.val); + result->err = fabs(g.err/g.val) * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + UNDERFLOW_ERROR(result); + } + } + else { + /* It is hard to control it here. + * We can only exponentiate the + * logarithm and eat the loss of + * precision. + */ + gsl_sf_result lng; + double sgn; + int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn); + int stat_e = gsl_sf_exp_mult_err_e(lng.val, lng.err, sgn, 0.0, result); + return GSL_ERROR_SELECT_2(stat_e, stat_lng); + } + } + else { + return gamma_xgthalf(x, result); + } +} + + +int +gsl_sf_gammastar_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x < 0.5) { + gsl_sf_result lg; + const int stat_lg = gsl_sf_lngamma_e(x, &lg); + const double lx = log(x); + const double c = 0.5*(M_LN2+M_LNPI); + const double lnr_val = lg.val - (x-0.5)*lx + x - c; + const double lnr_err = lg.err + 2.0 * GSL_DBL_EPSILON *((x+0.5)*fabs(lx) + c); + const int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, result); + return GSL_ERROR_SELECT_2(stat_lg, stat_e); + } + else if(x < 2.0) { + const double t = 4.0/3.0*(x-0.5) - 1.0; + return cheb_eval_e(&gstar_a_cs, t, result); + } + else if(x < 10.0) { + const double t = 0.25*(x-2.0) - 1.0; + gsl_sf_result c; + cheb_eval_e(&gstar_b_cs, t, &c); + result->val = c.val/(x*x) + 1.0 + 1.0/(12.0*x); + result->err = c.err/(x*x); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x < 1.0/GSL_ROOT4_DBL_EPSILON) { + return gammastar_ser(x, result); + } + else if(x < 1.0/GSL_DBL_EPSILON) { + /* Use Stirling formula for Gamma(x). + */ + const double xi = 1.0/x; + result->val = 1.0 + xi/12.0*(1.0 + xi/24.0*(1.0 - xi*(139.0/180.0 + 571.0/8640.0*xi))); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + result->val = 1.0; + result->err = 1.0/x; + return GSL_SUCCESS; + } +} + + +int +gsl_sf_gammainv_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if (x <= 0.0 && x == floor(x)) { /* negative integer */ + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } else if(x < 0.5) { + gsl_sf_result lng; + double sgn; + int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn); + if(stat_lng == GSL_EDOM) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(stat_lng != GSL_SUCCESS) { + result->val = 0.0; + result->err = 0.0; + return stat_lng; + } + else { + return gsl_sf_exp_mult_err_e(-lng.val, lng.err, sgn, 0.0, result); + } + } + else { + gsl_sf_result g; + int stat_g = gamma_xgthalf(x, &g); + if(stat_g == GSL_EOVRFLW) { + UNDERFLOW_ERROR(result); + } + else { + result->val = 1.0/g.val; + result->err = fabs(g.err/g.val) * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + CHECK_UNDERFLOW(result); + return GSL_SUCCESS; + } + } +} + + +int +gsl_sf_lngamma_complex_e(double zr, double zi, gsl_sf_result * lnr, gsl_sf_result * arg) +{ + if(zr <= 0.5) { + /* Transform to right half plane using reflection; + * in fact we do a little better by stopping at 1/2. + */ + double x = 1.0-zr; + double y = -zi; + gsl_sf_result a, b; + gsl_sf_result lnsin_r, lnsin_i; + + int stat_l = lngamma_lanczos_complex(x, y, &a, &b); + int stat_s = gsl_sf_complex_logsin_e(M_PI*zr, M_PI*zi, &lnsin_r, &lnsin_i); + + if(stat_s == GSL_SUCCESS) { + int stat_r; + lnr->val = M_LNPI - lnsin_r.val - a.val; + lnr->err = lnsin_r.err + a.err + 2.0 * GSL_DBL_EPSILON * fabs(lnr->val); + arg->val = -lnsin_i.val - b.val; + arg->err = lnsin_i.err + b.err + 2.0 * GSL_DBL_EPSILON * fabs(arg->val); + stat_r = gsl_sf_angle_restrict_symm_e(&(arg->val)); + return GSL_ERROR_SELECT_2(stat_r, stat_l); + } + else { + DOMAIN_ERROR_2(lnr,arg); + } + } + else { + /* otherwise plain vanilla Lanczos */ + return lngamma_lanczos_complex(zr, zi, lnr, arg); + } +} + + +int gsl_sf_taylorcoeff_e(const int n, const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x < 0.0 || n < 0) { + DOMAIN_ERROR(result); + } + else if(n == 0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(n == 1) { + result->val = x; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(x == 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else { + const double log2pi = M_LNPI + M_LN2; + const double ln_test = n*(log(x)+1.0) + 1.0 - (n+0.5)*log(n+1.0) + 0.5*log2pi; + + if(ln_test < GSL_LOG_DBL_MIN+1.0) { + UNDERFLOW_ERROR(result); + } + else if(ln_test > GSL_LOG_DBL_MAX-1.0) { + OVERFLOW_ERROR(result); + } + else { + double product = 1.0; + int k; + for(k=1; k<=n; k++) { + product *= (x/k); + } + result->val = product; + result->err = n * GSL_DBL_EPSILON * product; + CHECK_UNDERFLOW(result); + return GSL_SUCCESS; + } + } +} + + +int gsl_sf_fact_e(const unsigned int n, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(n < 18) { + result->val = fact_table[n].f; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(n <= GSL_SF_FACT_NMAX){ + result->val = fact_table[n].f; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + + +int gsl_sf_doublefact_e(const unsigned int n, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(n < 26) { + result->val = doub_fact_table[n].f; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(n <= GSL_SF_DOUBLEFACT_NMAX){ + result->val = doub_fact_table[n].f; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + + +int gsl_sf_lnfact_e(const unsigned int n, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(n <= GSL_SF_FACT_NMAX){ + result->val = log(fact_table[n].f); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + gsl_sf_lngamma_e(n+1.0, result); + return GSL_SUCCESS; + } +} + + +int gsl_sf_lndoublefact_e(const unsigned int n, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(n <= GSL_SF_DOUBLEFACT_NMAX){ + result->val = log(doub_fact_table[n].f); + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(GSL_IS_ODD(n)) { + gsl_sf_result lg; + gsl_sf_lngamma_e(0.5*(n+2.0), &lg); + result->val = 0.5*(n+1.0) * M_LN2 - 0.5*M_LNPI + lg.val; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err; + return GSL_SUCCESS; + } + else { + gsl_sf_result lg; + gsl_sf_lngamma_e(0.5*n+1.0, &lg); + result->val = 0.5*n*M_LN2 + lg.val; + result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err; + return GSL_SUCCESS; + } +} + + +int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(m > n) { + DOMAIN_ERROR(result); + } + else if(m == n || m == 0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else { + gsl_sf_result nf; + gsl_sf_result mf; + gsl_sf_result nmmf; + if(m*2 > n) m = n-m; + gsl_sf_lnfact_e(n, &nf); + gsl_sf_lnfact_e(m, &mf); + gsl_sf_lnfact_e(n-m, &nmmf); + result->val = nf.val - mf.val - nmmf.val; + result->err = nf.err + mf.err + nmmf.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + + +int gsl_sf_choose_e(unsigned int n, unsigned int m, gsl_sf_result * result) +{ + if(m > n) { + DOMAIN_ERROR(result); + } + else if(m == n || m == 0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if (n <= GSL_SF_FACT_NMAX) { + result->val = (fact_table[n].f / fact_table[m].f) / fact_table[n-m].f; + result->err = 6.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } else { + if(m*2 < n) m = n-m; + + if (n - m < 64) /* compute product for a manageable number of terms */ + { + double prod = 1.0; + unsigned int k; + + for(k=n; k>=m+1; k--) { + double tk = (double)k / (double)(k-m); + if(tk > GSL_DBL_MAX/prod) { + OVERFLOW_ERROR(result); + } + prod *= tk; + } + result->val = prod; + result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs(n-m); + return GSL_SUCCESS; + } + else + { + gsl_sf_result lc; + const int stat_lc = gsl_sf_lnchoose_e (n, m, &lc); + const int stat_e = gsl_sf_exp_err_e(lc.val, lc.err, result); + return GSL_ERROR_SELECT_2(stat_lc, stat_e); + } + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_fact(const unsigned int n) +{ + EVAL_RESULT(gsl_sf_fact_e(n, &result)); +} + +double gsl_sf_lnfact(const unsigned int n) +{ + EVAL_RESULT(gsl_sf_lnfact_e(n, &result)); +} + +double gsl_sf_doublefact(const unsigned int n) +{ + EVAL_RESULT(gsl_sf_doublefact_e(n, &result)); +} + +double gsl_sf_lndoublefact(const unsigned int n) +{ + EVAL_RESULT(gsl_sf_lndoublefact_e(n, &result)); +} + +double gsl_sf_lngamma(const double x) +{ + EVAL_RESULT(gsl_sf_lngamma_e(x, &result)); +} + +double gsl_sf_gamma(const double x) +{ + EVAL_RESULT(gsl_sf_gamma_e(x, &result)); +} + +double gsl_sf_gammastar(const double x) +{ + EVAL_RESULT(gsl_sf_gammastar_e(x, &result)); +} + +double gsl_sf_gammainv(const double x) +{ + EVAL_RESULT(gsl_sf_gammainv_e(x, &result)); +} + +double gsl_sf_taylorcoeff(const int n, const double x) +{ + EVAL_RESULT(gsl_sf_taylorcoeff_e(n, x, &result)); +} + +double gsl_sf_choose(unsigned int n, unsigned int m) +{ + EVAL_RESULT(gsl_sf_choose_e(n, m, &result)); +} + +double gsl_sf_lnchoose(unsigned int n, unsigned int m) +{ + EVAL_RESULT(gsl_sf_lnchoose_e(n, m, &result)); +} |