The principal expressions defining acosh and acos are such that
acosh(z) = ±i acos(z)
where the + is only true on the Im(z)>0 half of the complex plane
(and partly on Im(z)==0 depending on number representation).
fix the comment without expanding on the details.
gcc-12 with -frounding-mode will do inexact constant conversions at
runtime according to the runtime rounding mode.
in the math library we want constants to be rounding mode independent
so this patch fixes cases where new runtime conversions happen with
gcc-12.
fortunately this only affects two minor cases, the fix uses global
initializers where rounding mode does not apply.
after the patch the same amount of conversions happen with gcc-12 as
with gcc-11.
these are specified to use the sign of the imaginary part of the input
as the sign of zero in the result, but wrongly copied the sign of the
real part.
This makes it easier to build musl math code with a compiler that
does not support complex types (tcc) and in general more sensible
factorization of the internal headers.
the catan implementation from OpenBSD includes a FIXME-annotated
"overflow" branch that produces a meaningless and incorrect
large-magnitude result. it was reachable via three paths,
corresponding to gotos removed by this commit, in order:
1. pure imaginary argument with imaginary component greater than 1 in
magnitude. this case does not seem at all exceptional and is
handled (at least with the quality currently expected from our
complex math functions) by the existing non-exceptional code path.
2. arguments on the unit circle, including the pure-real argument 1.0.
these are not exceptional except for ±i, which should produce
results with infinite imaginary component and which lead to
computation of atan2(±0,0) in the existing non-exceptional code
path. such calls to atan2() however are well-defined by POSIX.
3. the specific argument +i. this route should be unreachable due to
the above (2), but subtle rounding effects might have made it
possible in rare cases. continuing on the non-exceptional code path
in this case would lead to computing the (real) log of an infinite
argument, then producing a NAN when multiplying it by I.
for now, remove the exceptional code paths entirely. replace the
multiplication by I with construction of a complex number using the
CMPLX macro so that the NAN issue (3) prevented cannot arise.
with these changes, catan should give reasonably correct results for
real arguments, and should no longer give completely-wrong results for
pure-imaginary arguments outside the interval (-i,+i).
the factor of -i noted in the comment at the top of casin.c was
omitted from the actual code, yielding a result rotated 90 degrees and
propagating into errors in other functions defined in terms of casin.
implement multiplication by -i as a rotation of the real and imaginary
parts of the result, rather than by actual multiplication, since the
latter cannot be optimized without knowledge that the operand is
finite. here, the rotation is the actual intent, anyway.
thanks to the hard work of Szabolcs Nagy (nsz), identifying the best
(from correctness and license standpoint) implementations from freebsd
and openbsd and cleaning them up! musl should now fully support c99
float and long double math functions, and has near-complete complex
math support. tgmath should also work (fully on gcc-compatible
compilers, and mostly on any c99 compiler).
based largely on commit 0376d44a890fea261506f1fc63833e7a686dca19 from
nsz's libm git repo, with some additions (dummy versions of a few
missing long double complex functions, etc.) by me.
various cleanups still need to be made, including re-adding (if
they're correct) some asm functions that were dropped.
Rich Felker
Szabolcs Nagy
Michael Morrell