Use of C99 extra long double math functions after r236148

[Rainer Hurling]

On 11.07.2012 00:58 (UTC+2), David Schultz wrote:
> On Tue, Jul 10, 2012, Rainer Hurling wrote:
>> On 10.07.2012 17:11 (UTC+2), David Schultz wrote:
>>> On Tue, Jul 10, 2012, Rainer Hurling wrote:
>>>> On 10.07.2012 16:02 (UTC+2), Warner Losh wrote:
>>>>>
>>>>> On Jul 10, 2012, at 3:10 AM, Rainer Hurling wrote:
>>>>>> As far as I understand from discussions on R mailing list
>>>>>> (r-devel at r-project.org), they plan to reduce the emulation and/or
>>>>>> workaround of long and complex math functions for FreeBSD and other
>>>>>> systems with their next releases of R devel. So we could really need
>>>>>> some
>>>>>> progress with our C99 conform math functions ;-)
>>>>>
>>>>> Not having R would be a bit pain in my backside.  That's one of the
>>>>> practical considerations that I was talking about.  It is very real, and
>>>>> if I have to, I'll commit the #define junk I railed against to get it
>>>>> back.  Please, let's get some progress.  I have some time to help.
>>>>
>>>> Yes, thank you Warner, that is also my problem. As I wrote some weeks
>>>> ago (05/28/2012) when starting this thread, I am using FreeBSD as a
>>>> scientific desktop because of its good scaling properties. For some
>>>> years now, FreeBSD fits all our needs with R, SAGA GIS, PostgreSQL and
>>>> some more.
>>>>
>>>> If I would not be able to run upcoming versions of R on FreeBSD any
>>>> more, that would be really, really hard :-(
>>>
>>> Do you have a list of the essential functions here?  There are 17 long
>>> double functions and some complex functions missing, but only a
>>> handful of those are of general interest.  The reason I ask is that if
>>> R is just looking for a few missing functions that are already mostly
>>> implemented, then the best solution is probably to finish that work.
>>> But if it's expecting us to have something arcane like long double
>>> Bessel functions of the first kind, then we need to pursue a workaround
>>> in the short term.
>>>
>>
>> That is, what I found by grepping the sources of a recent R development
>> version:
>>
>> expl:   src/nmath/pnchisq.c

Many thanks to you three for implementing expl() with r238722 and r238724.

I am not a C programmer, but would like to ask if the following example 
is correct and suituable as a minimalistic test of this new C99 function?


//-----------------------------------
#include <stdio.h>
#include <math.h>

int main(void)
{
   double c = 2.0;
   long double d = 2.0;

   double e = exp(c);
   long double f = expl(d);

   printf("exp(%f)  is %.*f\n",  c, 90, e);
   printf("expl(%Lf) is %.*Lf\n", d, 90, f);

   return 0;
}
//-----------------------------------


Compiled with 'c99 -o math_expl math_expl.c -lm' and running afterwards 
it gives me:

exp(2.000000)  is 
7.389056098930650406941822438966482877731323242187500000000000000000000000000000000000000000
expl(2.000000) is 
7.389056098930650227397942675366948606097139418125152587890625000000000000000000000000000000


Already the sixteenth position after decimal point differs.

Please forgive me, if this is a really stupid approach for producing 
some expl() output.

uname -a: FreeBSD xxx.xxx.xxx 10.0-CURRENT FreeBSD 10.0-CURRENT #0 
r238740: Tue Jul 24 18:08:13 CEST 2012 
xxx at xxx.xxx.xxx:/usr/obj/usr/src/sys/XXX  amd64

Rainer


>> logl:   src/nmath/dnbeta.c
>>          src/nmath/pnbeta.c
>
> Bruce has versions of these that could be committed with some cleanup.
> It's a matter of sorting through about 1200 emails from him and 3
> source trees to find the most up-to-date patches, then cleaning them
> up and testing and committing them.  I have no time right now, but I
> will do at least the first step as soon as I can, and try to get the
> patches to someone willing to do the final few steps.
>
>> log10l: src/extra/trio/trio.c
>>
>> log1pl: src/nmath/pnbeta.c
>
> If Bruce doesn't already have implementations of these, they are easy
> wrappers around logl() or some internal k_logl in Bruce's implementation.
>
>> powl:   src/extra/trio/triostr.c
>>          src/extra/trio/trio.c
>>          src/main/format.c
>
> It's hard to do a good job on powl(), but the simple approach
> (exp(log(x)*y)) plus a few special cases may suffice for many uses.
>
>> NEWS:l2044
>> The C99 functions acosh, asinh, atanh, snprintf and vsnprintf are
>> now required.
>
> We have had them forever.
>
>> NEWS:l3032
>> The C99 double complex type is now required.
>> The C99 complex trigonometric functions (such as csin) are not
>> currently required (FreeBSD lacks most of them): substitutes are
>> used if they are missing.
>
> We have these (but not the inverse trig functions).
>
>> NEWS:l3277
>> Complex arithmetic (notably z^n for complex z and integer n) gave
>> incorrect results since R 2.10.0 on platforms without C99 complex
>> support.  This and some lesser issues in trigonometric functions
>> have been corrected.
>> Such platforms were rare (we know of Cygwin and FreeBSD).
>> However, because of new compiler optimizations in the way complex
>> arguments are handled, the same code was selected on x86_64 Linux
>> with gcc 4.5.x at the default -O2 optimization (but not at -O).
>
> Not sure if this is relevant.
>
>> BTW: There seems to be a discrepancy about missing functions listed in
>> http://wiki.freebsd.org/MissingMathStuff and in
>> http://svnweb.freebsd.org/base/head/lib/msun/src/math.h?r1=227472&r2=236148&pathrev=236148.
>> So the wiki is a bit outdated now?
>
> My list:
>
> REAL FUNCTIONS (17):
>
> long double     log2l(long double);
> long double     logl(long double);
> long double     log1pl(long double);
>    long double     acoshl(long double);
>    long double     asinhl(long double);
>    long double     atanhl(long double);
>    long double     log10l(long double);
>
> long double     expl(long double);
> long double     expm1l(long double);
>    long double     coshl(long double);
>    long double     sinhl(long double);
>    long double     tanhl(long double);
>    long double     erfcl(long double);
>    long double     erfl(long double);
>
> long double     powl(long double, long double);
>
> long double     lgammal(long double);
> long double     tgammal(long double);
>
>
> COMPLEX FUNCTIONS (37):
>
> long double complex cexpl(long double complex);
> long double complex ccosl(long double complex);
>    long double complex ccoshl(long double complex);
> long double complex csinl(long double complex);
>    long double complex csinhl(long double complex);
> long double complex ctanl(long double complex);
>    long double complex ctanhl(long double complex);
>
>
> float complex cacosf(float complex);
>    float complex cacoshf(float complex);
> double complex cacos(double complex);
>    double complex cacosh(double complex);
> long double complex cacosl(long double complex);
>    long double complex cacoshl(long double complex);
>
> float complex casinf(float complex);
>    float complex casinhf(float complex);
> double complex casin(double complex);
>    double complex casinh(double complex);
> long double complex casinl(long double complex);
>    long double complex casinhl(long double complex);
>
> float complex catanf(float complex);
>    float complex catanhf(float complex);
> double complex catan(double complex);
>    double complex catanh(double complex);
> long double complex catanl(long double complex);
>    long double complex catanhl(long double complex);
>
> float complex clogf(float complex);
> double complex clog(double complex);
> long double complex clogl(long double complex);
>
> float complex cpowf(float complex, float complex);
> double complex cpow(double complex, double complex);
> long double complex cpowl(long double complex, long double complex);
>