Use of C99 extra long double math functions after r236148

Bruce Evans brde at optusnet.com.au
Tue May 28 05:57:36 UTC 2013


On Mon, 27 May 2013, David Schultz wrote:

> I wrote some tests to cover the corner cases for the complex
> inverse trig functions. They don't find any nontrivial bugs in
> your implementations. :-) Now that you have a commit bit, would
> you like to commit your code, or shall I?
>
> Below is a diff of all the changes needed to integrate it. I have
> a short list of style fixes, but otherwise I think what you have
> is good:
>  - wrap lines to 80 chars, please
>  - spaces between operators
>  - "static inline", not "inline static"
>  - don't use "inline" on large functions

indent(1) fixes the spaces between operators fairly well, without finding
many other problems or adding many.  It didn't find [m]any long lines
(but it doesn't understand its own line length limit).

Here are my local patches.  Just a few that were not integrated by
Stephen after we stopped working on it last October.

@ diff -u2 catrigf.c~ catrigf.c
@ --- catrigf.c~	2012-09-22 21:13:50.000000000 +0000
@ +++ catrigf.c	2012-09-22 21:35:51.287614000 +0000
@ @@ -353,12 +353,7 @@
@  	}
@ 
@ -	if (ax == 1 && ay < FLT_EPSILON) {
@ -#if 0
@ -		if (ay > 2*FLT_MIN)
@ -			rx = - logf(ay/2) / 2;
@ -		else
@ -#endif
@ -			rx = - (logf(ay) - m_ln2) / 2;
@ -	} else
@ +	if (ax == 1 && ay < FLT_EPSILON)
@ +		rx = - (logf(ay) - m_ln2) / 2;
@ +	else
@  		rx = log1pf(4*ax / sum_squares(ax-1, ay)) / 4;
@

This is in catrig.c, but catrigf.c wasn't regenerated from catrig.c, and
the scripts for the generation and their support file are no longer in
stephen's public_html directory.

@ diff -u2 catrigl.c~ catrigl.c
@ --- catrigl.c~	2012-09-22 21:14:24.000000000 +0000
@ +++ catrigl.c	2013-05-26 08:46:10.423187000 +0000
@ @@ -50,4 +50,6 @@
@  #define signbit(x)	(__builtin_signbitl(x)) 
@ 
@ +long double atanhl(long double);
@ +
@  static const long double
@  A_crossover =		10,

catrigl.c depends on atanhl(), logl() and log1pl() existing.  Stephen
has a not-very-dummy version of s_atanhl.c in this public_html
directory.  This needs a more direct conversion from the fdlibm
e_atanhl.c to be of commit quality.  I recently started testing with
it, and use my own logl().  Previously this patch had to change the
atanhl() call to atanh() to for catrigl.c to be usable.  I haven't
tested the long double complex functions for anything except efficiency
and consistency with the plain double complex functions yet, so my
tests don't should any difference from switching to atanhl().  They
just show that atanhl() is consistent in its limited use in catrigl.c.
I also haven't tested atanhl() as a real function.

Strangely, catrigl.c gives complex acoshl() and asinhl() without needing
real acoshl() and asinhl().  The real inverse hyperbolic trig functions
seem to be just as easy as the real inverse trig functions, but you
only converted the latter from the fdlibm versions to create the long
double versions.  Hopefully they are all as easy to translate
e_atanhl.c.

@ @@ -60,6 +62,6 @@
@  #if LDBL_MANT_DIG == 64
@  static const union IEEEl2bits
@ -um_e =		LD80C(0xadf85458a2bb4a9b,  1, 0, 2.71828182845904523536e0L),
@ -um_ln2 =	LD80C(0xb17217f7d1cf79ac, -1, 0, 6.93147180559945309417e-1L);
@ +um_e =		LD80C(0xadf85458a2bb4a9b,  1, 2.71828182845904523536e+0L),
@ +um_ln2 =	LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
@  #define		m_e	um_e.e
@  #define		m_ln2	um_ln2.e

Keep up with API changes.

@ @@ -348,5 +350,5 @@
@ 
@  	if (y == 0 && ax <= 1)
@ -		return (cpackl(atanhl(x), y)); 	/* XXX need atanhl() */
@ +		return (cpackl(atanh(x), y)); 	/* XXX need atanhl() */
@ 
@  	if (x == 0)

The comment doesn't apply if this file is actually usable.  Don't forget
to remove it before committing.

@ @@ -369,12 +371,7 @@
@  	}
@ 
@ -	if (ax == 1 && ay < LDBL_EPSILON) {
@ -#if 0
@ -		if (ay > 2*LDBL_MIN)
@ -			rx = - logl(ay/2) / 2;
@ -		else
@ -#endif
@ -			rx = - (logl(ay) - m_ln2) / 2;
@ -	} else
@ +	if (ax == 1 && ay < LDBL_EPSILON)
@ +		rx = - (logl(ay) - m_ln2) / 2;
@ +	else
@  		rx = log1pl(4*ax / sum_squares(ax-1, ay)) / 4;
@

Should be obtained by regeneration, as for catrigf.c.

Back to your changes...  They mostly look good, as usual...

% Index: tools/regression/lib/msun/test-invctrig.c
% ===================================================================
% --- tools/regression/lib/msun/test-invctrig.c	(revision 0)
% +++ tools/regression/lib/msun/test-invctrig.c	(working copy)
% @@ -0,0 +1,467 @@
% ....
% +#pragma STDC FENV_ACCESS	ON
% +#pragma	STDC CX_LIMITED_RANGE	OFF

Heheh, style rules for #pragma.  I like the old rule which says that
it should be indented 6 feet under.  It is still almost useless, since
we don't even have any C99 compilers than implement the fenv pragmas
yet.

% +/*
% + * XXX gcc implements complex multiplication incorrectly. In
% + * particular, it implements it as if the CX_LIMITED_RANGE pragma
% + * were ON. Consequently, we need this function to form numbers
% + * such as x + INFINITY * I, since gcc evalutes INFINITY * I as
% + * NaN + INFINITY * I.
% + */
% +static inline long double complex
% +cpackl(long double x, long double y)
% +{
% +	long double complex z;
% +
% +	__real__ z = x;
% +	__imag__ z = y;
% +	return (z);
% +}

Why duplicate this?  I guess it is because math_private,h is hard to
include.  I use complicated conditionals (mostly switches on
$(uname -p) and $(hostname) in shell scripts to locate it when
compiling from external directories.

The tests seem to be compiled with -O0.  That tests a different
environment than the usual runtime one, and in particular misses seeing
most precision bugs.  I mostly test with -O (-O2 with gcc is slower
and even harder to debug, while with clang it makes little difference),
but switch to -O0 to debug.  -g -O is now almost unusable because -O
optimizes away dead variables and -g is broken in many cases (sometimes
it can't even show live variables).

Bruce


More information about the freebsd-numerics mailing list