lgamma(x) = -log(x) - log(1+x) + x*(1-g) + x**2*P(x) with g = 0.57...
being the Euler constant and P(x) a polynomial. Substitution of small
into the RHS shows that the last 3 terms are negligible in comparison to
the leading term. The choice of 3 may be conservative.
The value large=2**(p+3) is detemined from Stirling's approximation
lgamma(x) = x*(log(x)-1) - log(x)/2 + log(2*pi)/2 + P(1/x)/x
Again, substitution of large into the RHS reveals the last 3 terms
are negligible in comparison to the leading term.
Move the x=+-0 special case into the |x|<small block.
In the ld80 and ld128 implementaion, use fdlibm compatible comparisons
involving ix, lx, and llx. This replaces several floating point
comparisons (some involving fabsl()) and also fixes the special cases
x=1 and x=2.
While here
. Remove unnecessary parentheses.
. Fix/improve comments due to the above changes.
. Fix nearby whitespace.
* src/e_lgamma_r.c:
. Sort declaration.
. Remove unneeded explicit cast for type conversion.
. Replace a double literal constant by an integer literal constant.
* src/e_lgammaf_r.c:
. Sort declaration.
* ld128/e_lgammal_r.c:
. Replace a long double literal constant by a double literal constant.
* ld80/e_lgammal_r.c:
. Remove unused '#include float.h'
. Replace a long double literal constant by a double literal constant.
Requested by: bde
. Hook e_lgammal[_r].c to the build.
. Create man page links for lgammal[-r].3.
* Symbol.map:
. Sort lgammal to its rightful place.
. Add FBSD_1.4 section for the new lgamal_r symbol.
* ld128/e_lgammal_r.c:
. 128-bit implementataion of lgammal_r().
* ld80/e_lgammal_r.c:
. Intel 80-bit format implementation of lgammal_r().
* src/e_lgamma.c:
. Expose lgammal as a weak reference to lgamma for platforms
where long double is mapped to double.
* src/e_lgamma_r.c:
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
. Expose lgammal_r as a weak reference to lgamma_r for platforms
where long double is mapped to double.
* src/e_lgammaf_r.c:
. Fixed the Cygnus Support conversion of e_lgamma_r.c to float.
This includes the generation of new polynomial and rational
approximations with fewer terms. For each approximation, include
a comment on an estimate of the accuracy over the relevant domain.
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
This allows the removal of several explicit casts of double values
to float.
* src/e_lgammal.c:
. Wrapper for lgammal() about lgammal_r().
* src/imprecise.c:
. Remove the lgamma.
* src/math.h:
. Add a prototype for lgammal_r().
* man/lgamma.3:
. Document the new functions.
Reviewed by: bde
inf and raises the divided-by-zero exception. Compilers
constant fold one/zero to inf but do not raise the exception.
Introduce a volatile vzero to prevent the constant folding.
Move the declaration of zero into the main declaration block.
While here, fix a nearby disordering of 'lx,ix'
Discussed with: bde
Remove parentheses in a return statement to be consistent with the
rest of the file.
Rename sin_pi() in the float version to sin_pif().
Remove large comment that precedes sin_pif(). The comment
duplicates a comment in e_lgamma_r.c where the algorithm
is documented.
Requested by: bde
sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
the precision of x. The new argument reduction is an optimization
compared to the old code, and it removes a chunk of dead code.
Accuracy tests in the intervals (-21,-20), (-20,-19), ... (-1,0)
show no differences between the old and new code.
Obtained from: bde
By Richard Earnshaw at ARM
>
>GCC has for a number of years provides a set of pre-defined macros for
>use with determining the ISA and features of the target during
>pre-processing. However, the design was always somewhat cumbersome in
>that each new architecture revision created a new define and then
>removed the previous one. This meant that it was necessary to keep
>updating the support code simply to recognise a new architecture being
>added.
>
>The ACLE specification (ARM C Language Extentions)
>(http://infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.set.swdev/index.html)
>provides a much more suitable interface and GCC has supported this
>since gcc-4.8.
>
>This patch makes use of the ACLE pre-defines to map to the internal
>feature definitions. To support older versions of GCC a compatibility
>header is provided that maps the traditional pre-defines onto the new
>ACLE ones.
Stop using __FreeBSD_ARCH_armv6__ and switch to __ARM_ARCH >= 6 in the
couple of places in tree. clang already implements ACLE. Add a define
that says we implement version 1.1, even though the implementation
isn't quite complete.
and tgammal in libm. These functions are part of ISO/IEC 9899:1999
and their prototypes should have been moved into the appropriate
__ISO_C_VISIBLE >= 1999 section. After moving the prototypes,
remnants of r236148 can be removed.
PR: standards/191754
Reviewed by: bde
. Add s_erfl.c to building libm.
. Add MLINKS for erfl.3 and erfcl.3.
* Symbol.map:
. Move erfl and erfcl to their proper location.
* ld128/s_erfl.c:
. Implementations of erfl and erfcl in the IEEE 754 128-bit format.
* ld80/s_erfl.c:
. Implementations of erfl and erfcl in the Intel 80-bit format.
* man/erf.3:
. Document the new functions.
. While here, remove an incomplete sentence.
* src/imprecise.c:
. Remove the stupidity of mapping erfl and erfcl to erf and erfc.
* src/math.h:
. Move the declarations of erfl and erfcl to their proper place.
* src/s_erf.c:
. For architectures where double and long double are the same
floating point format, use weak references to map erfl to
erf and ercl to erfc.
Reviewed by: bde (many earlier versions)
* Update the domain and range of comments for the polynomial
approximations, including using the the correct variable names
(e.g., pp(x) instead of p(x)).
* Use hex values of the form 0x3e0375d4 instead of 0x1.06eba8p-3,
which was obtained from printf("%.6a").
* In the domain [0.84375, 1.25], qa(x) can be reduced from a 4th
order polynomial to 3rd order.
* In the domain [1.25,1/0.35], sa(x) can be reduced from a 4th
order polynomial to 3rd order.
* In the domain [1/0.35, 11], the 4th order polynomials rb(x) and
sb(x) can be reduced to 2nd and 3rd order, respectively.
from erronously constant folding expressions of the form
'1 - tiny'. This allows erf[f](x) to raise INEXACT.
* Use 0.5, 1, and 2, which are exactly representable in radix-2
floating point formats. This reduces diffs between s_erf[fl].c.
* While here, add a comment about efx and efx8.
This includes:
o All directories named *ia64*
o All files named *ia64*
o All ia64-specific code guarded by __ia64__
o All ia64-specific makefile logic
o Mention of ia64 in comments and documentation
This excludes:
o Everything under contrib/
o Everything under crypto/
o sys/xen/interface
o sys/sys/elf_common.h
Discussed at: BSDcan
* msun/man/sinh.3:
* msun/man/tanh.3:
. Fix grammar.
* msun/src/e_coshl.c:
* msun/src/e_sinhl.c:
. Fix comment.
* msun/src/s_tanhl.c:
. Remove unused variables.
. Fix location/indentation of comments.
. Use comparison involving ints instead of long double.
. Re-order polynomial evaluation on ld128 for |x| < 0.25.
For now, retain the older order in an "#if 0 ... #else" block.
. Use int comparison to short-circuit the |x| < 1.5 condition.
Requested by: bde
. Hook coshl, sinhl, and tanhl into libm.
. Create symbolic links for corresponding manpages.
. While here remove a nearby extraneous space.
* Symbol.map:
* src/math.h:
. Move coshl, sinhl, and tanhl to their proper locations.
* man/cosh.3:
* man/sinh.3:
* man/tanh.3:
. Update the manpages.
* src/e_cosh.c:
* src/e_sinh.c:
* src/s_tanh.c:
. Add weak reference for LBDL_MANT_DIG==53 targets.
* src/imprecise.c:
. Remove the coshl, sinhl, and tanhl kludge.
* src/e_coshl.c:
. ld80 and ld128 implementation of coshl().
* src/e_sinhl.c:
. ld80 and ld128 implementation of sinhl().
* src/s_tanhl.c:
. ld80 and ld128 implementation of tanhl().
Obtained from: bde (mostly), das and kargl
* ld128/k_expl.h:
. Split out a computational kernel,__k_expl(x, &hi, &lo, &k) from expl(x).
x must be finite and not tiny or huge. The kernel returns hi and lo
values for extra precision and an exponent k for a 2**k scale factor.
. Define additional kernels k_hexpl() and hexpl() that include a 1/2
scaling and are used by the hyperbolic functions.
* ld80/s_expl.c:
* ld128/s_expl.c:
. Use the __k_expl() kernel.
Obtained from: bde
* Use bit twiddling. This requires inclusion of math_private.h
and inclusion of float.h in s_roundl.c. Raise invalid exception.
* Use literal integer constants where possible. Let the compiler
do the appropriate conversion.
* In s_roundf.c, use an F suffix on float constants instead of
promoting float to double and then converting the result back
to float. In s_roundl.c, use an L suffix.
* In s_roundl.c, use the ENTERI and RETURNI macros. This requires
the inclusion of fpmath.h and on __i386__ class hardware ieeefp.h.
Reviewed by: bde
preprocessor) gives the following error:
--- Version.map ---
<stdin>:287:4: error: invalid preprocessing directive
# Implemented as weak aliases for imprecise versions
^
1 error generated.
Change the comment to a C-style one, to prevent this error.
Approved by: re (hrs)
These are weak and so can be replaced by other versions in applications
that choose to do so, and will give a linker warning when used so that
applications that rely on the extra precision can avoid them.
Note that since the C/C++ specs only guarantee that long double has
precision equal to double, code that actually relies on these functions
having greater precision is unportable at best and broken at worst.
. Use integer literal constants instead of double literal constants.
* s_erff.c:
. Use integer literal constants instead of casting double literal
constants to float.
. Update the threshold values from those carried over from erf() to
values appropriate for float.
. New sets of polynomial coefficients for the rational approximations.
These coefficients have little, but positive, effect on the maximum
error in ULP in the four intervals, but do improve the overall
speed of execution.
. Remove redundant GET_FLOAT_WORD(ix,x) as hx already contained the
contents that is packed into ix.
. Update the mask that is used to zero-out lower-order bits in x in
the intervals [1.25, 2.857143] and [2.857143, 12]. In tests on
amd64, this change improves the maximum error in ULP from 6.27739
and 63.8095 to 3.16774 and 2.92095 on these intervals for erffc().
Reviewed by: bde
