diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 9b64674bbe7d..4b164545f0b8 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -1,287 +1,287 @@
 #  @(#)Makefile 5.1beta 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.
 #  ====================================================
 #
 #
 
 PACKAGE=	clibs
 ARCH_SUBDIR=	${MACHINE_CPUARCH:S/i386/i387/}
 
 .include "${ARCH_SUBDIR}/Makefile.inc"
 
 .PATH:	${.CURDIR}/${ARCH_SUBDIR}
 .if ${MACHINE_CPUARCH} == "i386" || ${MACHINE_CPUARCH} == "amd64"
 .PATH:	${.CURDIR}/x86
 CFLAGS+=	-I${.CURDIR}/x86
 .endif
 
 # long double format
 .if ${LDBL_PREC} == 64
 .PATH:  ${.CURDIR}/ld80
 CFLAGS+=	-I${.CURDIR}/ld80
 .elif ${LDBL_PREC} == 113
 .PATH:  ${.CURDIR}/ld128
-CFLAGS+=	-I${.CURDIR}/ld128
+CFLAGS+=	-I${.CURDIR}/ld128 -I${SRCTOP}/contrib/arm-optimized-routines/math
 .endif
 
 CFLAGS+=	-I${.CURDIR}/${ARCH_SUBDIR}
 
 .include <bsd.compiler.mk>
 .if ${COMPILER_TYPE} == "clang" && ${COMPILER_VERSION} >= 100000 && \
     (${MACHINE_CPUARCH} == "amd64" || ${MACHINE_CPUARCH} == "i386")
 # When using clang with x86_64 CPUs that support AVX, some floating point
 # transformations may raise exceptions that would not have been raised by the
 # original code. To avoid this, use the -fp-exception-behavior=maytrap flag,
 # introduced in clang 10.0.0.
 # See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=254911
 CFLAGS+=	-ffp-exception-behavior=maytrap
 .endif
 
 # Tell the compiler we don't set errno in any of the math functions. This
 # stops gcc from trying to generate a function call to set errno when using
 # a builtin in the implementation
 CFLAGS+=	-fno-math-errno
 
 .PATH:	${.CURDIR}/bsdsrc
 .PATH:	${.CURDIR}/src
 .PATH:	${.CURDIR}/man
 
 LIB=	m
 SHLIBDIR?= /lib
 SHLIB_MAJOR= 5
 WARNS?=	1
 IGNORE_PRAGMA=
 COMMON_SRCS= b_tgamma.c \
 	e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
 	e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
 	e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \
 	e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
 	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
 	e_pow.c e_powf.c e_rem_pio2.c \
 	e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
 	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
 	k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
 	k_tan.c k_tanf.c \
 	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
 	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c s_clog.c s_clogf.c \
 	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
 	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
 	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
 	s_finite.c s_finitef.c \
 	s_floor.c s_floorf.c s_fma.c s_fmaf.c \
 	s_fmax.c s_fmaxf.c s_fmin.c \
 	s_fminf.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
 	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
 	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
 	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
 	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
 	s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \
 	s_nexttowardf.c s_remquo.c s_remquof.c \
 	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c \
 	s_sincos.c s_sincosf.c s_sinf.c \
 	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
 	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
 
 # IEEE-754 2008 and ISO/IEC TS 18661-4 half-cycle trignometric functions
 COMMON_SRCS+= s_cospi.c s_cospif.c \
 	s_sinpi.c s_sinpif.c \
 	s_tanpi.c s_tanpif.c
 
 # LIBCSRCDIR is defined lazily, so can't be used for exists() below
 LIBC_SRCTOP=${SRCTOP}/lib/libc
 # Location of fpmath.h and _fpmath.h
 
 # Pick the current architecture directory for libc. In general, this is named
 # MACHINE_CPUARCH, but some ABIs are different enough to require their own libc,
 # so allow a directory named MACHINE_ARCH to override this (though treat
 # powerpc64le and powerpc64 the same).
 # Note: This is copied from libc/Makefile
 M=${MACHINE_ARCH:S/powerpc64le/powerpc64/}
 .if exists(${LIBC_SRCTOP}/${M})
 LIBC_ARCH=${M}
 .else
 LIBC_ARCH=${MACHINE_CPUARCH}
 .endif
 
 CFLAGS+=	-I${.CURDIR}/src -I${LIBC_SRCTOP}/include \
 		-I${LIBC_SRCTOP}/${LIBC_ARCH}
 SYM_MAPS+=	${.CURDIR}/Symbol.map
 
 VERSION_DEF=	${LIBC_SRCTOP}/Versions.def
 SYMBOL_MAPS=	${SYM_MAPS}
 
 # C99 long double functions
 COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 .if ${LDBL_PREC} != 53
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	b_tgammal.c catrigl.c \
 	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
 	e_coshl.c e_fmodl.c e_hypotl.c \
 	e_lgammal.c e_lgammal_r.c e_powl.c \
 	e_remainderl.c e_sinhl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
 	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cexpl.c \
 	s_clogl.c s_cosl.c s_cospil.c s_cprojl.c \
 	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
 	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
 	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
 	s_scalbnl.c s_sinl.c s_sincosl.c s_sinpil.c \
 	s_tanhl.c s_tanl.c s_tanpil.c s_truncl.c w_cabsl.c
 # Work around this warning from gcc:
 #     lib/msun/ld80/e_powl.c:275:1: error: floating constant exceeds range of
 #     'long double' [-Werror=overflow]
 #     if( y >= LDBL_MAX )
 # See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=130067
 .if ${COMPILER_TYPE} == "gcc"
 CFLAGS.e_powl.c+= -Wno-error=overflow
 .endif
 .endif
 
 # C99 complex functions
 COMMON_SRCS+=	catrig.c catrigf.c \
 	s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \
 	s_cimag.c s_cimagf.c s_cimagl.c \
 	s_conj.c s_conjf.c s_conjl.c \
 	s_cpow.c s_cpowf.c s_cpowl.c \
 	s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c \
 	s_csinh.c s_csinhf.c s_ctanh.c s_ctanhf.c
 
 # FreeBSD's C library supplies these functions (but note we still have s_frexp.c
 # and s_isnan.c enabled above, so they are duplicated):
 #COMMON_SRCS+=	s_fabs.c s_frexp.c s_isnan.c s_ldexp.c s_modf.c
 
 # Exclude the generic versions of what we provide in the MD area.
 .if defined(ARCH_SRCS) && "${MK_MACHDEP_OPTIMIZATIONS}" != no
 .for i in ${ARCH_SRCS}
 COMMON_SRCS:=  ${COMMON_SRCS:N${i:R}.c}
 .endfor
 .endif
 
 SRCS=	${COMMON_SRCS} ${ARCH_SRCS}
 
 INCS+=	fenv.h math.h
 
 MAN=	acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
 	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
 	cimag.3 clog.3 copysign.3 cos.3 cosh.3 cospi.3 \
 	cpow.3 csqrt.3 erf.3 \
 	exp.3 fabs.3 fdim.3 \
 	feclearexcept.3 feenableexcept.3 fegetenv.3 \
 	fegetround.3 fenv.3 floor.3 \
 	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
 	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
 	nextafter.3 remainder.3 rint.3 \
 	round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
 	sinh.3 sinpi.3 sqrt.3 tan.3 tanh.3 tanpi.3 trunc.3 \
 	complex.3
 
 MLINKS+=acos.3 acosf.3 acos.3 acosl.3
 MLINKS+=acosh.3 acoshf.3 acosh.3 acoshl.3
 MLINKS+=asin.3 asinf.3 asin.3 asinl.3
 MLINKS+=asinh.3 asinhf.3 asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 atan.3 atanl.3
 MLINKS+=atanh.3 atanhf.3 atanh.3 atanhl.3
 MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \
 	atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3
 MLINKS+=cacos.3 cacosf.3 cacos.3 cacosl.3 \
 	cacos.3 cacosh.3 cacos.3 cacoshf.3 cacos.3 cacoshl.3 \
 	cacos.3 casin.3 cacos.3 casinf.3 cacos.3 casinl.3 \
 	cacos.3 casinh.3 cacos.3 casinhf.3 cacos.3 casinhl.3 \
 	cacos.3 catan.3 cacos.3 catanf.3 cacos.3 catanl.3 \
 	cacos.3 catanh.3 cacos.3 catanhf.3 cacos.3 catanhl.3
 MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 ccos.3 ctan.3 ccos.3 ctanf.3
 MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \
 	ccosh.3 ctanh.3 ccosh.3 ctanhf.3
 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3
 MLINKS+=cexp.3 cexpf.3 cexp.3 cexpl.3
 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
 	cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \
 	cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \
 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
 MLINKS+=clog.3 clogf.3 clog.3 clogl.3
 MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
 MLINKS+=cos.3 cosf.3 cos.3 cosl.3
 MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
 MLINKS+=cospi.3 cospif.3 cospi.3 cospil.3
 MLINKS+=cpow.3 cpowf.3 cpow.3 cpowl.3
 MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \
 	exp.3 powl.3 exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 \
 	exp.3 expl.3
 MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
 MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3
 MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \
 	feclearexcept.3 fesetexceptflag.3 feclearexcept.3 fetestexcept.3
 MLINKS+=feenableexcept.3 fedisableexcept.3 feenableexcept.3 fegetexcept.3
 MLINKS+=fegetenv.3 feholdexcept.3 fegetenv.3 fesetenv.3 \
 	fegetenv.3 feupdateenv.3
 MLINKS+=fegetround.3 fesetround.3
 MLINKS+=floor.3 floorf.3 floor.3 floorl.3
 MLINKS+=fma.3 fmaf.3 fma.3 fmal.3
 MLINKS+=fmax.3 fmaxf.3 fmax.3 fmaxl.3 \
 	fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3
 MLINKS+=fmod.3 fmodf.3 fmod.3 fmodl.3
 MLINKS+=hypot.3 cabs.3 hypot.3 cabsf.3 hypot.3 cabsl.3 \
 	hypot.3 hypotf.3 hypot.3 hypotl.3
 MLINKS+=ieee_test.3 scalb.3 ieee_test.3 scalbf.3
 MLINKS+=ieee_test.3 significand.3 ieee_test.3 significandf.3
 MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \
 	ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3
 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \
 	lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \
 	lgamma.3 tgamma.3 lgamma.3 tgammaf.3
 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \
 	log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \
 	log.3 logf.3 log.3 logl.3 \
 	log.3 log2.3 log.3 log2f.3 log.3 log2l.3
 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 llrintl.3 \
 	lrint.3 lrintf.3 lrint.3 lrintl.3
 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \
 	lround.3 lroundf.3 lround.3 lroundl.3
 MLINKS+=nan.3 nanf.3 nan.3 nanl.3
 MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3
 MLINKS+=nextafter.3 nexttoward.3 nextafter.3 nexttowardf.3
 MLINKS+=nextafter.3 nexttowardl.3
 MLINKS+=remainder.3 remainderf.3 remainder.3 remainderl.3 \
 	remainder.3 remquo.3 remainder.3 remquof.3 remainder.3 remquol.3
 MLINKS+=rint.3 rintf.3 rint.3 rintl.3 \
 	rint.3 nearbyint.3 rint.3 nearbyintf.3 rint.3 nearbyintl.3
 MLINKS+=round.3 roundf.3 round.3 roundl.3
 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 sin.3 sinl.3
 MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
 MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
 MLINKS+=sinpi.3 sinpif.3 sinpi.3 sinpil.3
 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
 	sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 tan.3 tanl.3
 MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
 MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
 .include <src.opts.mk>
 
 HAS_TESTS=
 SUBDIR.${MK_TESTS}+= tests
 
 .include <bsd.compiler.mk>
 .if ${COMPILER_TYPE} == "clang" && ${COMPILER_VERSION} >= 120000
 # Silence '#pragma FENV_ACCESS' is not supported on this target - ignored
 CWARNFLAGS+=	-Wno-error=ignored-pragmas
 .endif
 
 .include <bsd.lib.mk>
diff --git a/lib/msun/ld128/b_tgammal.c b/lib/msun/ld128/b_tgammal.c
index 4bae4f3aded6..6df7264a4c9e 100644
--- a/lib/msun/ld128/b_tgammal.c
+++ b/lib/msun/ld128/b_tgammal.c
@@ -1,55 +1,12 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
  *
- * Copyright (c) 2013 David Chisnall
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
+ * Copyright (c) 2024 The FreeBSD Foundation
  */
 
-#include <float.h>
-#include <math.h>
-
 /*
- * If long double is not the same size as double, then these will lose
- * precision and we should emit a warning whenever something links against
- * them.
+ * This is a pure C function generously donated by ARM.
+ * See src/contrib/arm-optimized-routines/math/tgamma128.[ch].
  */
-#if (LDBL_MANT_DIG > 53)
-#define WARN_IMPRECISE(x) \
-	__warn_references(x, # x " has lower than advertised precision");
-#else
-#define WARN_IMPRECISE(x)
-#endif
-/*
- * Declare the functions as weak variants so that other libraries providing
- * real versions can override them.
- */
-#define	DECLARE_WEAK(x)\
-	__weak_reference(imprecise_## x, x);\
-	WARN_IMPRECISE(x)
-
-#define DECLARE_IMPRECISE(f) \
-	long double imprecise_ ## f ## l(long double v) { return f(v); }\
-	DECLARE_WEAK(f ## l)
-
-DECLARE_IMPRECISE(tgamma);
+#define tgamma128 tgammal
+#include "tgamma128.c"
diff --git a/lib/msun/man/lgamma.3 b/lib/msun/man/lgamma.3
index cf80c9e7f365..41d680e41bda 100644
--- a/lib/msun/man/lgamma.3
+++ b/lib/msun/man/lgamma.3
@@ -1,216 +1,213 @@
 .\" Copyright (c) 1985, 1991 Regents of the University of California.
 .\" All rights reserved.
 .\"
 .\" Redistribution and use in source and binary forms, with or without
 .\" modification, are permitted provided that the following conditions
 .\" are met:
 .\" 1. Redistributions of source code must retain the above copyright
 .\"    notice, this list of conditions and the following disclaimer.
 .\" 2. Redistributions in binary form must reproduce the above copyright
 .\"    notice, this list of conditions and the following disclaimer in the
 .\"    documentation and/or other materials provided with the distribution.
 .\" 3. Neither the name of the University nor the names of its contributors
 .\"    may be used to endorse or promote products derived from this software
 .\"    without specific prior written permission.
 .\"
 .\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 .\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 .\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 .\" ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 .\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 .\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 .\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 .\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 .\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 .\" SUCH DAMAGE.
 .\"
-.\"     from: @(#)lgamma.3	6.6 (Berkeley) 12/3/92
-.\"
-.Dd December 8, 2017
+.Dd April 12, 2024
 .Dt LGAMMA 3
 .Os
 .Sh NAME
 .Nm lgamma ,
 .Nm lgamma_r ,
 .Nm lgammaf ,
 .Nm lgammaf_r ,
 .Nm lgammal ,
 .Nm lgammal_r ,
 .Nm gamma ,
 .Nm gamma_r ,
 .Nm gammaf ,
 .Nm gammaf_r ,
 .Nm tgamma ,
 .Nm tgammaf ,
 .Nm tgammal ,
 .Nd log gamma functions, gamma function
 .Sh LIBRARY
 .Lb libm
 .Sh SYNOPSIS
 .In math.h
 .Ft extern int
 .Fa signgam ;
 .sp
 .Ft double
 .Fn lgamma "double x"
 .Ft double
 .Fn lgamma_r "double x" "int *signgamp"
 .Ft float
 .Fn lgammaf "float x"
 .Ft float
 .Fn lgammaf_r "float x" "int *signgamp"
 .Ft "long double"
 .Fn lgammal "long double x"
 .Ft "long double"
 .Fn lgammal_r "long double x" "int *signgamp"
 .Ft double
 .Fn gamma "double x"
 .Ft double
 .Fn gamma_r "double x" "int *signgamp"
 .Ft float
 .Fn gammaf "float x"
 .Ft float
 .Fn gammaf_r "float x" "int *signgamp"
 .Ft "long double"
 .Fn tgamma "double x"
 .Ft float
 .Fn tgammaf "float x"
 .Ft "long double"
 .Fn tgammal "long double x"
 .Sh DESCRIPTION
 .Fn lgamma x ,
 .Fn lgammaf x ,
 and
 .Fn lgammal x
 .if t \{\
 return ln\||\(*G(x)| where
 .Bd -unfilled -offset indent
 \(*G(x) = \(is\d\s8\z0\s10\u\u\s8\(if\s10\d t\u\s8x\-1\s10\d e\u\s8\-t\s10\d dt	for x > 0 and
 \(*G(x) = \(*p/(\(*G(1\-x)\|sin(\(*px))	for x < 1.
 .Ed
 .\}
 .if n \
 return ln\||\(*G(x)|.
 The external integer
 .Fa signgam
 returns the sign of \(*G(x).
 .Pp
 .Fn lgamma_r x signgamp ,
 .Fn lgammaf_r x signgamp ,
 and
 .Fn lgammal_r x signgamp
 provide the same functionality as
 .Fn lgamma x ,
 .Fn lgammaf x ,
 and
 .Fn lgammal x ,
 but the caller must provide an integer to store the sign of \(*G(x).
 .Pp
 The
 .Fn tgamma x ,
 .Fn tgammaf x ,
 and
 .Fn tgammal x
 functions return \(*G(x), with no effect on
 .Fa signgam .
 .Pp
 .Fn gamma ,
 .Fn gammaf ,
 .Fn gamma_r ,
 and
 .Fn gammaf_r
 are deprecated aliases for
 .Fn lgamma ,
 .Fn lgammaf ,
 .Fn lgamma_r ,
 and
 .Fn lgammaf_r ,
 respectively.
 .Sh IDIOSYNCRASIES
 Do not use the expression
 .Dq Li signgam\(**exp(lgamma(x))
 to compute g := \(*G(x).
 Instead use a program like this (in C):
 .Bd -literal -offset indent
 lg = lgamma(x); g = signgam\(**exp(lg);
 .Ed
 .Pp
 Only after
 .Fn lgamma
 or
 .Fn lgammaf
 has returned can signgam be correct.
 .Pp
 For arguments in its range,
 .Fn tgamma
 is preferred, as for positive arguments
 it is accurate to within one unit in the last place.
 Exponentiation of
 .Fn lgamma
 will lose up to 10 significant bits.
 .Sh RETURN VALUES
 .Fn gamma ,
 .Fn gammaf ,
 .Fn gammal ,
 .Fn gamma_r ,
 .Fn gammaf_r ,
 .Fn gammal_r ,
 .Fn lgamma ,
 .Fn lgammaf ,
 .Fn lgammal ,
 .Fn lgamma_r ,
 .Fn lgammaf_r ,
 and
 .Fn lgammal_r
 return appropriate values unless an argument is out of range.
 Overflow will occur for sufficiently large positive values, and
 non-positive integers.
 For large non-integer negative values,
 .Fn tgamma
 will underflow.
-.Sh BUGS
-To conform with newer C/C++ standards, a stub implementation for
-.Nm tgammal
-was committed to the math library, where
-.Nm tgammal
-is mapped to
-.Nm tgamma .
-Thus, the numerical accuracy is at most that of the 53-bit double
-precision implementation.
 .Sh SEE ALSO
 .Xr math 3
 .Sh STANDARDS
 The
 .Fn lgamma ,
 .Fn lgammaf ,
 .Fn lgammal ,
 .Fn tgamma ,
 .Fn tgammaf ,
 and
 .Fn tgammal
 functions are expected to conform to
 .St -isoC-99 .
 .Sh HISTORY
 The
 .Fn lgamma
 function appeared in
 .Bx 4.3 .
 The
 .Fn gamma
 function appeared in
 .Bx 4.4
 as a function which computed \(*G(x).
 This version was used in
 .Fx 1.1 .
 The name
 .Fn gamma
 was originally dedicated to the
 .Fn lgamma
 function,
 and that usage was restored by switching to Sun's fdlibm in
 .Fx 1.1.5 .
 The
 .Fn tgamma
 function appeared in
 .Fx 5.0 .
+The 128-bit
+.Ft long double
+version of
+.Fn tgammal
+replaced the 80-bit stub version in
+version in
+.Fx 14.1 ,
+thanks to an appropriate implementation from Arm.