Let's get moving

[Bruce Evans]

On Mon, 1 Jun 2015, Montgomery-Smith, Stephen wrote:

> Hey people,
>
> I was looking at https://wiki.freebsd.org/Numerics, and we are stalling.
> For example, we have working code for clog and clogf - why don't we
> commit it?

I thought that you might commit these before your complex inverse
hyperbolic functions :-).

catrig.c and catrigf.c have been committed.  These pass my accuracy
tests, but I'm not happy with their integration.  Bugs start in their
file names.  All other files in msun/src have names with prefixes like
e_ or s_.  Comments are stripped in the committed version of catrigf.c,
and your template for maintaining the comments was not committed.

catrigl.c was not committed.  I only have incomplete minor patches 
for this.  Most were in patches that I sent you in 2012 to merge:

X --- stephen-catrigl.c	2012-09-22 21:14:24.000000000 +0000
X +++ catrigl.c	2013-06-08 10:47:56.000000000 +0000
X @@ -60,6 +69,6 @@
X  #if LDBL_MANT_DIG == 64
X  static const union IEEEl2bits
X -um_e =		LD80C(0xadf85458a2bb4a9b,  1, 0, 2.71828182845904523536e0L),
X -um_ln2 =	LD80C(0xb17217f7d1cf79ac, -1, 0, 6.93147180559945309417e-1L);
X +um_e =		LD80C(0xadf85458a2bb4a9b,  1, 2.71828182845904523536e+0L),
X +um_ln2 =	LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
X  #define		m_e	um_e.e
X  #define		m_ln2	um_ln2.e

Catch up with API changes in LD80C().  These have been in -current since
2012 so are standard now.

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

In 2012, this used our hacked up versions of atanhl().  A production
version was committed in 2013.

Someone churned the API in math_private.h so that cpackl() no longer
exists.  I don't used the churned parts, so the above still compiles
for me.

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

Minor cleanups.  IIRC, this was not in the template, and the float,
double and long versions had difference here from being manually
edited.  The committed float and double versions have this and some
other cleanups, especially in the float version where generation
from the template made a mess.

The main things missing here are:
- no support for i387 mode switching (ENTERI()).  Not very easy to add
   cleanly since there are many entry and exit points.
- not up to date with cpackl().

My version of clog*() is a little more finished.  It has the i387 mode
switching support for clogl(), but is not up to date with cpack*().  It
needs splitting into separate files for each precision.  It handles
comment duplication not very well by just duplicating comments and
manually maintaining consistency.

catrig*() can also use clog*() when it exists, but have fairly good
private versions for the special cases needed.

> One of the reasons I bring this up is I would like to create a port for
> http://librsb.sourceforge.net/, and the only barriers are the lack of
> cpowf and cpow.  Right now I intend to commit it with a rather
> simplistic patch:
>
> +complex float clogf(complex float a) {
> +       return logf(cabsf(a)) + I*cargf(a);
> +}
> +
> +complex float cpowf(complex float a, complex float b) {
> +       return cexpf(b*clogf(a));
> +}
> +
> +complex double clog(complex double a) {
> +       return log(cabs(a)) + I*carg(a);
> +}
> +
> +complex double cpow(complex double a, complex double b) {
> +       return cexp(b*clog(a));
> +}

I think a lot of ports already have similar hacks.

libm has similar hacks for powl() and tgammal(), to support C++.  It
doesn't bother with similar hacks for complex functions.

C99 and at least old POSIX don't have complex Bessel or gamma functions,
so we're closer to having all standard real functions that complex ones.
Only powl() and cpowl() are standard, unsupported, not written, and
difficult.

Bruce