glibc/sysdeps/libm-ieee754/e_lgamma_r.c
Ulrich Drepper d705269e75 Update.
1997-06-19 19:38  Ulrich Drepper  <drepper@cygnus.com>

	* features.h: Define __STDC_IEC_559__ and _STDC_IEC_559_COMPLEX__.

	* elf/dl-minimal.c (__dcgettext): Remove assertion.

	* inet/rcmd.c: Correct a few typos.  Reported by Eric Troan.

	* manual/Makefile (distribute): Add dir.
	* manual/dir: New file.

	* math/Makefile (libm-support): Rename s_rinttol, s_rinttoll,
	s_roundtol, and s_roundtoll to s_lrint, s_llrint, s_lround,
	and s_llround respectively.
	(libm-calls): Add e_gamma_r.
	* math/libm-test.c (check_int_exc): New function.
	(signbit_test): Remove test for sign of NaN.
	(gamma_test): Clear exception after test of existence.
	Correct tests and and epsilons.
	(lgamma_test): Likewise.
	(ilogb_test): Correct all tests.
	(scalb_test): Rewrite.
	(rinttol_test): Rename to lrint_test and correct tests.
	(rinttoll_test): Likewise.
	(roundtol_test): Likewise.
	(roundtoll_test): Likewise.
	(main): Call lrint/lround functions instead of rinttol/roundtol.
	* math/math.h: Change prototypes for rinttol/roundtol.
	* math/mathcalls.h: Rearrange prototypes according to ISO C9X draft.
	* sysdeps/generic/mathbits.h: Define FP_ILOGB0 and FP_ILOGBNAN.
	* sysdeps/i386/mathbits.h: Likewise.
	* sysdeps/libm-i387/e_scalb.S: Handle special cases correctly.
	* sysdeps/libm-i387/e_scalbf.S: Likewise.
	* sysdeps/libm-i387/e_scalbl.S: Likewise.
	* sysdeps/libm-i387/s_asinh.S: Handle -inf correctly.
	* sysdeps/libm-i387/s_asinhf.S: Likewise.
	* sysdeps/libm-i387/s_asinhl.S: Likewise.
	* sysdeps/libm-i387/s_ilogb.S: Optimize.
	* sysdeps/libm-i387/s_ilogbf.S: Likewise.
	* sysdeps/libm-i387/s_ilogbl.S: Likewise.
	* sysdeps/libm-i387/s_rinttol.S: Rename to...
	* sysdeps/libm-i387/s_lrint.S: ...this.
	* sysdeps/libm-i387/s_rinttoll.S: Rename to...
	* sysdeps/libm-i387/s_llrint.S: ...this.
	* sysdeps/libm-i387/s_remquo.S: Correctly set sign of remainder.
	* sysdeps/libm-i387/s_remquof.S: Likewise.
	* sysdeps/libm-i387/s_remquol.S: Likewise.
	* sysdeps/libm-i387/e_gamma_r.c: New file.  Implementation of gamma
	function according to ISO C.
	* sysdeps/libm-i387/e_gammaf_r.c: New file.
	* sysdeps/libm-i387/e_gammal_r.c: New file.
	* sysdeps/libm-i387/e_lgamma_r.c: Don't let optimize compile the
	generation of exceptions away.
	* sysdeps/libm-i387/e_lgammaf_r.c: Likewise.
	* sysdeps/libm-i387/k_standard.c: Correct return value for infinity
	points of gamma function when not SVID mode.
	* sysdeps/libm-i387/s_rinttoll.c: Renamed to...
	* sysdeps/libm-i387/s_llrint.c: ...this.
	* sysdeps/libm-i387/s_rinttol.c: Renamed to...
	* sysdeps/libm-i387/s_lrint.c: ...this.
	* sysdeps/libm-i387/s_roundtoll.c: Renamed to...
	* sysdeps/libm-i387/s_llround.c: ...this.
	* sysdeps/libm-i387/s_roundtol.c: Renamed to...
	* sysdeps/libm-i387/s_lround.c: ..this.
	* sysdeps/libm-i387/s_scalbn.c: Change second parameter according to
	ISO C.
	* sysdeps/libm-i387/s_scalbnf.c: Likewise.
	* sysdeps/libm-i387/s_scalbnl.c: Likewise.
	* sysdeps/libm-i387/w_gamma.c: Call __ieee754_gamma_r if library
	mode is not _SVID_.
	* sysdeps/libm-i387/w_gammaf.c: Likewise.
	* sysdeps/libm-i387/w_gammal.c: Likewise.
	* sysdeps/m68k/fpu/__math.h: Rename __rinttol to __lrint and
	rinttol to lrint.
	* sysdeps/m68k/fpu/s_rinttol.c: Renamed to...
	* sysdeps/m68k/fpu/s_lrint.c: ...this.
	* sysdeps/m68k/fpu/s_rinttoll.c: Renamed to...
	* sysdeps/m68k/fpu/s_llrint.c: ...this.

	* md5-crypt/Makefile: Link md5test program with md5.o.

	* stdio-common/temptest.c: Don't use __stdio_gen_tempname which
	is not exported by the libc.so.

	* stdio-common/vfscanf.c: Correct scanning of strings after last
	change.

	* sysdeps/unix/sysv/linux/i386/sysdep.S: Use .comm to define errno.

1997-06-19 07:37  H.J. Lu  <hjl@gnu.ai.mit.edu>

	* time/tzfile.c (__tzfile_read): Store getc () return in int.

1997-06-13  Miles Bader  <miles@gnu.ai.mit.edu>

	* argp/argp-parse.c (argp_version_parser): Include `(PROGRAM ERROR)' in
	the no-version error text to indicate that something's fucked.
	[!_] (N_): New macro.
	(argp_default_options, argp_version_options): Wrap doc strings in N_().

1997-06-12  Miles Bader  <miles@gnu.ai.mit.edu>

	* argp/argp-parse.c (parser_parse_opt): Detect and report unhandled
	options here.
	(parser_parse_arg): Handle ARGP_KEY_ARGS here.
	Adjust NEXT pointer back if we fail to parse anything.
	(parser_parse_next): Simplify arg code.  Leave state NEXT frobbing
	to parser_parse_arg.

1997-06-11  Miles Bader  <miles@gnu.ai.mit.edu>

	* argp/argp.h (ARGP_KEY_ARGS, ARGP_KEY_FINI): New macros.
	* argp/argp-parse.c (parser_finalize): Do another pass over the
	parsers with ARGP_KEY_FINI.

1997-06-18  Miles Bader  <miles@gnu.ai.mit.edu>

	* string/Makefile (routines): Add argz-replace.

1997-06-16 00:16  Miles Bader  <miles@gnu.ai.mit.edu>

	* manual/string.texi (Argz Functions): Document argz_replace.

1997-06-12  Miles Bader  <miles@gnu.ai.mit.edu>

	* string/argz.h (__argz_replace, argz_replace): New declarations.
	* string/argz-replace.c: New file.

1997-06-11  Miles Bader  <miles@gnu.ai.mit.edu>

	* manual/argp.texi (Argp Special Keys): Document ARGP_KEY_ARGS.

1997-06-16 23:08  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* Makeconfig (libc-map): Remove definition.
	* Makerules (libc-map): Define it here, using the full name.
	(load-map-file): Remove case for empty $(..).
	($(common-objpfx)libc.so): Revert last change.

1997-06-17 22:18  Mark Kettenis  <kettenis@phys.uva.nl>

	* login/programs/utmpd.c (main): Improve signal handling.

	* login/programs/request.c (do_setutent, do_updwtmp):
	Allow arbitrary length filenames.
	* login/programs/utmpd.h (setutent_request, updwtmp_request):
	Get rid of fixed length file field.
	* login/utmp_daemon.c (do_setutent, do_updwtmp):
	Allow arbitrary	length filenames.

	* login/programs/request.c (do_pututline):
	Don't fail if connection->position is -1 on entry.

1997-06-15 16:32  Mark Kettenis  <kettenis@phys.uva.nl>

	* login/utmp_file.c (updwtmp_file): Use the same method for
	appending an entry as in pututline_file.

1997-06-11 18:59  Mark Kettenis  <kettenis@phys.uva.nl>

	* login/utmp_file.c (getutent_r_file):
	Use read lock instead of write lock.
	(getutline_r_file, internal_getut_r): Lock utmp file.
	(updwtmp_file): Use fcntl to lock file instead of flock.

1997-06-18 00:11  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* stdio-common/vfscanf.c (inchar, ungetc): Don't count EOF as
	character read in/put back.
	* stdio-common/tstscanf.c: Add test case for this.

1997-06-17 22:17  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* libc.map: Add more symbols.

1997-06-18 12:01  Ulrich Drepper  <drepper@cygnus.com>

	* manual/Makefile (install): Make sure `dir' file exists if we use
	install-info.

1997-06-17 19:32  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* manual/Makefile (info): Depend on dir-add.info.

1997-06-17 17:19  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* elf/ldd.bash.in: Fix spacing in message.

1997-06-17 14:28  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* elf/do-rel.h (elf_dynamic_do_rel): Always use version
	information if available.

1997-06-17 11:34  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sunrpc/Makefile ($(rpcsvc:%.x=$(objpfx)rpcsvc/%.h)): Make
	command non-empty to force make to recheck modification time.
	($(rpcsvc:%.x=$(objpfx)x%.c)): Likewise.

1997-06-17 00:26  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/stub/e_acoshl.c: Set errno to ENOSYS.
	* sysdeps/stub/e_acosl.c: Likewise.
	* sysdeps/stub/e_asinl.c: Likewise.
	* sysdeps/stub/e_atan2l.c: Likewise.
	* sysdeps/stub/e_expl.c: Likewise.
	* sysdeps/stub/e_fmodl.c: Likewise.
	* sysdeps/stub/e_j0l.c: Likewise.
	* sysdeps/stub/e_j1l.c: Likewise.
	* sysdeps/stub/e_jnl.c: Likewise.
	* sysdeps/stub/e_lgammal_r.c: Likewise.
	* sysdeps/stub/e_log10l.c: Likewise.
	* sysdeps/stub/e_logl.c: Likewise.
	* sysdeps/stub/e_powl.c: Likewise.
	* sysdeps/stub/e_rem_pio2l.c: Likewise.
	* sysdeps/stub/e_sqrtl.c: Likewise.
	* sysdeps/stub/k_cosl.c: Likewise.
	* sysdeps/stub/k_rem_pio2l.c: Likewise.
	* sysdeps/stub/k_sinl.c: Likewise.
	* sysdeps/stub/k_tanl.c: Likewise.
	* sysdeps/stub/s_atanl.c: Likewise.
	* sysdeps/stub/s_erfl.c: Likewise.
	* sysdeps/stub/s_exp2.c: Likewise.
	* sysdeps/stub/s_exp2f.c: Likewise.
	* sysdeps/stub/s_exp2l.c: Likewise.
	* sysdeps/stub/s_expm1l.c: Likewise.
	* sysdeps/stub/s_log1pl.c: Likewise.
	* sysdeps/stub/s_log2l.c: Likewise.

1997-06-18 11:46  Ulrich Drepper  <drepper@cygnus.com>

	* sysdeps/libm-ieee754/w_gamma.c: If _LIB_VERSION is _SVID_ compute
	result as before last change.
	* sysdeps/libm-ieee754/w_gammaf.c: Likewise.
	* sysdeps/libm-ieee754/w_gammal.c: Likewise.

1997-06-16 23:37  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/libm-ieee754/s_remquo.c: Fix off-by-one when computing
	quotient.
	* sysdeps/libm-ieee754/s_remquof.c: Likewise.
	* sysdeps/libm-ieee754/s_remquol.c: Likewise.

	* sysdeps/m68k/fpu/s_remquo.c: Remove FIXME and special case for
	quotient.

	* sysdeps/libm-ieee754/w_gamma.c: Add missing call to exp
	function.  Don't use global signgam.
	* sysdeps/libm-ieee754/w_gammaf.c: Likewise.
	* sysdeps/libm-ieee754/w_gammal.c: Likewise.

	* math/Makefile (libm-calls): Remove w_gamma_r.
	* sysdeps/libm-ieee754/w_gamma_r.c: Remove file.
	* sysdeps/libm-ieee754/w_gammaf_r.c: Remove file.
	* sysdeps/libm-ieee754/w_gammal_r.c: Remove file.

	* math/libm-test.c (atanh_test): Declare x only if needed.
	(signbit_test): Fix typo.
	(gamma_test): Check whether function is implemented.  Add
	epsilons.
	(lgamma_test): Likewise.
	(fmod_test): Add epsilons.
	(exp2_test): Use right function for existence test.

1997-06-07 09:20  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/powerpc/Dist: Add fenv_const.c.
	* sysdeps/unix/sysv/linux/Dist: Add net/if_slip.h.
	* sysdeps/unix/sysv/linux/powerpc/Dist: Add init-first.h and
	syscall.h.
	* sysdeps/unix/sysv/linux/sparc/Dist: Add init-first.h.
	* string/Makefile (distribute): Add tst-svc.expect.
	* nis/Makefile (distribute): Add nis_intern.h and Banner.
	* elf/Makefile (distribute): Add dl-hash.h.
	* Rules (subdir_echo-distinfo): Add headers from $(distribute).

	* login/Makefile (others): Add utmpdump.
	* sysdeps/unix/sysv/linux/Makefile (sysdep_headers): Add
	net/if_slip.h.

	* manual/Makefile (dir-add.texi): Also look in indirectly included
	files.

1997-06-16 23:15  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/unix/sysv/linux/configure.in: Remove commands to
	generate stdio_lim.h.
	* sysdeps/unix/sysv/linux/mk-stdiolim.c: Remove.
	* sysdeps/unix/sysv/linux/Makefile: Add rules to generate
	stdio_lim.h here.
	(common-generated): Add generated files.
	(inhibit-stdio_lim): Define.
	* sysdeps/posix/Makefile [$(inhibit-stdio_lim)=yes]: Disable rules
	to generate stdio_lim.h.

	* sysdeps/unix/sysv/linux/Makefile: Suppress inclusion of
	dependecy files if no_deps is set.
	($(objpfx)syscall-%.d): Add header file as target to dependency
	generation.

1997-06-14 19:19  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/unix/sysv/linux/arm/Dist: New file.

1997-06-14 17:59  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* extra-lib.mk (others): Depend on versioned shared library, not
	the unversioned one.
	* Makerules (build-shlib): Don't make the version link here.
	($(common-objpfx)libc.so$(libc.so-version)): New rule for libc
	version link.

1997-06-16 03:07  Ulrich Drepper  <drepper@cygnus.com>

	* sysdeps/libm-ieee754/k_standard.c: Undo change of Tue Aug 6
	01:13:56 1996.

	* argp/argp-help.c (argp_args_usage): Supply correct argp to
	filter_doc.
	* argp/argp-help.c (hol_add_cluster): Initialize CL->depth.
	* argp/argp-help.c (_help): Supply STATE to argp_args_usage.
	* argp/argp.h (ARGP_KEY_HELP_ARGS_DOC): New macro.
	* argp/argp-fmtstream.c: Add casts to prevent warnings.
	* argp/argp.h (OPTION_NO_USAGE): New macro.
	* argp/argp-help.c (usage_long_opt, usage_argful_short_opt,
	* argp/argp-fmtstream.c (__argp_fmtstream_update): Account for case
	* argp/argp-help.c <stddef.h>: New include.
	* argp/argp.h (argp_state_help, __argp_state_help, argp_usage,
	* argp/argp.h (argp_program_bug_address): Make const.
	* argp/argp-parse.c (argp_default_parser): Set STATE->name for OPT_PROGNAME.
	* argp/argp-help.c (__argp_error, __argp_failure, __argp_state_help):
	* argp/argp-parse.c (parser_init): Set PARSER->state.flags.
1997-06-19 19:11:22 +00:00

315 lines
11 KiB
C

/* @(#)er_lgamma.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $";
#endif
/* __ieee754_lgamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3)
* = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
* 2. Polynomial approximation of lgamma around its
* minimun ymin=1.461632144968362245 to maintain monotonicity.
* On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
* Let z = x-ymin;
* lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
* where
* poly(z) is a 14 degree polynomial.
* 2. Rational approximation in the primary interval [2,3]
* We use the following approximation:
* s = x-2.0;
* lgamma(x) = 0.5*s + s*P(s)/Q(s)
* with accuracy
* |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
* Our algorithms are based on the following observation
*
* zeta(2)-1 2 zeta(3)-1 3
* lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
* 2 3
*
* where Euler = 0.5771... is the Euler constant, which is very
* close to 0.5.
*
* 3. For x>=8, we have
* lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
* (better formula:
* lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
* Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
* where
* |w - f(z)| < 2**-58.74
*
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
* Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
*
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1)=lgamma(2)=0
* lgamma(x) ~ -log(x) for tiny x
* lgamma(0) = lgamma(inf) = inf
* lgamma(-integer) = +-inf
*
*/
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const double
#else
static double
#endif
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
/* tt = -(tail of tf) */
tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
#ifdef __STDC__
static const double zero= 0.00000000000000000000e+00;
#else
static double zero= 0.00000000000000000000e+00;
#endif
#ifdef __STDC__
static double sin_pi(double x)
#else
static double sin_pi(x)
double x;
#endif
{
double y,z;
int n,ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
y = -x; /* x is assume negative */
/*
* argument reduction, make sure inexact flag not raised if input
* is an integer
*/
z = __floor(y);
if(z!=y) { /* inexact anyway */
y *= 0.5;
y = 2.0*(y - __floor(y)); /* y = |x| mod 2.0 */
n = (int) (y*4.0);
} else {
if(ix>=0x43400000) {
y = zero; n = 0; /* y must be even */
} else {
if(ix<0x43300000) z = y+two52; /* exact */
GET_LOW_WORD(n,z);
n &= 1;
y = n;
n<<= 2;
}
}
switch (n) {
case 0: y = __kernel_sin(pi*y,zero,0); break;
case 1:
case 2: y = __kernel_cos(pi*(0.5-y),zero); break;
case 3:
case 4: y = __kernel_sin(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cos(pi*(y-1.5),zero); break;
default: y = __kernel_sin(pi*(y-2.0),zero,0); break;
}
return -y;
}
#ifdef __STDC__
double __ieee754_lgamma_r(double x, int *signgamp)
#else
double __ieee754_lgamma_r(x,signgamp)
double x; int *signgamp;
#endif
{
double t,y,z,nadj,p,p1,p2,p3,q,r,w;
int i,hx,lx,ix;
EXTRACT_WORDS(hx,lx,x);
/* purge off +-inf, NaN, +-0, and negative arguments */
*signgamp = 1;
if ((unsigned int) hx==0xfff00000&&lx==0)
return x-x;
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x*x;
if((ix|lx)==0) return one/fabs(x);
if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx<0) {
*signgamp = -1;
return -__ieee754_log(-x);
} else return -__ieee754_log(x);
}
if(hx<0) {
if(ix>=0x43300000) /* |x|>=2**52, must be -integer */
return x/zero;
t = sin_pi(x);
if(t==zero) return one/fabsf(t); /* -integer */
nadj = __ieee754_log(pi/fabs(t*x));
if(t<zero) *signgamp = -1;
x = -x;
}
/* purge off 1 and 2 */
if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
/* for x < 2.0 */
else if(ix<0x40000000) {
if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -__ieee754_log(x);
if(ix>=0x3FE76944) {y = one-x; i= 0;}
else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
else {y = x; i=2;}
} else {
r = zero;
if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
else {y=x-one;i=2;}
}
switch(i) {
case 0:
z = y*y;
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
p = y*p1+p2;
r += (p-0.5*y); break;
case 1:
z = y*y;
w = z*y;
p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
p = z*p1-(tt-w*(p2+y*p3));
r += (tf + p); break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += (-0.5*y + p1/p2);
}
}
else if(ix<0x40200000) { /* x < 8.0 */
i = (int)x;
t = zero;
y = x-(double)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = half*y+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+6.0); /* FALLTHRU */
case 6: z *= (y+5.0); /* FALLTHRU */
case 5: z *= (y+4.0); /* FALLTHRU */
case 4: z *= (y+3.0); /* FALLTHRU */
case 3: z *= (y+2.0); /* FALLTHRU */
r += __ieee754_log(z); break;
}
/* 8.0 <= x < 2**58 */
} else if (ix < 0x43900000) {
t = __ieee754_log(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
} else
/* 2**58 <= x <= inf */
r = x*(__ieee754_log(x)-one);
if(hx<0) r = nadj - r;
return r;
}