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.

1 files changed, 1417 insertions, 0 deletions

diff --git a/gsl-1.9/specfunc/coulomb.c b/gsl-1.9/specfunc/coulomb.c
new file mode 100644
index 0000000..7c68076
--- /dev/null
+++ b/gsl-1.9/specfunc/coulomb.c
@@ -0,0 +1,1417 @@
+/* specfunc/coulomb.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 */
+
+/* Evaluation of Coulomb wave functions F_L(eta, x), G_L(eta, x),
+ * and their derivatives. A combination of Steed's method, asymptotic
+ * results, and power series.
+ *
+ * Steed's method:
+ *  [Barnett, CPC 21, 297 (1981)]
+ * Power series and other methods:
+ *  [Biedenharn et al., PR 97, 542 (1954)]
+ *  [Bardin et al., CPC 3, 73 (1972)]
+ *  [Abad+Sesma, CPC 71, 110 (1992)]
+ */
+#include <config.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_psi.h>
+#include <gsl/gsl_sf_airy.h>
+#include <gsl/gsl_sf_pow_int.h>
+#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_sf_coulomb.h>
+
+#include "error.h"
+
+/* the L=0 normalization constant
+ * [Abramowitz+Stegun 14.1.8]
+ */
+static
+double
+C0sq(double eta)
+{
+  double twopieta = 2.0*M_PI*eta;
+
+  if(fabs(eta) < GSL_DBL_EPSILON) {
+    return 1.0;
+  }
+  else if(twopieta > GSL_LOG_DBL_MAX) {
+    return 0.0;
+  }
+  else {
+    gsl_sf_result scale;
+    gsl_sf_expm1_e(twopieta, &scale);
+    return twopieta/scale.val;
+  }
+}
+
+
+/* the full definition of C_L(eta) for any valid L and eta
+ * [Abramowitz and Stegun 14.1.7]
+ * This depends on the complex gamma function. For large
+ * arguments the phase of the complex gamma function is not
+ * very accurately determined. However the modulus is, and that
+ * is all that we need to calculate C_L.
+ *
+ * This is not valid for L <= -3/2  or  L = -1.
+ */
+static
+int
+CLeta(double L, double eta, gsl_sf_result * result)
+{
+  gsl_sf_result ln1; /* log of numerator Gamma function */
+  gsl_sf_result ln2; /* log of denominator Gamma function */
+  double sgn = 1.0;
+  double arg_val, arg_err;
+
+  if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
+    gsl_sf_lngamma_e(L+1.0, &ln1);
+  }
+  else {
+    gsl_sf_result p1;                 /* phase of numerator Gamma -- not used */
+    gsl_sf_lngamma_complex_e(L+1.0, eta, &ln1, &p1); /* should be ok */
+  }
+
+  gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
+  if(L < -1.0) sgn = -sgn;
+
+  arg_val  = L*M_LN2 - 0.5*eta*M_PI + ln1.val - ln2.val;
+  arg_err  = ln1.err + ln2.err;
+  arg_err += GSL_DBL_EPSILON * (fabs(L*M_LN2) + fabs(0.5*eta*M_PI));
+  return gsl_sf_exp_err_e(arg_val, arg_err, result);
+}
+
+
+int
+gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
+{
+  /* CHECK_POINTER(result) */
+
+  if(lam <= -1.0) {
+    DOMAIN_ERROR(result);
+  }
+  else if(fabs(lam) < GSL_DBL_EPSILON) {
+    /* saves a calculation of complex_lngamma(), otherwise not necessary */
+    result->val = sqrt(C0sq(eta));
+    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
+    return GSL_SUCCESS;
+  }
+  else {
+    return CLeta(lam, eta, result);
+  }
+}
+
+
+/* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
+ */
+int
+gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
+{
+  int k;
+  gsl_sf_result cl_0;
+  gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
+  cl[0] = cl_0.val;
+
+  for(k=1; k<=kmax; k++) {
+    double L = lam_min + k;
+    cl[k] = cl[k-1] * hypot(L, eta)/(L*(2.0*L+1.0));
+  }
+
+  return GSL_SUCCESS;
+}
+
+
+/* Evaluate the series for Phi_L(eta,x) and Phi_L*(eta,x)
+ * [Abramowitz+Stegun 14.1.5]
+ * [Abramowitz+Stegun 14.1.13]
+ *
+ * The sequence of coefficients A_k^L is
+ * manifestly well-controlled for L >= -1/2
+ * and eta < 10.
+ *
+ * This makes sense since this is the region
+ * away from threshold, and you expect
+ * the evaluation to become easier as you
+ * get farther from threshold.
+ *
+ * Empirically, this is quite well-behaved for
+ *   L >= -1/2
+ *   eta < 10
+ *   x   < 10
+ */
+#if 0
+static
+int
+coulomb_Phi_series(const double lam, const double eta, const double x,
+                   double * result, double * result_star)
+{
+  int kmin =   5;
+  int kmax = 200;
+  int k;
+  double Akm2 = 1.0;
+  double Akm1 = eta/(lam+1.0);
+  double Ak;
+
+  double xpow = x;
+  double sum  = Akm2 + Akm1*x;
+  double sump = (lam+1.0)*Akm2 + (lam+2.0)*Akm1*x;
+  double prev_abs_del   = fabs(Akm1*x);
+  double prev_abs_del_p = (lam+2.0) * prev_abs_del;
+
+  for(k=2; k<kmax; k++) {
+    double del;
+    double del_p;
+    double abs_del;
+    double abs_del_p;
+
+    Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));
+
+    xpow *= x;
+    del   = Ak*xpow;
+    del_p = (k+lam+1.0)*del;
+    sum  += del;
+    sump += del_p;
+
+    abs_del   = fabs(del);
+    abs_del_p = fabs(del_p);
+
+    if(          abs_del/(fabs(sum)+abs_del)          < GSL_DBL_EPSILON
+       &&   prev_abs_del/(fabs(sum)+prev_abs_del)     < GSL_DBL_EPSILON
+       &&      abs_del_p/(fabs(sump)+abs_del_p)       < GSL_DBL_EPSILON
+       && prev_abs_del_p/(fabs(sump)+prev_abs_del_p)  < GSL_DBL_EPSILON
+       && k > kmin
+       ) break;
+
+    /* We need to keep track of the previous delta because when
+     * eta is near zero the odd terms of the sum are very small
+     * and this could lead to premature termination.
+     */
+    prev_abs_del   = abs_del;
+    prev_abs_del_p = abs_del_p;
+
+    Akm2 = Akm1;
+    Akm1 = Ak;
+  }
+
+  *result      = sum;
+  *result_star = sump;
+
+  if(k==kmax) {
+    GSL_ERROR ("error", GSL_EMAXITER);
+  }
+  else {
+    return GSL_SUCCESS;
+  }
+}
+#endif /* 0 */
+
+
+/* Determine the connection phase, phi_lambda.
+ * See coulomb_FG_series() below. We have
+ * to be careful about sin(phi)->0. Note that
+ * there is an underflow condition for large 
+ * positive eta in any case.
+ */
+static
+int
+coulomb_connection(const double lam, const double eta,
+                   double * cos_phi, double * sin_phi)
+{
+  if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
+    *cos_phi = 1.0;
+    *sin_phi = 0.0;
+    GSL_ERROR ("error", GSL_EUNDRFLW);
+  }
+  else if(eta > -GSL_LOG_DBL_EPSILON/(4.0*M_PI)) {
+    const double eps = 2.0 * exp(-2.0*M_PI*eta);
+    const double tpl = tan(M_PI * lam);
+    const double dth = eps * tpl / (tpl*tpl + 1.0);
+    *cos_phi = -1.0 + 0.5 * dth*dth;
+    *sin_phi = -dth;
+    return GSL_SUCCESS;
+  }
+  else {
+    double X   = tanh(M_PI * eta) / tan(M_PI * lam);
+    double phi = -atan(X) - (lam + 0.5) * M_PI;
+    *cos_phi = cos(phi);
+    *sin_phi = sin(phi);
+    return GSL_SUCCESS;
+  }
+}
+
+
+/* Evaluate the Frobenius series for F_lam(eta,x) and G_lam(eta,x).
+ * Homegrown algebra. Evaluates the series for F_{lam} and
+ * F_{-lam-1}, then uses
+ *    G_{lam} = (F_{lam} cos(phi) - F_{-lam-1}) / sin(phi)
+ * where
+ *    phi = Arg[Gamma[1+lam+I eta]] - Arg[Gamma[-lam + I eta]] - (lam+1/2)Pi
+ *        = Arg[Sin[Pi(-lam+I eta)] - (lam+1/2)Pi
+ *        = atan2(-cos(lam Pi)sinh(eta Pi), -sin(lam Pi)cosh(eta Pi)) - (lam+1/2)Pi
+ *
+ *        = -atan(X) - (lam+1/2) Pi,  X = tanh(eta Pi)/tan(lam Pi)
+ *
+ * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
+ */
+static
+int
+coulomb_FG_series(const double lam, const double eta, const double x,
+                  gsl_sf_result * F, gsl_sf_result * G)
+{
+  const int max_iter = 800;
+  gsl_sf_result ClamA;
+  gsl_sf_result ClamB;
+  int stat_A = CLeta(lam, eta, &ClamA);
+  int stat_B = CLeta(-lam-1.0, eta, &ClamB);
+  const double tlp1 = 2.0*lam + 1.0;
+  const double pow_x = pow(x, lam);
+  double cos_phi_lam;
+  double sin_phi_lam;
+
+  double uA_mm2 = 1.0;                  /* uA sum is for F_{lam} */
+  double uA_mm1 = x*eta/(lam+1.0);
+  double uA_m;
+  double uB_mm2 = 1.0;                  /* uB sum is for F_{-lam-1} */
+  double uB_mm1 = -x*eta/lam;
+  double uB_m;
+  double A_sum = uA_mm2 + uA_mm1;
+  double B_sum = uB_mm2 + uB_mm1;
+  double A_abs_del_prev = fabs(A_sum);
+  double B_abs_del_prev = fabs(B_sum);
+  gsl_sf_result FA, FB;
+  int m = 2;
+
+  int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);
+
+  if(stat_conn == GSL_EUNDRFLW) {
+    F->val = 0.0;  /* FIXME: should this be set to Inf too like G? */
+    F->err = 0.0;
+    OVERFLOW_ERROR(G);
+  }
+
+  while(m < max_iter) {
+    double abs_dA;
+    double abs_dB;
+    uA_m = x*(2.0*eta*uA_mm1 - x*uA_mm2)/(m*(m+tlp1));
+    uB_m = x*(2.0*eta*uB_mm1 - x*uB_mm2)/(m*(m-tlp1));
+    A_sum += uA_m;
+    B_sum += uB_m;
+    abs_dA = fabs(uA_m);
+    abs_dB = fabs(uB_m);
+    if(m > 15) {
+      /* Don't bother checking until we have gone out a little ways;
+       * a minor optimization. Also make sure to check both the
+       * current and the previous increment because the odd and even
+       * terms of the sum can have very different behaviour, depending
+       * on the value of eta.
+       */
+      double max_abs_dA = GSL_MAX(abs_dA, A_abs_del_prev);
+      double max_abs_dB = GSL_MAX(abs_dB, B_abs_del_prev);
+      double abs_A = fabs(A_sum);
+      double abs_B = fabs(B_sum);
+      if(   max_abs_dA/(max_abs_dA + abs_A) < 4.0*GSL_DBL_EPSILON
+         && max_abs_dB/(max_abs_dB + abs_B) < 4.0*GSL_DBL_EPSILON
+         ) break;
+    }
+    A_abs_del_prev = abs_dA;
+    B_abs_del_prev = abs_dB;
+    uA_mm2 = uA_mm1;
+    uA_mm1 = uA_m;
+    uB_mm2 = uB_mm1;
+    uB_mm1 = uB_m;
+    m++;
+  }
+
+  FA.val = A_sum * ClamA.val * pow_x * x;
+  FA.err = fabs(A_sum) * ClamA.err * pow_x * x + 2.0*GSL_DBL_EPSILON*fabs(FA.val);
+  FB.val = B_sum * ClamB.val / pow_x;
+  FB.err = fabs(B_sum) * ClamB.err / pow_x + 2.0*GSL_DBL_EPSILON*fabs(FB.val);
+
+  F->val = FA.val;
+  F->err = FA.err;
+
+  G->val = (FA.val * cos_phi_lam - FB.val)/sin_phi_lam;
+  G->err = (FA.err * fabs(cos_phi_lam) + FB.err)/fabs(sin_phi_lam);
+
+  if(m >= max_iter)
+    GSL_ERROR ("error", GSL_EMAXITER);
+  else
+    return GSL_ERROR_SELECT_2(stat_A, stat_B);
+}
+
+
+/* Evaluate the Frobenius series for F_0(eta,x) and G_0(eta,x).
+ * See [Bardin et al., CPC 3, 73 (1972), (14)-(17)];
+ * note the misprint in (17): nu_0=1 is correct, not nu_0=0.
+ */
+static
+int
+coulomb_FG0_series(const double eta, const double x,
+                   gsl_sf_result * F, gsl_sf_result * G)
+{
+  const int max_iter = 800;
+  const double x2  = x*x;
+  const double tex = 2.0*eta*x;
+  gsl_sf_result C0;
+  int stat_CL = CLeta(0.0, eta, &C0);
+  gsl_sf_result r1pie;
+  int psi_stat = gsl_sf_psi_1piy_e(eta, &r1pie);
+  double u_mm2 = 0.0;  /* u_0 */
+  double u_mm1 = x;    /* u_1 */
+  double u_m;
+  double v_mm2 = 1.0;                               /* nu_0 */
+  double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val);   /* nu_1 */
+  double v_m;
+  double u_sum = u_mm2 + u_mm1;
+  double v_sum = v_mm2 + v_mm1;
+  double u_abs_del_prev = fabs(u_sum);
+  double v_abs_del_prev = fabs(v_sum);
+  int m = 2;
+  double u_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(u_sum);
+  double v_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(v_sum);
+  double ln2x = log(2.0*x);
+
+  while(m < max_iter) {
+    double abs_du;
+    double abs_dv;
+    double m_mm1 = m*(m-1.0);
+    u_m = (tex*u_mm1 - x2*u_mm2)/m_mm1;
+    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*eta*(2*m-1)*u_m)/m_mm1;
+    u_sum += u_m;
+    v_sum += v_m;
+    abs_du = fabs(u_m);
+    abs_dv = fabs(v_m);
+    u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
+    v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
+    if(m > 15) {
+      /* Don't bother checking until we have gone out a little ways;
+       * a minor optimization. Also make sure to check both the
+       * current and the previous increment because the odd and even
+       * terms of the sum can have very different behaviour, depending
+       * on the value of eta.
+       */
+      double max_abs_du = GSL_MAX(abs_du, u_abs_del_prev);
+      double max_abs_dv = GSL_MAX(abs_dv, v_abs_del_prev);
+      double abs_u = fabs(u_sum);
+      double abs_v = fabs(v_sum);
+      if(   max_abs_du/(max_abs_du + abs_u) < 40.0*GSL_DBL_EPSILON
+         && max_abs_dv/(max_abs_dv + abs_v) < 40.0*GSL_DBL_EPSILON
+         ) break;
+    }
+    u_abs_del_prev = abs_du;
+    v_abs_del_prev = abs_dv;
+    u_mm2 = u_mm1;
+    u_mm1 = u_m;
+    v_mm2 = v_mm1;
+    v_mm1 = v_m;
+    m++;
+  }
+
+  F->val  = C0.val * u_sum;
+  F->err  = C0.err * fabs(u_sum);
+  F->err += fabs(C0.val) * u_sum_err;
+  F->err += 2.0 * GSL_DBL_EPSILON * fabs(F->val);
+
+  G->val  = (v_sum + 2.0*eta*u_sum * ln2x) / C0.val;
+  G->err  = (fabs(v_sum) + fabs(2.0*eta*u_sum * ln2x)) / fabs(C0.val) * fabs(C0.err/C0.val);
+  G->err += (v_sum_err + fabs(2.0*eta*u_sum_err*ln2x)) / fabs(C0.val);
+  G->err += 2.0 * GSL_DBL_EPSILON * fabs(G->val);
+
+  if(m == max_iter)
+    GSL_ERROR ("error", GSL_EMAXITER);
+  else
+    return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
+}
+
+
+/* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
+ * Homegrown algebra.
+ */
+static
+int
+coulomb_FGmhalf_series(const double eta, const double x,
+                       gsl_sf_result * F, gsl_sf_result * G)
+{
+  const int max_iter = 800;
+  const double rx  = sqrt(x);
+  const double x2  = x*x;
+  const double tex = 2.0*eta*x;
+  gsl_sf_result Cmhalf;
+  int stat_CL = CLeta(-0.5, eta, &Cmhalf);
+  double u_mm2 = 1.0;                      /* u_0 */
+  double u_mm1 = tex * u_mm2;              /* u_1 */
+  double u_m;
+  double v_mm2, v_mm1, v_m;
+  double f_sum, g_sum;
+  double tmp1;
+  gsl_sf_result rpsi_1pe;
+  gsl_sf_result rpsi_1p2e;
+  int m = 2;
+
+  gsl_sf_psi_1piy_e(eta,     &rpsi_1pe);
+  gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);
+
+  v_mm2 = 2.0*M_EULER - M_LN2 - rpsi_1pe.val + 2.0*rpsi_1p2e.val;
+  v_mm1 = tex*(v_mm2 - 2.0*u_mm2);
+
+  f_sum = u_mm2 + u_mm1;
+  g_sum = v_mm2 + v_mm1;
+
+  while(m < max_iter) {
+    double m2 = m*m;
+    u_m = (tex*u_mm1 - x2*u_mm2)/m2;
+    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*m*u_m)/m2;
+    f_sum += u_m;
+    g_sum += v_m;
+    if(   f_sum != 0.0
+       && g_sum != 0.0
+       && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;
+    u_mm2 = u_mm1;
+    u_mm1 = u_m;
+    v_mm2 = v_mm1;
+    v_mm1 = v_m;
+    m++;
+  }
+  
+  F->val = Cmhalf.val * rx * f_sum;
+  F->err = Cmhalf.err * fabs(rx * f_sum) + 2.0*GSL_DBL_EPSILON*fabs(F->val);
+
+  tmp1 = f_sum*log(x);
+  G->val = -rx*(tmp1 + g_sum)/Cmhalf.val;
+  G->err = fabs(rx)*(fabs(tmp1) + fabs(g_sum))/fabs(Cmhalf.val) * fabs(Cmhalf.err/Cmhalf.val);
+
+  if(m == max_iter)
+    GSL_ERROR ("error", GSL_EMAXITER);
+  else
+    return stat_CL;
+}
+
+
+/* Evolve the backwards recurrence for F,F'.
+ *
+ *    F_{lam-1}  = (S_lam F_lam + F_lam') / R_lam
+ *    F_{lam-1}' = (S_lam F_{lam-1} - R_lam F_lam)
+ * where
+ *    R_lam = sqrt(1 + (eta/lam)^2)
+ *    S_lam = lam/x + eta/lam
+ *
+ */
+static
+int
+coulomb_F_recur(double lam_min, int kmax,
+                double eta, double x,
+                double F_lam_max, double Fp_lam_max,
+                double * F_lam_min, double * Fp_lam_min
+                )
+{
+  double x_inv = 1.0/x;
+  double fcl = F_lam_max;
+  double fpl = Fp_lam_max;
+  double lam_max = lam_min + kmax;
+  double lam = lam_max;
+  int k;
+
+  for(k=kmax-1; k>=0; k--) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el  + lam*x_inv;
+    double fc_lm1;
+    fc_lm1 = (fcl*sl + fpl)/rl;
+    fpl    =  fc_lm1*sl - fcl*rl;
+    fcl    =  fc_lm1;
+    lam -= 1.0;
+  }
+
+  *F_lam_min  = fcl;
+  *Fp_lam_min = fpl;  
+  return GSL_SUCCESS;
+}
+
+
+/* Evolve the forward recurrence for G,G'.
+ *
+ *   G_{lam+1}  = (S_lam G_lam - G_lam')/R_lam
+ *   G_{lam+1}' = R_{lam+1} G_lam - S_lam G_{lam+1}
+ *
+ * where S_lam and R_lam are as above in the F recursion.
+ */
+static
+int
+coulomb_G_recur(const double lam_min, const int kmax,
+                const double eta, const double x,
+                const double G_lam_min, const double Gp_lam_min,
+                double * G_lam_max, double * Gp_lam_max
+                )
+{
+  double x_inv = 1.0/x;
+  double gcl = G_lam_min;
+  double gpl = Gp_lam_min;
+  double lam = lam_min + 1.0;
+  int k;
+
+  for(k=1; k<=kmax; k++) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el + lam*x_inv;
+    double gcl1 = (sl*gcl - gpl)/rl;
+    gpl   = rl*gcl - sl*gcl1;
+    gcl   = gcl1;
+    lam += 1.0;
+  }
+  
+  *G_lam_max  = gcl;
+  *Gp_lam_max = gpl;
+  return GSL_SUCCESS;
+}
+
+
+/* Evaluate the first continued fraction, giving
+ * the ratio F'/F at the upper lambda value.
+ * We also determine the sign of F at that point,
+ * since it is the sign of the last denominator
+ * in the continued fraction.
+ */
+static
+int
+coulomb_CF1(double lambda,
+            double eta, double x,
+            double * fcl_sign,
+            double * result,
+            int * count
+            )
+{
+  const double CF1_small = 1.e-30;
+  const double CF1_abort = 1.0e+05;
+  const double CF1_acc   = 2.0*GSL_DBL_EPSILON;
+  const double x_inv     = 1.0/x;
+  const double px        = lambda + 1.0 + CF1_abort;
+
+  double pk = lambda + 1.0;
+  double F  = eta/pk + pk*x_inv;
+  double D, C;
+  double df;
+
+  *fcl_sign = 1.0;
+  *count = 0;
+
+  if(fabs(F) < CF1_small) F = CF1_small;
+  D = 0.0;
+  C = F;
+
+  do {
+    double pk1 = pk + 1.0;
+    double ek  = eta / pk;
+    double rk2 = 1.0 + ek*ek;
+    double tk  = (pk + pk1)*(x_inv + ek/pk1);
+    D   =  tk - rk2 * D;
+    C   =  tk - rk2 / C;
+    if(fabs(C) < CF1_small) C = CF1_small;
+    if(fabs(D) < CF1_small) D = CF1_small;
+    D = 1.0/D;
+    df = D * C;
+    F  = F * df;
+    if(D < 0.0) {
+      /* sign of result depends on sign of denominator */
+      *fcl_sign = - *fcl_sign;
+    }
+    pk = pk1;
+    if( pk > px ) {
+      *result = F;
+      GSL_ERROR ("error", GSL_ERUNAWAY);
+    }
+    ++(*count);
+  }
+  while(fabs(df-1.0) > CF1_acc);
+  
+  *result = F;
+  return GSL_SUCCESS;
+}
+
+
+#if 0
+static
+int
+old_coulomb_CF1(const double lambda,
+                double eta, double x,
+                double * fcl_sign,
+                double * result
+                )
+{
+  const double CF1_abort = 1.e5;
+  const double CF1_acc   = 10.0*GSL_DBL_EPSILON;
+  const double x_inv     = 1.0/x;
+  const double px        = lambda + 1.0 + CF1_abort;
+  
+  double pk = lambda + 1.0;
+  
+  double D;
+  double df;
+
+  double F;
+  double p;
+  double pk1;
+  double ek;
+  
+  double fcl = 1.0;
+
+  double tk;
+
+  while(1) {
+    ek = eta/pk;
+    F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;
+    pk1 = pk + 1.0;
+    if(fabs(eta*x + pk*pk1) > CF1_acc) break;
+    fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));
+    pk = 2.0 + pk;
+  }
+
+  D  = 1.0/((pk + pk1)*(x_inv + ek/pk1));
+  df = -fcl*(1.0 + ek*ek)*D;
+  
+  if(fcl != 1.0) fcl = -1.0;
+  if(D    < 0.0) fcl = -fcl;
+  
+  F = F + df;
+
+  p = 1.0;
+  do {
+    pk = pk1;
+    pk1 = pk + 1.0;
+    ek  = eta / pk;
+    tk  = (pk + pk1)*(x_inv + ek/pk1);
+    D   =  tk - D*(1.0+ek*ek);
+    if(fabs(D) < sqrt(CF1_acc)) {
+      p += 1.0;
+      if(p > 2.0) {
+        printf("HELP............\n");
+      }
+    }
+    D = 1.0/D;
+    if(D < 0.0) {
+      /* sign of result depends on sign of denominator */
+      fcl = -fcl;
+    }
+    df = df*(D*tk - 1.0);
+    F  = F + df;
+    if( pk > px ) {
+      GSL_ERROR ("error", GSL_ERUNAWAY);
+    }
+  }
+  while(fabs(df) > fabs(F)*CF1_acc);
+  
+  *fcl_sign = fcl;
+  *result = F;
+  return GSL_SUCCESS;
+}
+#endif /* 0 */
+
+
+/* Evaluate the second continued fraction to 
+ * obtain the ratio
+ *    (G' + i F')/(G + i F) := P + i Q
+ * at the specified lambda value.
+ */
+static
+int
+coulomb_CF2(const double lambda, const double eta, const double x,
+            double * result_P, double * result_Q, int * count
+            )
+{
+  int status = GSL_SUCCESS;
+
+  const double CF2_acc   = 4.0*GSL_DBL_EPSILON;
+  const double CF2_abort = 2.0e+05;
+
+  const double wi    = 2.0*eta;
+  const double x_inv = 1.0/x;
+  const double e2mm1 = eta*eta + lambda*(lambda + 1.0);
+  
+  double ar = -e2mm1;
+  double ai =  eta;
+
+  double br =  2.0*(x - eta);
+  double bi =  2.0;
+
+  double dr =  br/(br*br + bi*bi);
+  double di = -bi/(br*br + bi*bi);
+
+  double dp = -x_inv*(ar*di + ai*dr);
+  double dq =  x_inv*(ar*dr - ai*di);
+
+  double A, B, C, D;
+
+  double pk =  0.0;
+  double P  =  0.0;
+  double Q  =  1.0 - eta*x_inv;
+
+  *count = 0;
+ 
+  do {
+    P += dp;
+    Q += dq;
+    pk += 2.0;
+    ar += pk;
+    ai += wi;
+    bi += 2.0;
+    D  = ar*dr - ai*di + br;
+    di = ai*dr + ar*di + bi;
+    C  = 1.0/(D*D + di*di);
+    dr =  C*D;
+    di = -C*di;
+    A  = br*dr - bi*di - 1.;
+    B  = bi*dr + br*di;
+    C  = dp*A  - dq*B;
+    dq = dp*B  + dq*A;
+    dp = C;
+    if(pk > CF2_abort) {
+      status = GSL_ERUNAWAY;
+      break;
+    }
+    ++(*count);
+  }
+  while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);
+
+  if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {
+    status = GSL_ELOSS;
+  }
+
+  *result_P = P;
+  *result_Q = Q;
+  return status;
+}
+
+
+/* WKB evaluation of F, G. Assumes  0 < x < turning point.
+ * Overflows are trapped, GSL_EOVRFLW is signalled,
+ * and an exponent is returned such that:
+ *
+ *   result_F = fjwkb * exp(-exponent)
+ *   result_G = gjwkb * exp( exponent)
+ *
+ * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]
+ *
+ * Unfortunately, this is not very accurate in general. The
+ * test cases typically have 3-4 digits of precision. One could
+ * argue that this is ok for general use because, for instance,
+ * F is exponentially small in this region and so the absolute
+ * accuracy is still roughly acceptable. But it would be better
+ * to have a systematic method for improving the precision. See
+ * the Abad+Sesma method discussion below.
+ */
+static
+int
+coulomb_jwkb(const double lam, const double eta, const double x,
+             gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,
+             double * exponent)
+{
+  const double llp1      = lam*(lam+1.0) + 6.0/35.0;
+  const double llp1_eff  = GSL_MAX(llp1, 0.0);
+  const double rho_ghalf = sqrt(x*(2.0*eta - x) + llp1_eff);
+  const double sinh_arg  = sqrt(llp1_eff/(eta*eta+llp1_eff)) * rho_ghalf / x;
+  const double sinh_inv  = log(sinh_arg + hypot(1.0,sinh_arg));
+
+  const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);
+
+  const double zeta_half = pow(3.0*phi/2.0, 1.0/3.0);
+  const double prefactor = sqrt(M_PI*phi*x/(6.0 * rho_ghalf));
+  
+  double F = prefactor * 3.0/zeta_half;
+  double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */
+  double F_exp;
+  double G_exp;
+  
+  const double airy_scale_exp = phi;
+  gsl_sf_result ai;
+  gsl_sf_result bi;
+  gsl_sf_airy_Ai_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &ai);
+  gsl_sf_airy_Bi_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &bi);
+  F *= ai.val;
+  G *= bi.val;
+  F_exp = log(F) - airy_scale_exp;
+  G_exp = log(G) + airy_scale_exp;
+
+  if(G_exp >= GSL_LOG_DBL_MAX) {
+    fjwkb->val = F;
+    gjwkb->val = G;
+    fjwkb->err = 1.0e-3 * fabs(F); /* FIXME: real error here ... could be smaller */
+    gjwkb->err = 1.0e-3 * fabs(G);
+    *exponent = airy_scale_exp;
+    GSL_ERROR ("error", GSL_EOVRFLW);
+  }
+  else {
+    fjwkb->val = exp(F_exp);
+    gjwkb->val = exp(G_exp);
+    fjwkb->err = 1.0e-3 * fabs(fjwkb->val);
+    gjwkb->err = 1.0e-3 * fabs(gjwkb->val);
+    *exponent = 0.0;
+    return GSL_SUCCESS;
+  }
+}
+
+
+/* Asymptotic evaluation of F and G below the minimal turning point.
+ *
+ * This is meant to be a drop-in replacement for coulomb_jwkb().
+ * It uses the expressions in [Abad+Sesma]. This requires some
+ * work because I am not sure where it is valid. They mumble
+ * something about |x| < |lam|^(-1/2) or 8|eta x| > lam when |x| < 1.
+ * This seems true, but I thought the result was based on a uniform
+ * expansion and could be controlled by simply using more terms.
+ */
+#if 0
+static
+int
+coulomb_AS_xlt2eta(const double lam, const double eta, const double x,
+                   gsl_sf_result * f_AS, gsl_sf_result * g_AS,
+                   double * exponent)
+{
+  /* no time to do this now... */
+}
+#endif /* 0 */
+
+
+
+/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
+
+int
+gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
+                            const double lam_F,
+                            const int  k_lam_G,      /* lam_G = lam_F - k_lam_G */
+                            gsl_sf_result * F, gsl_sf_result * Fp,
+                            gsl_sf_result * G, gsl_sf_result * Gp,
+                            double * exp_F, double * exp_G)
+{
+  const double lam_G = lam_F - k_lam_G;
+
+  if(x < 0.0 || lam_F <= -0.5 || lam_G <= -0.5) {
+    GSL_SF_RESULT_SET(F,  0.0, 0.0);
+    GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
+    GSL_SF_RESULT_SET(G,  0.0, 0.0);
+    GSL_SF_RESULT_SET(Gp, 0.0, 0.0);
+    *exp_F = 0.0;
+    *exp_G = 0.0;
+    GSL_ERROR ("domain error", GSL_EDOM);
+  }
+  else if(x == 0.0) {
+    gsl_sf_result C0;
+    CLeta(0.0, eta, &C0);
+    GSL_SF_RESULT_SET(F,  0.0, 0.0);
+    GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
+    GSL_SF_RESULT_SET(G,  0.0, 0.0); /* FIXME: should be Inf */
+    GSL_SF_RESULT_SET(Gp, 0.0, 0.0); /* FIXME: should be Inf */
+    *exp_F = 0.0;
+    *exp_G = 0.0;
+    if(lam_F == 0.0){
+      GSL_SF_RESULT_SET(Fp, C0.val, C0.err);
+    }
+    if(lam_G == 0.0) {
+      GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val));
+    }
+    GSL_ERROR ("domain error", GSL_EDOM);
+    /* After all, since we are asking for G, this is a domain error... */
+  }
+  else if(x < 1.2 && 2.0*M_PI*eta < 0.9*(-GSL_LOG_DBL_MIN) && fabs(eta*x) < 10.0) {
+    /* Reduce to a small lambda value and use the series
+     * representations for F and G. We cannot allow eta to
+     * be large and positive because the connection formula
+     * for G_lam is badly behaved due to an underflow in sin(phi_lam) 
+     * [see coulomb_FG_series() and coulomb_connection() above].
+     * Note that large negative eta is ok however.
+     */
+    const double SMALL = GSL_SQRT_DBL_EPSILON;
+    const int N    = (int)(lam_F + 0.5);
+    const int span = GSL_MAX(k_lam_G, N);
+    const double lam_min = lam_F - N;    /* -1/2 <= lam_min < 1/2 */
+    double F_lam_F, Fp_lam_F;
+    double G_lam_G, Gp_lam_G;
+    double F_lam_F_err, Fp_lam_F_err;
+    double Fp_over_F_lam_F;
+    double F_sign_lam_F;
+    double F_lam_min_unnorm, Fp_lam_min_unnorm;
+    double Fp_over_F_lam_min;
+    gsl_sf_result F_lam_min;
+    gsl_sf_result G_lam_min, Gp_lam_min;
+    double F_scale;
+    double Gerr_frac;
+    double F_scale_frac_err;
+    double F_unnorm_frac_err;
+
+    /* Determine F'/F at lam_F. */
+    int CF1_count;
+    int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
+
+    int stat_ser;
+    int stat_Fr;
+    int stat_Gr;
+
+    /* Recurse down with unnormalized F,F' values. */
+    F_lam_F  = SMALL;
+    Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
+    if(span != 0) {
+      stat_Fr = coulomb_F_recur(lam_min, span, eta, x,
+                                F_lam_F, Fp_lam_F,
+                                &F_lam_min_unnorm, &Fp_lam_min_unnorm
+                                );
+    }
+    else {
+      F_lam_min_unnorm  =  F_lam_F;
+      Fp_lam_min_unnorm = Fp_lam_F;
+      stat_Fr = GSL_SUCCESS;
+    }
+
+    /* Determine F and G at lam_min. */
+    if(lam_min == -0.5) {
+      stat_ser = coulomb_FGmhalf_series(eta, x, &F_lam_min, &G_lam_min);
+    }
+    else if(lam_min == 0.0) {
+      stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min);
+    }
+    else if(lam_min == 0.5) {
+      /* This cannot happen. */
+      F->val  = F_lam_F;
+      F->err  = 2.0 * GSL_DBL_EPSILON * fabs(F->val);
+      Fp->val = Fp_lam_F;
+      Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val);
+      G->val  = G_lam_G;
+      G->err  = 2.0 * GSL_DBL_EPSILON * fabs(G->val);
+      Gp->val = Gp_lam_G;
+      Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val);
+      *exp_F = 0.0;
+      *exp_G = 0.0;
+      GSL_ERROR ("error", GSL_ESANITY);
+    }
+    else {
+      stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min);
+    }
+
+    /* Determine remaining quantities. */
+    Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
+    Gp_lam_min.val  = Fp_over_F_lam_min*G_lam_min.val - 1.0/F_lam_min.val;
+    Gp_lam_min.err  = fabs(Fp_over_F_lam_min)*G_lam_min.err;
+    Gp_lam_min.err += fabs(1.0/F_lam_min.val) * fabs(F_lam_min.err/F_lam_min.val);
+    F_scale     = F_lam_min.val / F_lam_min_unnorm;
+
+    /* Apply scale to the original F,F' values. */
+    F_scale_frac_err  = fabs(F_lam_min.err/F_lam_min.val);
+    F_unnorm_frac_err = 2.0*GSL_DBL_EPSILON*(CF1_count+span+1);
+    F_lam_F     *= F_scale;
+    F_lam_F_err  = fabs(F_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
+    Fp_lam_F    *= F_scale;
+    Fp_lam_F_err = fabs(Fp_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
+
+    /* Recurse up to get the required G,G' values. */
+    stat_Gr = coulomb_G_recur(lam_min, GSL_MAX(N-k_lam_G,0), eta, x,
+                              G_lam_min.val, Gp_lam_min.val,
+                              &G_lam_G, &Gp_lam_G
+                              );
+
+    F->val  = F_lam_F;
+    F->err  = F_lam_F_err;
+    F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F);
+
+    Fp->val  = Fp_lam_F;
+    Fp->err  = Fp_lam_F_err;
+    Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F);
+
+    Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val);
+
+    G->val  = G_lam_G;
+    G->err  = Gerr_frac * fabs(G_lam_G);
+    G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val);
+
+    Gp->val  = Gp_lam_G;
+    Gp->err  = Gerr_frac * fabs(Gp->val);
+    Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val);
+
+    *exp_F = 0.0;
+    *exp_G = 0.0;
+
+    return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr);
+  }
+  else if(x < 2.0*eta) {
+    /* Use WKB approximation to obtain F and G at the two
+     * lambda values, and use the Wronskian and the
+     * continued fractions for F'/F to obtain F' and G'.
+     */
+    gsl_sf_result F_lam_F, G_lam_F;
+    gsl_sf_result F_lam_G, G_lam_G;
+    double exp_lam_F, exp_lam_G;
+    int stat_lam_F;
+    int stat_lam_G;
+    int stat_CF1_lam_F;
+    int stat_CF1_lam_G;
+    int CF1_count;
+    double Fp_over_F_lam_F;
+    double Fp_over_F_lam_G;
+    double F_sign_lam_F;
+    double F_sign_lam_G;
+
+    stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F);
+    if(k_lam_G == 0) {
+      stat_lam_G = stat_lam_F;
+      F_lam_G = F_lam_F;
+      G_lam_G = G_lam_F;
+      exp_lam_G = exp_lam_F;
+    }
+    else {
+      stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G);
+    }
+
+    stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
+    if(k_lam_G == 0) {
+      stat_CF1_lam_G  = stat_CF1_lam_F;
+      F_sign_lam_G    = F_sign_lam_F;
+      Fp_over_F_lam_G = Fp_over_F_lam_F;
+    }
+    else {
+      stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count);
+    }
+
+    F->val = F_lam_F.val;
+    F->err = F_lam_F.err;
+
+    G->val = G_lam_G.val;
+    G->err = G_lam_G.err;
+
+    Fp->val  = Fp_over_F_lam_F * F_lam_F.val;
+    Fp->err  = fabs(Fp_over_F_lam_F) * F_lam_F.err;
+    Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val);
+
+    Gp->val  = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val;
+    Gp->err  = fabs(Fp_over_F_lam_G) * G_lam_G.err;
+    Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val);
+
+    *exp_F = exp_lam_F;
+    *exp_G = exp_lam_G;
+
+    if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) {
+      GSL_ERROR ("overflow", GSL_EOVRFLW);
+    }
+    else {
+      return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G);
+    }
+  }
+  else {
+    /* x > 2 eta, so we know that we can find a lambda value such
+     * that x is above the turning point. We do this, evaluate
+     * using Steed's method at that oscillatory point, then
+     * use recursion on F and G to obtain the required values.
+     *
+     * lam_0   = a value of lambda such that x is below the turning point
+     * lam_min = minimum of lam_0 and the requested lam_G, since
+     *           we must go at least as low as lam_G
+     */
+    const double SMALL = GSL_SQRT_DBL_EPSILON;
+    const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta));
+    const int N = ceil(lam_F - C + 0.5);
+    const double lam_0   = lam_F - GSL_MAX(N, 0);
+    const double lam_min = GSL_MIN(lam_0, lam_G);
+    double F_lam_F, Fp_lam_F;
+    double G_lam_G, Gp_lam_G;
+    double F_lam_min_unnorm, Fp_lam_min_unnorm;
+    double F_lam_min, Fp_lam_min;
+    double G_lam_min, Gp_lam_min;
+    double Fp_over_F_lam_F;
+    double Fp_over_F_lam_min;
+    double F_sign_lam_F, F_sign_lam_min;
+    double P_lam_min, Q_lam_min;
+    double alpha;
+    double gamma;
+    double F_scale;
+
+    int CF1_count;
+    int CF2_count;
+    int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
+    int stat_CF2;
+    int stat_Fr;
+    int stat_Gr;
+
+    int F_recur_count;
+    int G_recur_count;
+
+    double err_amplify;
+
+    F_lam_F  = F_sign_lam_F * SMALL;  /* unnormalized */
+    Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
+
+    /* Backward recurrence to get F,Fp at lam_min */
+    F_recur_count = GSL_MAX(k_lam_G, N);
+    stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x,
+                              F_lam_F, Fp_lam_F,
+                              &F_lam_min_unnorm, &Fp_lam_min_unnorm
+                              );
+    Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
+
+    /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */
+    stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count);
+    alpha = Fp_over_F_lam_min - P_lam_min;
+    gamma = alpha/Q_lam_min;
+
+    F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ;
+
+    F_lam_min  = F_sign_lam_min / sqrt(alpha*alpha/Q_lam_min + Q_lam_min);
+    Fp_lam_min = Fp_over_F_lam_min * F_lam_min;
+    G_lam_min  = gamma * F_lam_min;
+    Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min;
+
+    /* Apply scale to values of F,Fp at lam_F (the top). */
+    F_scale = F_lam_min / F_lam_min_unnorm;    
+    F_lam_F  *= F_scale;
+    Fp_lam_F *= F_scale;
+
+    /* Forward recurrence to get G,Gp at lam_G (the top). */
+    G_recur_count = GSL_MAX(N-k_lam_G,0);
+    stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x,
+                              G_lam_min, Gp_lam_min,
+                              &G_lam_G, &Gp_lam_G
+                              );
+
+    err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1;
+
+    F->val  = F_lam_F;
+    F->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val);
+
+    Fp->val = Fp_lam_F;
+    Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val);
+
+    G->val  = G_lam_G;
+    G->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val);
+
+    Gp->val = Gp_lam_G;
+    Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val);
+
+    *exp_F = 0.0;
+    *exp_G = 0.0;
+
+    return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr);
+  }
+}
+
+
+int
+gsl_sf_coulomb_wave_F_array(double lam_min, int kmax,
+                                 double eta, double x, 
+                                 double * fc_array,
+                                 double * F_exp)
+{
+  if(x == 0.0) {
+    int k;
+    *F_exp = 0.0;
+    for(k=0; k<=kmax; k++) {
+      fc_array[k] = 0.0;
+    }
+    if(lam_min == 0.0){
+      gsl_sf_result f_0;
+      CLeta(0.0, eta, &f_0);
+      fc_array[0] = f_0.val;
+    }
+    return GSL_SUCCESS;
+  }
+  else {
+    const double x_inv = 1.0/x;
+    const double lam_max = lam_min + kmax;
+    gsl_sf_result F, Fp;
+    gsl_sf_result G, Gp;
+    double G_exp;
+
+    int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0,
+                                              &F, &Fp, &G, &Gp, F_exp, &G_exp);
+
+    double fcl  = F.val;
+    double fpl = Fp.val;
+    double lam = lam_max;
+    int k;
+
+    fc_array[kmax] = F.val;
+
+    for(k=kmax-1; k>=0; k--) {
+      double el = eta/lam;
+      double rl = hypot(1.0, el);
+      double sl = el  + lam*x_inv;
+      double fc_lm1 = (fcl*sl + fpl)/rl;
+      fc_array[k]   = fc_lm1;
+      fpl           =  fc_lm1*sl - fcl*rl;
+      fcl           =  fc_lm1;
+      lam -= 1.0;
+    }
+
+    return stat_FG;
+  }
+}
+
+
+int
+gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax,
+                                  double eta, double x,
+                                  double * fc_array, double * gc_array,
+                                  double * F_exp, double * G_exp)
+{
+  const double x_inv = 1.0/x;
+  const double lam_max = lam_min + kmax;
+  gsl_sf_result F, Fp;
+  gsl_sf_result G, Gp;
+
+  int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
+                                            &F, &Fp, &G, &Gp, F_exp, G_exp);
+
+  double fcl  = F.val;
+  double fpl = Fp.val;
+  double lam = lam_max;
+  int k;
+
+  double gcl, gpl;
+
+  fc_array[kmax] = F.val;
+
+  for(k=kmax-1; k>=0; k--) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el  + lam*x_inv;
+    double fc_lm1;
+    fc_lm1 = (fcl*sl + fpl)/rl;
+    fc_array[k] = fc_lm1;
+    fpl         =  fc_lm1*sl - fcl*rl;
+    fcl         =  fc_lm1;
+    lam -= 1.0;
+  }
+
+  gcl = G.val;
+  gpl = Gp.val;
+  lam = lam_min + 1.0;
+
+  gc_array[0] = G.val;
+
+  for(k=1; k<=kmax; k++) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el + lam*x_inv;
+    double gcl1 = (sl*gcl - gpl)/rl;
+    gc_array[k] = gcl1;
+    gpl         = rl*gcl - sl*gcl1;
+    gcl         = gcl1;
+    lam += 1.0;
+  }
+
+  return stat_FG;
+}
+
+
+int
+gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax,
+                                   double eta, double x,
+                                   double * fc_array, double * fcp_array,
+                                   double * gc_array, double * gcp_array,
+                                   double * F_exp, double * G_exp)
+
+{
+  const double x_inv = 1.0/x;
+  const double lam_max = lam_min + kmax;
+  gsl_sf_result F, Fp;
+  gsl_sf_result G, Gp;
+
+  int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
+                                            &F, &Fp, &G, &Gp, F_exp, G_exp);
+
+  double fcl  = F.val;
+  double fpl = Fp.val;
+  double lam = lam_max;
+  int k;
+
+  double gcl, gpl;
+
+  fc_array[kmax]  = F.val;
+  fcp_array[kmax] = Fp.val;
+
+  for(k=kmax-1; k>=0; k--) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el  + lam*x_inv;
+    double fc_lm1;
+    fc_lm1 = (fcl*sl + fpl)/rl;
+    fc_array[k]  = fc_lm1;
+    fpl          = fc_lm1*sl - fcl*rl;
+    fcp_array[k] = fpl;
+    fcl          =  fc_lm1;
+    lam -= 1.0;
+  }
+
+  gcl = G.val;
+  gpl = Gp.val;
+  lam = lam_min + 1.0;
+
+  gc_array[0]  = G.val;
+  gcp_array[0] = Gp.val;
+
+  for(k=1; k<=kmax; k++) {
+    double el = eta/lam;
+    double rl = hypot(1.0, el);
+    double sl = el + lam*x_inv;
+    double gcl1 = (sl*gcl - gpl)/rl;
+    gc_array[k]  = gcl1;
+    gpl          = rl*gcl - sl*gcl1;
+    gcp_array[k] = gpl;
+    gcl          = gcl1;
+    lam += 1.0;
+  }
+
+  return stat_FG;
+}
+
+
+int
+gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
+                                    double eta, double x,
+                                    double * fc_array,
+                                    double * F_exp)
+{
+  if(x < 0.0 || lam_min < -0.5) {
+    GSL_ERROR ("domain error", GSL_EDOM);
+  }
+  else if(x < 10.0/GSL_DBL_MAX) {
+    int k;
+    for(k=0; k<=kmax; k++) {
+      fc_array[k] = 0.0;
+    }
+    if(lam_min == 0.0) {
+      fc_array[0] = sqrt(C0sq(eta));
+    }
+    *F_exp = 0.0;
+    if(x == 0.0)
+      return GSL_SUCCESS;
+    else
+      GSL_ERROR ("underflow", GSL_EUNDRFLW);
+  }
+  else {
+    int k;
+    int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax,
+                                                  eta, x, 
+                                                  fc_array,
+                                                  F_exp);
+
+    for(k=0; k<=kmax; k++) {
+      fc_array[k] = fc_array[k] / x;
+    }
+    return stat_F;
+  }
+}
+
+