libm

libm

[Szabolcs Nagy]

i've looked into libm implementations
to figure out what's best for musl

some resources:
    http://nsz.repo.hu/libm.txt

there are various issues:

* reproducible semantics:

linux uses double extended precision by default on i386
which does not give the same results as standard 64bit
ieee754 (sse2, x86_64, arm)

bsd uses plain double precision by default on i386
(so at least when bsd functions are ported the precision
claims should be reevaluated)

most of the recent advances in efficient and correct
function evaluation assume strict ieee754 64bit arithmetics
(so they can do formal proofs etc, some algorithms assume
fma instruction which is not available on musl's targets)

portable semantics are not easy to guarantee (extended
precision, rounding modes, compiler optimizations)
(especially if efficiency is important as well)

* tradeoffs:

modern libms (libmcr, libultim, crlibm) try to guarantee
correctness but arg reduction and faithful (or even
correct) rounding is hard to do and hard to verify

some libms prefer simple implementation (eg the go math
package sincos is from cephes and starts to get inaccurate
when |x| > 2^30, but the implementation is much simpler
than the arg reduction in fdlibm)

in theory one should prefer correctness instead of size
and speed, but there are significant differences in code
size and implementation effort
the speed of a correct implementation can be good in
average case, but not in worstcase

* libms in practice:

many functions are the same in glibc and the various
bsd libcs (they are mostly based on fdlibm, but glibc
64bit double precision functions are from the gpl
licensed libultim)

the extended precision algorithms are reused across
different libcs as well, but there are 80bit vs 128bit
differences. the correctness of extended precision
algorithms are much less studied (eg there are no
correctly rounded versions or worstcases for 2pi arg
reductions)

most of the complex functions are simple once elementary
functions are there (bsd libcs has bsd licensed
implementations of these)


conclusion:

the simplest approach for musl at this point is to
reuse the already available math functions
(+run the available tests on them)

this can be done by diffing the various bsd and glibc
functions and choosing the best one (usually they are
the same except code organization and bit manipulation
details)

most of the task is understanding floating point
semantics well (not the mathematics of the algorithms)
and doing code cleanups

code and ideas from crlibm might be possible to use
but that would take much more effort and knowledge
(assuming we want small code size)

designing new algorithms for elementary functions seems
to require a huge effort and the best tradoffs are not
even clear

Re: libm

[Rich Felker]

On Mon, Jan 23, 2012 at 05:41:52PM +0100, Szabolcs Nagy wrote:
> i've looked into libm implementations
> to figure out what's best for musl

Thanks!!

> * reproducible semantics:

Unless others are deeply committed to having results match other
platforms/libcs bit-for-bit, this is not a major concern for me.
I'm happy with bit-exact just for non-transcendental (i.e. sqrt) and
1ulp (or even several ulp) max error for transcendental functions.

> * tradeoffs:
> 
> modern libms (libmcr, libultim, crlibm) try to guarantee
> correctness but arg reduction and faithful (or even
> correct) rounding is hard to do and hard to verify

Correct arg reduction is important to me. Incorrect arg reduction is
equivalent to billions of billions of ulp error in the final results.

> * libms in practice:
> 
> many functions are the same in glibc and the various
> bsd libcs (they are mostly based on fdlibm, but glibc
> 64bit double precision functions are from the gpl
> licensed libultim)

You really mean double, not extended? That's odd since fdlibm covers
double.. I wonder when/why they switched.

> the extended precision algorithms are reused across
> different libcs as well, but there are 80bit vs 128bit
> differences. the correctness of extended precision
> algorithms are much less studied (eg there are no
> correctly rounded versions or worstcases for 2pi arg
> reductions)

Any ideas how the different ones evolved (separately written or common
ancestor code, etc.) and where we should look to pull code from?

> most of the complex functions are simple once elementary
> functions are there (bsd libcs has bsd licensed
> implementations of these)

Well they don't have any more difficult numeric analysis issues, but
there are quite a few corner cases with complex functions which many
platforms have gotten wrong (issues with branch cuts and signed zeros,
etc.).

> conclusion:
> 
> the simplest approach for musl at this point is to
> reuse the already available math functions
> (+run the available tests on them)

I generally agree, provided no major bugs are found - and even if they
are, they should be fixable.

> this can be done by diffing the various bsd and glibc
> functions and choosing the best one (usually they are
> the same except code organization and bit manipulation
> details)

That was my impression too.

> code and ideas from crlibm might be possible to use
> but that would take much more effort and knowledge
> (assuming we want small code size)

Yes, crlibm looks very interesting but also very large...

> designing new algorithms for elementary functions seems
> to require a huge effort and the best tradoffs are not
> even clear

Agreed. I think it's roughly as big a task as the rest of musl
combined, so it's best not to go there.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-01-23 12:07:15 -0500]:
> On Mon, Jan 23, 2012 at 05:41:52PM +0100, Szabolcs Nagy wrote:
> > * tradeoffs:
> > 
> > modern libms (libmcr, libultim, crlibm) try to guarantee
> > correctness but arg reduction and faithful (or even
> > correct) rounding is hard to do and hard to verify
> 
> Correct arg reduction is important to me. Incorrect arg reduction is
> equivalent to billions of billions of ulp error in the final results.
> 

yes, but trigonometric arg reduction is a real pain
and ppl shouldn't use large arguments there anyway..

so there are arguments both ways
but this is already reasonably
solved in fdlibm

> > * libms in practice:
> > 
> > many functions are the same in glibc and the various
> > bsd libcs (they are mostly based on fdlibm, but glibc
> > 64bit double precision functions are from the gpl
> > licensed libultim)
> 
> You really mean double, not extended? That's odd since fdlibm covers
> double.. I wonder when/why they switched.
> 

i forgot to mention that both the bsds and glibc have most of the
functions implemented in i387 asm (using the builtin instructions),
i'm not sure how far musl wants to take that

the generic 64bit double c code of glibc is from fdlibm and libultim
but libultim provides the important elementary functions (exp, log, pow,
sin, atan,..) while fdlibm the more obscure or simple ones (floor, rint,
fabs, cosh, tanh, j0, erf,..)

$ grep '* IBM' glibc/sysdeps/ieee754/dbl-64/*.c |wc -l
25
$ grep 'Copyright .* Sun' glibc/sysdeps/ieee754/dbl-64/*.c |wc -l
33

the libultim functions provide correct rounding but can be slow
(they use multiprecision arithmetic for various things and
eg for arg reduction they increase the precision until it's enough)

> > the extended precision algorithms are reused across
> > different libcs as well, but there are 80bit vs 128bit
> > differences. the correctness of extended precision
> > algorithms are much less studied (eg there are no
> > correctly rounded versions or worstcases for 2pi arg
> > reductions)
> 
> Any ideas how the different ones evolved (separately written or common
> ancestor code, etc.) and where we should look to pull code from?
> 

most of the available *l (and *f) code is based on fdlibm
but ported by someone to long double (or float)

openbsd and freebsd mostly use the same code by the same contributors
as far as i can see, but they have separate 80bit and 128bit variants
for important elementary functions

glibc is sometimes the same as the bsds but often not (eg expl and
sinl are clearly different)

i don't think there are other available widely used implementations
so we have openbsd, freebsd and glibc
(mathcw has long double functions but i don't see source code and
it does not seem to be active)

Re: libm

[Szabolcs Nagy]

* Szabolcs Nagy <nsz@port70.net> [2012-01-23 20:12:15 +0100]:
> * Rich Felker <dalias@aerifal.cx> [2012-01-23 12:07:15 -0500]:
> > On Mon, Jan 23, 2012 at 05:41:52PM +0100, Szabolcs Nagy wrote:
> > > many functions are the same in glibc and the various
> > > bsd libcs (they are mostly based on fdlibm, but glibc
> > > 64bit double precision functions are from the gpl
> > > licensed libultim)
> > 
> > You really mean double, not extended? That's odd since fdlibm covers
> > double.. I wonder when/why they switched.
> > 
> 

sorry libultim is lgpl
the source tree contains a gpl license for some reason, but all
source comments say lgpl, both in glibc and the original files

meanwhile i compared some implementations
openbsd, freebsd and glibc are worth to look at
(netbsd and bionic libms seem to be less maintained)

Re: libm

[Pascal Cuoq]

[-- Attachment #1: Type: text/plain, Size: 1496 bytes --]

On Fri, Jan 27, 2012 at 5:02 PM, Szabolcs Nagy <nsz@port70.net> wrote:

> meanwhile i compared some implementations
> openbsd, freebsd and glibc are worth to look at
>

Sorry to intrude in your conversation, but for reasons of
my own I was recently looking for an aesthetically pleasing
single-precision implementation of trigonometric functions,
and some of the single-precisions functions out there are
so ugly that I would seriously consider implementing
these in musl by a call to their respective double-precision
counterparts, despite the disadvantages. At least it's compact.

Consider:
http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cosf.c?rev=1.3;content-type=text%2Fplain

The code written at Sun and that many libms reuse must
have been double-precision only, and someone had to
make a single-precision version of them at some point.
Without the analysis that led to the original code, one
cannot blame Ian Lance Taylor for staying close
to the double-precision function:
http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cos.c?rev=1.4;content-type=text%2Fplain

But the degree of the polynomial in the single-precision version
is about twice what is necessary, and the coefficients
look like coefficients of the Taylor expansion rounded to
nearest float, when they should be Chebyshev coefficients.

Criticism is easy and art is difficult. I don't know how
to write a better single-precision function, I just know
that this one, for instance, cannot be good.

Pascal

[-- Attachment #2: Type: text/html, Size: 2316 bytes --]

Re: libm

[Rich Felker]

On Fri, Jan 27, 2012 at 08:01:18PM +0100, Pascal Cuoq wrote:
> On Fri, Jan 27, 2012 at 5:02 PM, Szabolcs Nagy <nsz@port70.net> wrote:
> 
> > meanwhile i compared some implementations
> > openbsd, freebsd and glibc are worth to look at
> >
> 
> Sorry to intrude in your conversation, but for reasons of
> my own I was recently looking for an aesthetically pleasing
> single-precision implementation of trigonometric functions,
> and some of the single-precisions functions out there are
> so ugly that I would seriously consider implementing
> these in musl by a call to their respective double-precision
> counterparts, despite the disadvantages. At least it's compact.

*Nod*, this is the way I'm leaning for single-precision anyway. The
only possible reason to reconsider is that some systems (i.e. some
ARMs) have slow software emulation of double precision but native
single precision, or so I'm told...

> Consider:
> http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cosf.c?rev=1.3;content-type=text%2Fplain
> 
> The code written at Sun and that many libms reuse must
> have been double-precision only, and someone had to
> make a single-precision version of them at some point.
> Without the analysis that led to the original code, one
> cannot blame Ian Lance Taylor for staying close
> to the double-precision function:
> http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cos.c?rev=1.4;content-type=text%2Fplain
> 
> But the degree of the polynomial in the single-precision version
> is about twice what is necessary, and the coefficients
> look like coefficients of the Taylor expansion rounded to
> nearest float, when they should be Chebyshev coefficients.

Agreed, this is ridiculous.

> Criticism is easy and art is difficult. I don't know how
> to write a better single-precision function, I just know
> that this one, for instance, cannot be good.

Well I think it's fair to judge when they could have just wrapped the
double precision versions. Granted this can result in 1ulp error, but
it's not a big deal and within the tolerances required by the
standard...and if their implementations of the single-precision
versions are as bad as you say, I can't imagine they're always
correctly rounded anyway..

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-01-27 14:34:24 -0500]:
> On Fri, Jan 27, 2012 at 08:01:18PM +0100, Pascal Cuoq wrote:
> > On Fri, Jan 27, 2012 at 5:02 PM, Szabolcs Nagy <nsz@port70.net> wrote:
> > 
> > > meanwhile i compared some implementations
> > > openbsd, freebsd and glibc are worth to look at
> > >
> > 
> > Sorry to intrude in your conversation, but for reasons of
> > my own I was recently looking for an aesthetically pleasing
> > single-precision implementation of trigonometric functions,
> > and some of the single-precisions functions out there are
> > so ugly that I would seriously consider implementing
> > these in musl by a call to their respective double-precision
> > counterparts, despite the disadvantages. At least it's compact.
> 
> *Nod*, this is the way I'm leaning for single-precision anyway. The
> only possible reason to reconsider is that some systems (i.e. some
> ARMs) have slow software emulation of double precision but native
> single precision, or so I'm told...
> 

yes, this is good for code sizeand maintainability

some code will expect faster execution when it uses
the float versions so some benchmarking will be needed
for arm to see the tradeoffs

> > Consider:
> > http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cosf.c?rev=1.3;content-type=text%2Fplain
> > 
> > The code written at Sun and that many libms reuse must
> > have been double-precision only, and someone had to
> > make a single-precision version of them at some point.
> > Without the analysis that led to the original code, one
> > cannot blame Ian Lance Taylor for staying close
> > to the double-precision function:
> > http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_cos.c?rev=1.4;content-type=text%2Fplain
> > 
> > But the degree of the polynomial in the single-precision version
> > is about twice what is necessary, and the coefficients
> > look like coefficients of the Taylor expansion rounded to
> > nearest float, when they should be Chebyshev coefficients.
> 
> Agreed, this is ridiculous.

bruce d evans in freebsd fixed these

http://svnweb.FreeBSD.org/base/stable/9/lib/msun/src/k_cosf.c?view=log
see rev 152951 and 152869

openbsd is a bit behind freebsd libm
improvements in some places
but freebsd misses some long double
functions implemented in openbsd
that's why i recommended looking into
both of them

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-01-23 12:07:15 -0500]:
> On Mon, Jan 23, 2012 at 05:41:52PM +0100, Szabolcs Nagy wrote:
> > i've looked into libm implementations
> > to figure out what's best for musl
> 
> > the extended precision algorithms are reused across
...
> 
> Any ideas how the different ones evolved (separately written or common
> ancestor code, etc.) and where we should look to pull code from?
> 

meanwhile i looked more into libm design issues

here are some questions i come up with:
for background issues see
http://nsz.repo.hu/libm

Code organization:
(ldX is X bit long double)

Do we want ld128?
Should we try to use common code for ld80 and ld128?
How to do ld64: wrap double functions or alias them?
How to tie the ld* code to the arch in the build system?
Make complex optional?
Keep complex in math/ or cmath/?

Workarounds:

Use STDC pragmas (eventhough gcc does not support them)?
Use volatile consistently to avoid evaluation precision and const folding issues?
Use special compiler flags against unwanted optimization (-frounding-math, -ffloat-store)?
Do inline/macro optimization for small functions? (isnan, isinf, signbit, creal, cimag,..)
In complex code prefer creal(), cimag() or a union to (un)pack re,im?

Code cleanups:

Keep diffability with freebsd/openbsd code or reformat the code for clarity?
Keep e_, s_, k_ fdlibm source file name prefixes?
Should 0x1p0 float format be preferred over decimal format?
Should isnan, signbit,.. be preferred over inplace bithacks?
Is unpacking a double into 2 int32_t ok (signed int representation)?
Is unpacking the mantissa of ld80 into an int64_t ok?

Re: libm

[Rich Felker]

On Mon, Feb 27, 2012 at 10:02:53PM +0100, Szabolcs Nagy wrote:
> * Rich Felker <dalias@aerifal.cx> [2012-01-23 12:07:15 -0500]:
> > On Mon, Jan 23, 2012 at 05:41:52PM +0100, Szabolcs Nagy wrote:
> > > i've looked into libm implementations
> > > to figure out what's best for musl
> > 
> > > the extended precision algorithms are reused across
> ....
> > 
> > Any ideas how the different ones evolved (separately written or common
> > ancestor code, etc.) and where we should look to pull code from?
> > 
> 
> meanwhile i looked more into libm design issues
> 
> here are some questions i come up with:
> for background issues see
> http://nsz.repo.hu/libm
> 
> Code organization:
> (ldX is X bit long double)
> 
> Do we want ld128?

Only for systems where long double is 128-bit, and so far, we don't
support any such systems. I'd say it's very low priority right now.

> Should we try to use common code for ld80 and ld128?

This would be nice but I doubt it's feasible, especially without a lot
of divergence from upstream.

> How to do ld64: wrap double functions or alias them?

Wrap. Alias is non-conformant; distinct functions visible to the
application (and not in the reserved namespace) must not have the same
address. I don't mind bending this rule for non-standard functions
that just exist as weak aliases for compatibility with legacy systems,
since conformant applications cannot use them anyway, but for standard
functions this is really a conformance issue. Hopefully the wrapers
will compile to a single tail-call jump instruction anyway.

> How to tie the ld* code to the arch in the build system?

Just put all the files in place, and put #if LDBL_MANT_DIG==... etc.
in the files to make them compile to empty .o files on platforms where
they're not needed.

> Make complex optional?

At this point I don't see a good reason to make anything optional. For
static linking it won't get pulled in unnecessarily anyway, and for
dynamic linking, the whole libc.so is still very small. If size
becomes an issue in the future, I'd like a more generic way of
choosing what symbols to link into the .so at link time rather than
adding lots of #ifdef to the source.

> Keep complex in math/ or cmath/?

I'm fairly indifferent to the issue, especially since with the above
comments in mind it shouldn't affect the build system. If math/ is too
cluttered already, cmath/ may be nice.

> Workarounds:
> 
> Use STDC pragmas (eventhough gcc does not support them)?

There may be another -f option we need to make gcc extra-attentive to
standards for floating point...

> Use volatile consistently to avoid evaluation precision and const
> folding issues?

This might be a nice extra precaution, and might be best if we can't
find any other consistently-working fix.

> Use special compiler flags against unwanted optimization
> (-frounding-math, -ffloat-store)?

I think these are what we want. My understanding is that they behave
similarly to the STDC pragmas. If we add both the STDC pragmas (which
gcc will ignore) and -frounding-math and -ffloat-store, then either
GCC or a C99 conformant compiler should generate correct code.

> Do inline/macro optimization for small functions? (isnan, isinf,
> signbit, creal, cimag,..)
> In complex code prefer creal(), cimag() or a union to (un)pack re,im?

My preference would be to use the standard function-like-macro names
(isnan, creal, etc.) internally. They should expand to something just
as efficient as the union unpacking stuff used internally, but
actually be portable and clear to somebody reading the code.

Of course if it's a lot of gratuitous changes from upstream I'd be
willing to leave it alone, as long as the upstream code is not
invoking UB.

> Code cleanups:
> 
> Keep diffability with freebsd/openbsd code or reformat the code for
> clarity?

> Keep e_, s_, k_ fdlibm source file name prefixes?

I'm fairly indifferent on the matter.

> Should 0x1p0 float format be preferred over decimal format?

I would prefer hex floats so the exact value is clear and so there's
no danger of toolchain bugs causing the constants to be miscompiled.

> Should isnan, signbit,.. be preferred over inplace bithacks?

Yes, but see above..

> Is unpacking a double into 2 int32_t ok (signed int representation)?

Personally I find it really ugly, especially since it introduces an
endian dependency. And there may be UB issues with signed overflow
unless the code is very careful. uint64_t is the ideal unpacking
format for double.

> Is unpacking the mantissa of ld80 into an int64_t ok?

Yes, and this has an unavoidable endian issue... Actually uint64_t
would be better. A valid, finite, normal ld80 will always be negative
as a signed integer, by the way (because the leading 1 is mandatory
but not implicit).

Rich

Re: libm

[Szabolcs Nagy]

see
http://nsz.repo.hu/libm

i put together a libm using the freebsd libm
and some code from openbsd

i did various useless code formatting cleanups
as the code was in bad shape (we will have to
rewrite some parts anyway, so it is not that
important to keep fdlibm diffability, instead
i recorded the date of the freebsd and openbsd
source tree checkout so we can see if they fixed
something since then)

the // FIXME and TODO comments are from me
(i did not mark all the issues, i got tired after a while)


notes:

the source contains parts of musl to make testing/building
easy

i modified math.h and added complex.h and tgmath.h
(although i'm not sure about the later, it uses c11)

i already rewrote some simple functions when they were
incorrect (ilogbl) or overcomplicated (fmax, logb)

internal stuff is in libm.h which includes math.h
and float.h so every .c just includes libm.h

long double is handled in an ugly way now:
ldhack.h contains all the ld macros and unions
that is used in the *l.c functions
for my convenience i defined LD64, LD80 and LD128
there and just used these with #if directive

complex is missing

tgammaf, tgamma is missing
(the only 3clause bsd licensed code)
(cephes lib has another implementation)

i added frexp from musl stdlib for completeness

i did some testing of the double functions
using crlibm and ucbtest testvectors
and so far things seem fine
(pow has some issues in certain corner cases
and it seems x87 hardware sqrt is not correctly
rounded eventhough ieee requires it to be,
or maybe it's another double rounding issue)


questions:

keep frexp in stdlib?

how to do the long double ifdefs?

check x87 fpu precision setting from ld80 code?
(or we assume it's extended precision)

preferred way of handling consts?
(freebsd code uses 1.0f, (float)1, 0x1p0, static const float one=1;..)
(i'm not sure why 'one' or 'zero' is used instead of 1 or 0
maybe it's another compiler float-store bug workaround)
i assume (float)1.0 is the same as 1f

i wonder if it's ok to use a polyeval(coeffs, x) function instead
of writing out the polynomial (some openbsd long double code uses
such thing but not everywhere)


todo:
- move the code to musl tree so the changes are recorded there
- fix int32_t issues (there are many)
- define isnan, signbit, etc macros to be efficient and
change bithacks into isnan etc in the code when appropriate
- fix extended precision issues (probably only fma)

there are lower priority tasks:
use scalbn instead of ldexp in internal code (ldexp is a wrapper)
clean up hard to follow control flows
fix invalid comments
clean up long double code
use consts consistently
make proper libm test suit

Re: libm

[Rich Felker]

On Sat, Mar 03, 2012 at 11:57:58PM +0100, Szabolcs Nagy wrote:
> see
> http://nsz.repo.hu/libm

Some initial questions...

- Why are there two versions (one "slow") of rem_pio2?

- Is __invtrigl.c used? It clashes with the app-reserved namespace.

> i did various useless code formatting cleanups
> as the code was in bad shape (we will have to
> rewrite some parts anyway, so it is not that
> important to keep fdlibm diffability, instead
> i recorded the date of the freebsd and openbsd
> source tree checkout so we can see if they fixed
> something since then)

Sounds fine to me.

> i modified math.h and added complex.h and tgmath.h
> (although i'm not sure about the later, it uses c11)

There's an old gcc way to do tgmath.h without c11, but it's full of
gcc-specific hacks involving __typeof__. If gcc4 allows _Generic even
with -std=c99 etc. then I think _Generic is the best.

> questions:
> 
> keep frexp in stdlib?

It's ok to move it. I had it there from back when there was a separate
libm because printf depends on it.

> how to do the long double ifdefs?

#if LDBL_MANT_DIG == ...

> check x87 fpu precision setting from ld80 code?
> (or we assume it's extended precision)

All existing code (well, mainly macro definitions) in musl assumes its
set to extended, and we provide no code to change it.

> preferred way of handling consts?
> (freebsd code uses 1.0f, (float)1, 0x1p0, static const float one=1;..)
> (i'm not sure why 'one' or 'zero' is used instead of 1 or 0
> maybe it's another compiler float-store bug workaround)
> i assume (float)1.0 is the same as 1f

I like 1.0f, etc. BTW are you sure 1f is valid C? Certainly 1L is not
a valid alternative for 1.0L when you want long double..

> i wonder if it's ok to use a polyeval(coeffs, x) function instead
> of writing out the polynomial (some openbsd long double code uses
> such thing but not everywhere)

I suspect the functions that aren't using it right now depend on a
particular evaluation order to avoid loss of precision. I could be
wrong though; that's just my best guess.

> todo:
> - move the code to musl tree so the changes are recorded there

Very interested in doing this sooner rather than later.

> - fix int32_t issues (there are many)
> - define isnan, signbit, etc macros to be efficient and
> change bithacks into isnan etc in the code when appropriate

Try something like:

union __float_rep { float __x; uint32_t __r; };
#define __FLOAT_REP(x) ((union __float_rep){x}.__r)

#define isnan(x) \
	sizeof(x)==sizeof(float) ? (__FLOAT_REP(x)&0x7fffffff)>0x7f800000 : \
	sizeof(x)==sizeof(double) ? (__DOUBLE_REP(x)&....

Does that work? Of course it would work to just let it get converted
up implicitly to long double and just have one case, but that would be
more expensive I think.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-04 01:53:40 -0500]:
> On Sat, Mar 03, 2012 at 11:57:58PM +0100, Szabolcs Nagy wrote:
> > see
> > http://nsz.repo.hu/libm
> 
> Some initial questions...
> 
> - Why are there two versions (one "slow") of rem_pio2?
> 
originally code they are __ieee754_rem_pio2 and __kernel_rem_pio2
maybe slow is not the best name

the slow version is only called (by the non-slow one)
for large numbers (>2^20) when multi precision arithmetic
is needed using many bits of 2/pi
(there is a large ipio2 array in that file)

maybe it should be called __rem_pio2_large ?

> - Is __invtrigl.c used? It clashes with the app-reserved namespace.
> 
fixed
it is used by long double inv trig functions

> gcc-specific hacks involving __typeof__. If gcc4 allows _Generic even
> with -std=c99 etc. then I think _Generic is the best.
> 
no, gcc does not support it
(i have a gcc 4.7 built from svn in october and it has -std=c1x,
but _Generic is not supported)

i wrote it using generic since i couldnt come up with
any way to solve it in c99

> > how to do the long double ifdefs?
> 
> #if LDBL_MANT_DIG == ...
> 
ok will do that

> > check x87 fpu precision setting from ld80 code?
> > (or we assume it's extended precision)
> 
> All existing code (well, mainly macro definitions) in musl assumes its
> set to extended, and we provide no code to change it.
> 
ok

> > preferred way of handling consts?
> > (freebsd code uses 1.0f, (float)1, 0x1p0, static const float one=1;..)
> > (i'm not sure why 'one' or 'zero' is used instead of 1 or 0
> > maybe it's another compiler float-store bug workaround)
> > i assume (float)1.0 is the same as 1f
> 
> I like 1.0f, etc. BTW are you sure 1f is valid C? Certainly 1L is not
> a valid alternative for 1.0L when you want long double..
> 
ok
(1f might be wrong, i just wrote that in the mail)

> > i wonder if it's ok to use a polyeval(coeffs, x) function instead
> > of writing out the polynomial (some openbsd long double code uses
> > such thing but not everywhere)
> 
> I suspect the functions that aren't using it right now depend on a
> particular evaluation order to avoid loss of precision. I could be
> wrong though; that's just my best guess.
> 
it's just a minor detail, but this was one of the parts where
things could be a bit nicer and more consistent (probably)

(coeffs in array, and using some static inline eval etc,
it might not be a good idea for some technical reason,
but i dont yet see why)

> > todo:
> > - move the code to musl tree so the changes are recorded there
> 
> Very interested in doing this sooner rather than later.
> 
ok
complex is missing but it is not hard to add
(if you want to wait for it)

the parts i wrote may need some review or testing
(it can be done after merging of course)

(it would be nice if musl had a dedicated testing
package where one could add test functions
and easy to run)

> > - fix int32_t issues (there are many)
> > - define isnan, signbit, etc macros to be efficient and
> > change bithacks into isnan etc in the code when appropriate
> 
> Try something like:
> 
> union __float_rep { float __x; uint32_t __r; };
> #define __FLOAT_REP(x) ((union __float_rep){x}.__r)
> 
> #define isnan(x) \
> 	sizeof(x)==sizeof(float) ? (__FLOAT_REP(x)&0x7fffffff)>0x7f800000 : \
> 	sizeof(x)==sizeof(double) ? (__DOUBLE_REP(x)&....
> 
> Does that work? Of course it would work to just let it get converted
> up implicitly to long double and just have one case, but that would be
> more expensive I think.
>
ok so using compound literal
i was thinking about this
something like this will work

the difficult part is figuring out if it can be used internally

(many math functions start by turning the float/double
into unsigned int, then doing various checks to handle
special cases like nan, inf, signbit, etc using bithacks

i thought most of them could be isnan() etc
eg compare the current fmax to the freebsd one

i'm not sure if this can be done in general
because some bit hacks will be needed anyway
and it might be better doing the conversion
once

many functions do checks like
  if(x < 12.34) ... else if(x < 78.9) ...
but using unsigned int arithmetics like
  if ((x&bitmask) < 0xabcdef) ...
which is not clear, but i guess this
cannot be changed because float compare
does a different thing
but it would be nice to reduce the bithacks
without changing the semantics or performance)

Re: libm

[Rich Felker]

On Sun, Mar 04, 2012 at 03:50:41PM +0100, Szabolcs Nagy wrote:
> > - Why are there two versions (one "slow") of rem_pio2?
> > 
> originally code they are __ieee754_rem_pio2 and __kernel_rem_pio2
> maybe slow is not the best name
> 
> the slow version is only called (by the non-slow one)
> for large numbers (>2^20) when multi precision arithmetic
> is needed using many bits of 2/pi
> (there is a large ipio2 array in that file)

Yes, that left me rather confused. The original names are confusing
too. In both cases it sounds like you have 2 versions of the same code
(fast/slow or different conformance profiles) rather than 2 different
functions. 

> maybe it should be called __rem_pio2_large ?

That's better. Not sure what the best name would be..

> > gcc-specific hacks involving __typeof__. If gcc4 allows _Generic even
> > with -std=c99 etc. then I think _Generic is the best.
> > 
> no, gcc does not support it
> (i have a gcc 4.7 built from svn in october and it has -std=c1x,
> but _Generic is not supported)
> 
> i wrote it using generic since i couldnt come up with
> any way to solve it in c99

There's no way to solve it in C99 matching the return types, but if
you don't care about return types, it's pretty easy (albeit ugly and
verbose) with ?: conditionals based on sizeof(x).

Since gcc still doesn't support _Generic, I think it might be
preferable to use the old approach (sizeof tests) and add (with a
conditionally defined macro) a cast to __typeof__((x)+1.0f) (or
__typeof__((x)+(y)+1.0f), etc. as needed) on the result that's only
enabled when __GNUC__ is defined. Then the fallback case would fail to
have the right result types, but at least the right functions would
get called.

> > > i wonder if it's ok to use a polyeval(coeffs, x) function instead
> > > of writing out the polynomial (some openbsd long double code uses
> > > such thing but not everywhere)
> > 
> > I suspect the functions that aren't using it right now depend on a
> > particular evaluation order to avoid loss of precision. I could be
> > wrong though; that's just my best guess.
> > 
> it's just a minor detail, but this was one of the parts where
> things could be a bit nicer and more consistent (probably)
> 
> (coeffs in array, and using some static inline eval etc,
> it might not be a good idea for some technical reason,
> but i dont yet see why)

If it's not going to change the behavior of the code and it's just a
matter of inlining/un-inlining, I'm fairly indifferent to how you want
to do it.

> > > todo:
> > > - move the code to musl tree so the changes are recorded there
> > 
> > Very interested in doing this sooner rather than later.
> > 
> ok
> complex is missing but it is not hard to add
> (if you want to wait for it)

Either way.

> the parts i wrote may need some review or testing
> (it can be done after merging of course)

Yes, I'm okay with that.

> (it would be nice if musl had a dedicated testing
> package where one could add test functions
> and easy to run)

There's the libc-testsuite git repo, and then cluts. Neither is
musl-specific. The former is much broader in scope; the latter is
probably more intensive in the few areas it tests.

> ok so using compound literal
> i was thinking about this
> something like this will work

OK.

> the difficult part is figuring out if it can be used internally
> 
> (many math functions start by turning the float/double
> into unsigned int, then doing various checks to handle
> special cases like nan, inf, signbit, etc using bithacks
> 
> i thought most of them could be isnan() etc
> eg compare the current fmax to the freebsd one
> 
> i'm not sure if this can be done in general
> because some bit hacks will be needed anyway
> and it might be better doing the conversion
> once

Yes, at least it would be a major undertaking to convert all the code
and ensure that no bugs are introduced. I think it's fine leaving it
how it is for now. I still want the fast isnan, etc. for applications
to use, though (and for use by other parts of libc like printf).

> many functions do checks like
>   if(x < 12.34) ... else if(x < 78.9) ...
> but using unsigned int arithmetics like
>   if ((x&bitmask) < 0xabcdef) ...
> which is not clear, but i guess this
> cannot be changed because float compare
> does a different thing
> but it would be nice to reduce the bithacks
> without changing the semantics or performance)

The semantics should only be changed if the goal is to consider
NAN/INF as "large" values too. Of course you also get abs() for free
using the representation and bitmask. Otherwise, for finite
comparisons comparing the representations and comparing the
floating-point values should have the same behavior.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-04 13:43:39 -0500]:
> On Sun, Mar 04, 2012 at 03:50:41PM +0100, Szabolcs Nagy wrote:
> > maybe it should be called __rem_pio2_large ?
> 
> That's better. Not sure what the best name would be..
> 
renamed

> > i wrote it using generic since i couldnt come up with
> > any way to solve it in c99
> 
> There's no way to solve it in C99 matching the return types, but if
> you don't care about return types, it's pretty easy (albeit ugly and
> verbose) with ?: conditionals based on sizeof(x).
> 
hm i haven't even thought about return types

my problem was calling the double version of a function when the
arg is integral and the float version when the arg is float.
(it can be done with sizeof if we can assume there is a larger
integer type than float, but i thought that was ugly)

even if we don't care about return types and int arguments
it seems tgmath is very ugly

> > complex is missing but it is not hard to add
> > (if you want to wait for it)
> 
> Either way.
> 
> > the parts i wrote may need some review or testing
> > (it can be done after merging of course)
> 
> Yes, I'm okay with that.
> 

i can do the merge and then you git pull
or you can do it from libm

the long double ifdefs needs to be fixed
and math.h sorted out

> There's the libc-testsuite git repo, and then cluts. Neither is
> musl-specific. The former is much broader in scope; the latter is
> probably more intensive in the few areas it tests.
> 

i know but i'm not content with either of those
imho there should be one libctest repo which is a bit more
organized than libc-testsuit and contains the cluts tests
as well

> > ok so using compound literal
> > i was thinking about this
> > something like this will work
> 
> OK.
> 

i saw that you removed the compound literal definition of
NAN, INFINITY etc from math.h

if we want to keep the headers c89 friendly then
isnan, isinf,.. cannot be defined without function call
maybe use static __isnanf,.. functions?

(i committed a compound literal based implementatoin for now)

> Yes, at least it would be a major undertaking to convert all the code
> and ensure that no bugs are introduced. I think it's fine leaving it
> how it is for now. I still want the fast isnan, etc. for applications
> to use, though (and for use by other parts of libc like printf).
> 

i agree

> The semantics should only be changed if the goal is to consider
> NAN/INF as "large" values too. Of course you also get abs() for free
> using the representation and bitmask. Otherwise, for finite
> comparisons comparing the representations and comparing the
> floating-point values should have the same behavior.
> 
hm but float compare will compile to different instruction
so there is at least performance difference

Re: libm

[Szabolcs Nagy]

* Szabolcs Nagy <nsz@port70.net> [2012-03-04 15:50:41 +0100]:
> * Rich Felker <dalias@aerifal.cx> [2012-03-04 01:53:40 -0500]:
> > On Sat, Mar 03, 2012 at 11:57:58PM +0100, Szabolcs Nagy wrote:
> > > how to do the long double ifdefs?
> > 
> > #if LDBL_MANT_DIG == ...
> > 
> ok will do that
> 

i use the following checks in all long double code now
LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384

the inplace union declarations for bit hacks are not ok
anymore because endianness is not checked

the union declarations should go to an arch specific .h
eg. arch/i386/longdouble.h

but before that long double bit hacks may need some
cleanups (freebsd union vs openbsd macros vs ..)

Re: libm

[Rich Felker]

On Mon, Mar 05, 2012 at 09:51:35AM +0100, Szabolcs Nagy wrote:
> > > i wrote it using generic since i couldnt come up with
> > > any way to solve it in c99
> > 
> > There's no way to solve it in C99 matching the return types, but if
> > you don't care about return types, it's pretty easy (albeit ugly and
> > verbose) with ?: conditionals based on sizeof(x).
> > 
> hm i haven't even thought about return types
> 
> my problem was calling the double version of a function when the
> arg is integral and the float version when the arg is float.
> (it can be done with sizeof if we can assume there is a larger
> integer type than float, but i thought that was ugly)

Hm? If you need a macro to test whether an argument is an integer or
floating point, try this:

#define __IS_FP(x) !!((1?1:(x))/2)

No assumptions at all.

> even if we don't care about return types and int arguments
> it seems tgmath is very ugly

Yes, it is ugly. I don't think there's any way around the fact that
tgmath.h is ugly.

> > > complex is missing but it is not hard to add
> > > (if you want to wait for it)
> > 
> > Either way.
> > 
> > > the parts i wrote may need some review or testing
> > > (it can be done after merging of course)
> > 
> > Yes, I'm okay with that.
> > 
> 
> i can do the merge and then you git pull
> or you can do it from libm
> 
> the long double ifdefs needs to be fixed
> and math.h sorted out

I can do it, but might be another day or two. Feel free to keep fixing
stuff up in the mean time if you have time. :)

> > There's the libc-testsuite git repo, and then cluts. Neither is
> > musl-specific. The former is much broader in scope; the latter is
> > probably more intensive in the few areas it tests.
> > 
> 
> i know but i'm not content with either of those
> imho there should be one libctest repo which is a bit more
> organized than libc-testsuit and contains the cluts tests
> as well
> 
> > > ok so using compound literal
> > > i was thinking about this
> > > something like this will work
> > 
> > OK.
> > 
> 
> i saw that you removed the compound literal definition of
> NAN, INFINITY etc from math.h
> 
> if we want to keep the headers c89 friendly then
> isnan, isinf,.. cannot be defined without function call
> maybe use static __isnanf,.. functions?

The removal has nothing to do with c89; it's actually the fact that
the compound literal definition was not a constant expression and the
standard requires these macros to expand to constant expressions.

As far as I know, gcc allows the compound literal even in old-standard
modes, probably because it was a feature of "GNU C" long before C99
adopted it. In any case I see no good way to define them without
compound literals except the function calls..

> > The semantics should only be changed if the goal is to consider
> > NAN/INF as "large" values too. Of course you also get abs() for free
> > using the representation and bitmask. Otherwise, for finite
> > comparisons comparing the representations and comparing the
> > floating-point values should have the same behavior.
> > 
> hm but float compare will compile to different instruction
> so there is at least performance difference

Yes. In principle it would be faster to do everything as floating
point so no store/load delays are needed. In relality I doubt it makes
much of a difference. In almost all of these functions the actual
computation time dominates.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-05 09:04:59 -0500]:
> On Mon, Mar 05, 2012 at 09:51:35AM +0100, Szabolcs Nagy wrote:
> Hm? If you need a macro to test whether an argument is an integer or
> floating point, try this:
> 
> #define __IS_FP(x) !!((1?1:(x))/2)
> 

oh nice
(and subtle)

> I can do it, but might be another day or two. Feel free to keep fixing
> stuff up in the mean time if you have time. :)
> 

ok

> > i saw that you removed the compound literal definition of
> > NAN, INFINITY etc from math.h
> > 
> 
> The removal has nothing to do with c89; it's actually the fact that
> the compound literal definition was not a constant expression and the
> standard requires these macros to expand to constant expressions.
> 

i see

> adopted it. In any case I see no good way to define them without
> compound literals except the function calls..
> 

btw while testing these macros i noticed that when
multiple classification macros are used fpclassify
gets called many times
(with the previous solution)

the extra calls could be optimized by adding
__attribute__((const)) to fpclassify
(resulted less calls, smaller code size etc)

> > hm but float compare will compile to different instruction
> > so there is at least performance difference
> 
> Yes. In principle it would be faster to do everything as floating
> point so no store/load delays are needed. In relality I doubt it makes
> much of a difference. In almost all of these functions the actual
> computation time dominates.
> 

i see

Re: libm

[Rich Felker]

On Mon, Mar 05, 2012 at 04:17:11PM +0100, Szabolcs Nagy wrote:
> * Rich Felker <dalias@aerifal.cx> [2012-03-05 09:04:59 -0500]:
> > On Mon, Mar 05, 2012 at 09:51:35AM +0100, Szabolcs Nagy wrote:
> > Hm? If you need a macro to test whether an argument is an integer or
> > floating point, try this:
> > 
> > #define __IS_FP(x) !!((1?1:(x))/2)
> 
> oh nice
> (and subtle)

Note that it avoids evaluating x too, so it's safe in macros that need
to avoid evaluating the argument more than once.

> > adopted it. In any case I see no good way to define them without
> > compound literals except the function calls..
> > 
> 
> btw while testing these macros i noticed that when
> multiple classification macros are used fpclassify
> gets called many times
> (with the previous solution)
> 
> the extra calls could be optimized by adding
> __attribute__((const)) to fpclassify
> (resulted less calls, smaller code size etc)

Indeed. I was wrongly thinking they'd need __attribute__((__pure__))
rather than const, which is trickier because it's only available in
newish gcc. But const has been around since before 2.95 even, so just
#ifdef __GNUC__ is a sufficient check for using it.

Alternatively I might just move the definitions of fpclassify (at
least for non-LD types) back into math.h as static inline. This would
let gcc optimize it down a lot more.

Rich

Re: libm

[Rich Felker]

On Mon, Mar 05, 2012 at 04:17:11PM +0100, Szabolcs Nagy wrote:
> * Rich Felker <dalias@aerifal.cx> [2012-03-05 09:04:59 -0500]:
> > On Mon, Mar 05, 2012 at 09:51:35AM +0100, Szabolcs Nagy wrote:
> > Hm? If you need a macro to test whether an argument is an integer or
> > floating point, try this:
> > 
> > #define __IS_FP(x) !!((1?1:(x))/2)
> > 
> 
> oh nice
> (and subtle)

I checked out your latest tgmath.h; it's surprisingly clean and
readable! A couple ideas..

1. __fun could be just fun (macro arguments are free of namespace
issues)

2. You could add something like:
#ifdef __GNUC__
#define __RETCAST(x) (__typeof__((x))
#else
#define __RETCAST(x)
#endif
Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
needed to make integer types result in a cast to double, though.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-09 05:22:39 -0500]:

> On Mon, Mar 05, 2012 at 04:17:11PM +0100, Szabolcs Nagy wrote:
> > * Rich Felker <dalias@aerifal.cx> [2012-03-05 09:04:59 -0500]:
> > > On Mon, Mar 05, 2012 at 09:51:35AM +0100, Szabolcs Nagy wrote:
> > > Hm? If you need a macro to test whether an argument is an integer or
> > > floating point, try this:
> > > 
> > > #define __IS_FP(x) !!((1?1:(x))/2)
> > > 
> > 
> > oh nice
> > (and subtle)
> 
> I checked out your latest tgmath.h; it's surprisingly clean and
> readable! A couple ideas..
> 
> 1. __fun could be just fun (macro arguments are free of namespace
> issues)
> 

hm true

> 2. You could add something like:
> #ifdef __GNUC__
> #define __RETCAST(x) (__typeof__((x))
> #else
> #define __RETCAST(x)
> #endif
> Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> needed to make integer types result in a cast to double, though.
> 

ok, that makes sense

actually i'm not sure when the extra () protection
is needed like __typeof__((x))
(the reason for it in (x)+(y) is clear)

is it because single argument macros might get
called with a comma expression?

so
#define A(x) B((x))
is ok

but
#define A(x,y) B((x),(y))
is redundant

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-09 05:22:39 -0500]:
> 2. You could add something like:
> #ifdef __GNUC__
> #define __RETCAST(x) (__typeof__((x))
> #else
> #define __RETCAST(x)
> #endif
> Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> needed to make integer types result in a cast to double, though.

hm the int->double and complex/real cases are tricky

i thought +1.0 or +I would solve these, but that's wrong

Re: libm

[Rich Felker]

On Fri, Mar 09, 2012 at 12:09:46PM +0100, Szabolcs Nagy wrote:
> * Rich Felker <dalias@aerifal.cx> [2012-03-09 05:22:39 -0500]:
> > 2. You could add something like:
> > #ifdef __GNUC__
> > #define __RETCAST(x) (__typeof__((x))
> > #else
> > #define __RETCAST(x)
> > #endif
> > Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> > needed to make integer types result in a cast to double, though.
> 
> hm the int->double and complex/real cases are tricky
> 
> i thought +1.0 or +I would solve these, but that's wrong

I think +0.0f might solve it. Isn't the promoted type of any integer
type with float double?

Rich

Re: libm

[Rich Felker]

On Fri, Mar 09, 2012 at 11:57:03AM +0100, Szabolcs Nagy wrote:
> > 2. You could add something like:
> > #ifdef __GNUC__
> > #define __RETCAST(x) (__typeof__((x))
> > #else
> > #define __RETCAST(x)
> > #endif
> > Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> > needed to make integer types result in a cast to double, though.
> > 
> 
> ok, that makes sense
> 
> actually i'm not sure when the extra () protection
> is needed like __typeof__((x))
> (the reason for it in (x)+(y) is clear)
> 
> is it because single argument macros might get
> called with a comma expression?
> 
> so
> #define A(x) B((x))
> is ok
> 
> but
> #define A(x,y) B((x),(y))
> is redundant

It's probably redundant. Only case I can think of where it might help
is improving(?) the error reporting when some invalid-as-macro-arg
expressions that use ?: as a grouping for a comma operator in the
second operand, or perhaps some fishy stuff with compound literals.

Rich

Re: libm

[Rich Felker]

On Fri, Mar 09, 2012 at 10:56:55AM -0500, Rich Felker wrote:
> > > Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> > > needed to make integer types result in a cast to double, though.
> > 
> > hm the int->double and complex/real cases are tricky
> > 
> > i thought +1.0 or +I would solve these, but that's wrong
> 
> I think +0.0f might solve it. Isn't the promoted type of any integer
> type with float double?

Nope, it results in float. glibc has a trick involving obscure
interactions of null pointer constants and the ?: operator that
generates the right type using __typeof__, which is no big deal
because this code is already only used when __typeof__ is available,
anyway.

There's a good blog article on it somewhere which I can't seem to find
at the moment...

Rich

Re: libm

[Rich Felker]

On Fri, Mar 09, 2012 at 12:02:54PM -0500, Rich Felker wrote:
> On Fri, Mar 09, 2012 at 10:56:55AM -0500, Rich Felker wrote:
> > > > Then add __RETCAST((x)), __RETCAST((x)+(y)), etc. Some trick will be
> > > > needed to make integer types result in a cast to double, though.
> > > 
> > > hm the int->double and complex/real cases are tricky
> > > 
> > > i thought +1.0 or +I would solve these, but that's wrong
> > 
> > I think +0.0f might solve it. Isn't the promoted type of any integer
> > type with float double?
> 
> Nope, it results in float. glibc has a trick involving obscure
> interactions of null pointer constants and the ?: operator that
> generates the right type using __typeof__, which is no big deal
> because this code is already only used when __typeof__ is available,
> anyway.
> 
> There's a good blog article on it somewhere which I can't seem to find
> at the moment...

Try this:

#define __RETCAST(x) (__typeof__(*(0 \
? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 \
: (__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))

In the outer conditional operator, the second operand is a pointer to
double if (void *)__IS_FP(x) is a null pointer constant (which is true
iff __IS_FP(x)==0) and otherwise is a pointer to void.

Conversely, the third operand of the outer conditional is a pointer to
__typeof__(x) iff __IS_FP(x)!=0, and otherwise a pointer to void.

Thus the outer conditional sees either "double *" and "void *", or
else "void *" and "__typeof__(x) *", resulting it in having type
"double *" or "__typeof__(x)", reflecting whether x was an integer or
floating point type.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-09 22:28:44 -0500]:
> Try this:
> 
> #define __RETCAST(x) (__typeof__(*(0 \
> ? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 \
> : (__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))
> 
> In the outer conditional operator, the second operand is a pointer to
> double if (void *)__IS_FP(x) is a null pointer constant (which is true
> iff __IS_FP(x)==0) and otherwise is a pointer to void.
> 
> Conversely, the third operand of the outer conditional is a pointer to
> __typeof__(x) iff __IS_FP(x)!=0, and otherwise a pointer to void.
> 
> Thus the outer conditional sees either "double *" and "void *", or
> else "void *" and "__typeof__(x) *", resulting it in having type
> "double *" or "__typeof__(x)", reflecting whether x was an integer or
> floating point type.
> 

of course there was another twist:
some complex functions return real value :)

i solved that with a __TO_REAL(x) macro using similar tricks
so __RETCAST(__TO_REAL(x)) gives correct return type


there are two remaining problems:

for two argument functions __RETCAST((x)+(y)) is not ok
if x is float, y is int then the return type should be double, not float
(this can be solved by a __RETCAST2(x,y) i guess)

if sizeof(long double) == sizeof(double) then the selected
function will be wrong
(it should not matter much, we can leave it as is or enforce double)

Re: libm

[Rich Felker]

On Sat, Mar 10, 2012 at 01:45:53PM +0100, Szabolcs Nagy wrote:
> there are two remaining problems:
> 
> for two argument functions __RETCAST((x)+(y)) is not ok
> if x is float, y is int then the return type should be double, not float
> (this can be solved by a __RETCAST2(x,y) i guess)

I don't see any better way right off...
Actually it's really ugly that C's promotions don't respect the size
of the integer type when combining with floating point types. x+0.0f
where x is an integer type is almost sure to destroy the value of x
except for very small values. C really should have made the promotion
be to the smallest floating point type that has at least as many bits
in its mantissa as value bits in the integer type, or long double if
no such type exists.

Or from the other side (with C as it is), programmers should be really
careful about implicit conversions to floating point types from
integer types.

> if sizeof(long double) == sizeof(double) then the selected
> function will be wrong
> (it should not matter much, we can leave it as is or enforce double)

It would probably be nice to put the case for double before long
double to avoid going through an extra wrapper function on LD64
systems, but it's not a big deal.

Rich

Re: libm

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2012-03-10 08:12:45 -0500]:

> On Sat, Mar 10, 2012 at 01:45:53PM +0100, Szabolcs Nagy wrote:
> > for two argument functions __RETCAST((x)+(y)) is not ok
> > if x is float, y is int then the return type should be double, not float
> > (this can be solved by a __RETCAST2(x,y) i guess)
> 
> I don't see any better way right off...

i committed a solution
handling all the combinations turned out to be quite tricky


> Actually it's really ugly that C's promotions don't respect the size
> of the integer type when combining with floating point types. x+0.0f

i'd say that the ugly thing is that they added tgmath.h
eventhough there is no good way to do type generics in c

the implicit type conversions and promotions are already
hard to tackle in c so one should be careful anyway
(and use explicit casts when not sure)

> It would probably be nice to put the case for double before long
> double to avoid going through an extra wrapper function on LD64
> systems, but it's not a big deal.

i made the long double conditions false whenever
sizeof(double) == sizeof(long double)
so on LD64 only double functions will be used