complex arc-trig etc

Bruce Evans brde at optusnet.com.au
Mon Dec 17 08:39:01 UTC 2012


On Sun, 16 Dec 2012, Stephen Montgomery-Smith wrote:

> Hey guys, my complex arc-trig functions at
> http://people.freebsd.org/~stephen/ have been sitting there a long time.
> Anyone want to commit them?

I was sort of waiting for the next rounds of changes:

1. update the template to give similar code for float and long double
2. make it work for i386 using ENTERI()
3. test (2)
4. fix some style bugs.  It would be easy for me to run it through
    indent(1) and get large diffs, but I want you to do it.
5. decide what to do for atanhl().  I would prefer a version that is a
    direct translation of e_atanh.c including style bugs and the filename
    prefix but not the comment, unless you change the algorithm.
6. decide what to do with the template and/or the generation from it
    that strips comments.

I haven't done (2), (3) or (6) for clog*() either.  For (6), instead of
a template and stripped comments, I have all the functions in the same
file and all comments duplicated, and have to split the file manually
to compare the functions.

I only have these local patches.  They correspond to (1) and (5).  Also,
keep up with the API change in LD80C().

% --- /home/stephen/public_html/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;
% 
% --- /home/stephen/public_html/catrigl.c	2012-09-22 21:14:24.000000000 +0000
% +++ catrigl.c	2012-09-28 19:39:47.799585000 +0000
% @@ -60,6 +60,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
% @@ -348,5 +348,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)
% @@ -369,12 +369,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;
%

You also need log*l() be committed before catrigl() can work.

Bruce


More information about the freebsd-numerics mailing list