long double versions of functions in gcc under Cygwin

[lcw1964]

Greetings, all,

I am trying to port a little bit of math code to gcc, that in the
original version used the long double version of several functions (in
particular, atanl, fabsl, and expl).

I get a complie-time "unidentifi ed reference" error to the expl()
calls, but gcc seems to digest atanl and fabsl just fine. Changing expl
to exp cures the compile time problem, but I get at best double
precision in the final results. I am assuming that the use of exp() vs.
expl() is the weak link.

The GCC documentation seems to imply that expl() is supported, but I
have no idea where to find it or how to link it in properly. For that
matter, I can't seem to find prototypes in math.h for fabsl or atanl,
and they don't make gcc cough at all.

I hope this tenderfoot can find some direction, or I may resort to
singing the praises of the egregiously un-portable lcc-win32 with its
impressive 100+ digit precision qfloat library ;)

cheers,

Les

p.s. I am trying to keep this simple, so if there is a solution within
the main gcc offerings without me having to turn to the GSL, I would
like to try that first.

[Keith Thompson]

"lcw1964" <le***********@ alumni.uwo.cawr ites:

jacob navia wrote:

[...]

>Obviously this is a bug. You could report it to them, maybe they
are interested in knowing about it. There must be some mailing
list in the cygwin docs. I was subscribed ages ago.

Under linux:
[root@gateway tmp]# cat texpl.c
#include <stdio.h>
#include <math.h>
int main(void)
{
long double n = expl(1.0L);
printf("%Lg\n", n);
return 0;
}
[root@gateway tmp]# gcc texpl.c -lm
[root@gateway tmp]# ./a.out
2.71828
[root@gateway tmp]#

this works, so it must be a bug in the cygwin environment/library.

Messr. Navia, that worked for me.

The key is the addition of the -lm parameter to the command line, which
I did not have before.

I suspect that I have made an embarassing beginner's error and have
stirred up a lot of hubbub unnecessarily!

I have no idea what those three little characters mean (-lm), but
something tells me that before I post my next newbie question I do a
little more homework first.

That's a reasonable conclusion, but I'm afraid it happens to be wrong.

The "-lm" option, in some (mostly Unix-like) implementations , tells
the linker to link in the math library. This is question 14.3 in the
comp.lang.c FAQ, <http://www.c-faq.com/>.

But that's not what's going on here. The following is mostly specific
to Cygwin, and therefore only marginally topical.

Cygwin's implementation happens to be smart enough that it doesn't
need the "-lm" option; it links the math library without being asked
if it needs it.

jacob's program works for me under Cygwin, but only because the
compiler itself is smart enough to replace the expression expl(1.0L)
with its value during compilation. If you change the line
long double n = expl(1.0L);
to
long double one = 1.0L;
long double n = expl(one)
then it fails. If you declare one as "const" and compile with "-O1"
or higher, then it works again.

The compiler is capable of evaluating expl() in some very limited
circumstances. If that fails, it generates a call to expl(), which
doesn't exist in the runtime library (regardless of whether you use
"-lm").

The point is that the compiler and the runtime libraries are two
different part of the implementation. If they're provided separately,
you can see odd behavior if one of them supports a given feature and
the other doesn't.

If "double" precision is good enough, use exp(). If not you'll need
to find another solution. (You might be able to use exp() to create a
close approximation to the correct long double result, but I don't
know how to refine that to the required precision.)

There are open-source C libraries including glibc. You might be able
to extract an expl() implementation from one of them. Google is your
friend.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[Keith Thompson]

"Dann Corbit" <dc*****@connx. comwrites:

"jaysome" <ja*****@spamco p.netwrote in message
news:58******** *************** *********@4ax.c om...

[...]

>There are 477 usages of "goto" in this source. That gives me a queasy
feeling.

I'm in the Knuth camp:
http://portal.acm.org/citation.cfm?i...FTOKEN=6184618

[...]

The link is to Knuth's article "Structured Programming with go to
Statements", from the December 1974 issue of ACM Computing Surveys.

In my humble opinion, most instances of goto statements point to a
missing feature in the language. In C, the major missing features are
named break (which could be used to break out of a specified loop
rather than the nearest enclosing one) and a decent exception
mechanism.

All code can be written with conditional and unconditional gotos; all
the structured control constructs are, in a sense, syntactic sugar.
In my opinion, C would benefit from just a couple more such constructs
that could nearly eliminate the need for gotos. (Of course it would
be a very long time before any such change could become sufficiently
widely supported that we could actually use it.)

In the meantime, here in the real world, a goto is sometimes the
cleanest available solution to a problem. (I haven't studied the code
in question, so I can't comment on whether it turned out to be the
cleanest solution 477 times.)

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[lcw1964]

pete wrote:

lcw1964 wrote:

expl()

#include <float.h>

long double fs_expl(long double x);
long double fs_logl(long double x);
long double fs_sqrtl(long double x);

long double fs_expl(long double x)
{
long unsigned n, square;
long double b, e;
static long double x_max, x_min;

if (1 x_max) {
x_max = fs_logl(LDBL_MA X);
x_min = fs_logl(LDBL_MI N);
}
if (x_max >= x && x >= x_min) {
for (square = 0; x 1; x /= 2) {
++square;
}
while (-1 x) {
++square;
x /= 2;
}
e = b = n = 1;
do {
b /= n++;
b *= x;
e += b;
b /= n++;
b *= x;
e += b;
} while (b LDBL_EPSILON / 4);
while (square-- != 0) {
e *= e;
}
} else {
e = x 0 ? LDBL_MAX : 0;
}
return e;
}

long double fs_logl(long double x)
{
long int n;
long double a, b, c, epsilon;
static long double A, B, C;

if (LDBL_MAX >= x && x 0) {
if (1 A) {
A = fs_sqrtl(2);
B = A / 2;
C = fs_logl(A);
}
for (n = 0; x A; x /= 2) {
++n;
}
while (B x) {
--n;
x *= 2;
}
a = (x - 1) / (x + 1);
x = C * n + a;
c = a * a;
n = 1;
epsilon = LDBL_EPSILON * x;
if (0 a) {
if (epsilon 0) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (epsilon b);
} else {
if (0 epsilon) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (b epsilon);
}
x *= 2;
} else {
x = -LDBL_MAX;
}
return x;
}

long double fs_sqrtl(long double x)
{
long int n;
long double a, b;

if (LDBL_MAX >= x && x 0) {
for (n = 0; x 2; x /= 4) {
++n;
}
while (0.5 x) {
--n;
x *= 4;
}
a = x;
b = (1 + x) / 2;
do {
x = b;
b = (a / x + x) / 2;
} while (x b);
while (n 0) {
x *= 2;
--n;
}
while (0 n) {
x /= 2;
++n;
}
} else {
if (x != 0) {
x = LDBL_MAX;
}
}
return x;
}

--
pete

Thanks for sharing this, Pete, but already I have hit a snag.

My configuration of gcc/Cygwin (which includes EVERYTHING since in the
download I didn't know what to exclude) does not define any of those
important long double constants in math.h. And as for float.h, the only
such file I seem to have is the one with the MinGW package that I am
specifically not using here--moreover, it doesn't seem to have any
prototypes or defines I need to refer to. Yes, I know I should do some
research and define those constants myself, since they are no doubt
platform dependent, but my point is that you have clearly composed
these routines in an implementation that may not generalize to other
situations without adjustments being made. This speaks to the very
issue of portability and standardization that is so important to this
group!

Thanks for sharing the routines though. I think I am figuring out that
if I want to compute these elementary math functions to the high degree
of precision I crave, I am going to have to go outside the gcc "box" to
get them, or try my hand at writing them myself!

boy, I have a lot to learn. I appreciate everyone's indulgence.

Les

[Keith Thompson]

"lcw1964" <le***********@ alumni.uwo.cawr ites:
[...]

Actually I spoke too soon!

The following variant of M. Navia's example generates the error too:

#include <stdio.h>
#include <math.h>
int main(void)
{ long double x;
x = 1.0L;
long double n = expl(x);
printf("%.20Lg\ n",n);
return 0;
}

See my other followups in this thread.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[lcw1964]

Keith Thompson wrote:

>
See my other followups in this thread.

--

Thanks, Mr. Thompson--it looks like while you were replying I was
figuring it out on my own concomitantly. Thanks for your consideration.

Les

[Keith Thompson]

"lcw1964" <le***********@ alumni.uwo.cawr ites:
[snip]

My configuration of gcc/Cygwin (which includes EVERYTHING since in the
download I didn't know what to exclude) does not define any of those
important long double constants in math.h. And as for float.h, the only
such file I seem to have is the one with the MinGW package that I am
specifically not using here--moreover, it doesn't seem to have any
prototypes or defines I need to refer to.

The constants LDBL_MAX and so forth are defined in <float.h>, not
<math.h>

Don't waste your time search for a file called "float.h". Just write
code that uses it, and see if it works. For example, this works for
me under Cygwin:

#include <float.h>
#include <stdio.h>
int main(void)
{
long double max = LDBL_MAX;
printf("max = %Lg\n", max);
return 0;
}

The manner in which the compiler finds and processes whatever file or
other entity corresponds to the file.h header is system-specific, and
for the most part you just shouldn't worry about it.

[...]

Thanks for sharing the routines though. I think I am figuring out that
if I want to compute these elementary math functions to the high degree
of precision I crave, I am going to have to go outside the gcc "box" to
get them, or try my hand at writing them myself!

Out of curiousity, why do you want such high precision? Are you sure
that double (as opposed to long double) won't meet your needs?

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[pete]

lcw1964 wrote:

>
pete wrote:

lcw1964 wrote:

expl()

#include <float.h>

long double fs_expl(long double x);
long double fs_logl(long double x);
long double fs_sqrtl(long double x);

long double fs_expl(long double x)
{
long unsigned n, square;
long double b, e;
static long double x_max, x_min;

if (1 x_max) {
x_max = fs_logl(LDBL_MA X);
x_min = fs_logl(LDBL_MI N);
}
if (x_max >= x && x >= x_min) {
for (square = 0; x 1; x /= 2) {
++square;
}
while (-1 x) {
++square;
x /= 2;
}
e = b = n = 1;
do {
b /= n++;
b *= x;
e += b;
b /= n++;
b *= x;
e += b;
} while (b LDBL_EPSILON / 4);
while (square-- != 0) {
e *= e;
}
} else {
e = x 0 ? LDBL_MAX : 0;
}
return e;
}

long double fs_logl(long double x)
{
long int n;
long double a, b, c, epsilon;
static long double A, B, C;

if (LDBL_MAX >= x && x 0) {
if (1 A) {
A = fs_sqrtl(2);
B = A / 2;
C = fs_logl(A);
}
for (n = 0; x A; x /= 2) {
++n;
}
while (B x) {
--n;
x *= 2;
}
a = (x - 1) / (x + 1);
x = C * n + a;
c = a * a;
n = 1;
epsilon = LDBL_EPSILON * x;
if (0 a) {
if (epsilon 0) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (epsilon b);
} else {
if (0 epsilon) {
epsilon = -epsilon;
}
do {
n += 2;
a *= c;
b = a / n;
x += b;
} while (b epsilon);
}
x *= 2;
} else {
x = -LDBL_MAX;
}
return x;
}

long double fs_sqrtl(long double x)
{
long int n;
long double a, b;

if (LDBL_MAX >= x && x 0) {
for (n = 0; x 2; x /= 4) {
++n;
}
while (0.5 x) {
--n;
x *= 4;
}
a = x;
b = (1 + x) / 2;
do {
x = b;
b = (a / x + x) / 2;
} while (x b);
while (n 0) {
x *= 2;
--n;
}
while (0 n) {
x /= 2;
++n;
}
} else {
if (x != 0) {
x = LDBL_MAX;
}
}
return x;
}

--
pete

Thanks for sharing this, Pete, but already I have hit a snag.

My configuration of gcc/Cygwin (which includes EVERYTHING since in the
download I didn't know what to exclude) does not define any of those
important long double constants in math.h.
And as for float.h, the only
such file I seem to have is the one with the MinGW package that I am
specifically not using here--moreover, it doesn't seem to have any
prototypes or defines I need to refer to. Yes, I know I should do some
research and define those constants myself, since they are no doubt
platform dependent, but my point is that you have clearly composed
these routines in an implementation that may not generalize to other
situations without adjustments being made.
This speaks to the very
issue of portability and standardization that is so important to this
group!

The code that I posted, is completely portable C code,
for both hosted and freestanding implementations of C.

--
pete

[lcw1964]

pete wrote:

The code that I posted, is completely portable C code,
for both hosted and freestanding implementations of C.

I am duly humbled. It compiles to an object file just beautifully.

Now, for my next trick, the newbie must learn the command line
parameters to create an executable from more than one object file!

Yes, Google is a great thing. And so are the FAQs, so I think I don't
need any more hand holding. I will take it from here.

As for my answer to Mr. Thompson's question of "why?", I am an
irritating dilettante who can only offer my version of Sir Edmund
Hilary's justification--"because" :)

Les

[Keith Thompson]

"lcw1964" <le***********@ alumni.uwo.cawr ites:
[...]

As for my answer to Mr. Thompson's question of "why?", I am an
irritating dilettante who can only offer my version of Sir Edmund
Hilary's justification--"because" :)

Good enough!

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[lcw1964]

lcw1964 wrote:

pete wrote:

The code that I posted, is completely portable C code,
for both hosted and freestanding implementations of C.

I am duly humbled. It compiles to an object file just beautifully.

Now, for my next trick, the newbie must learn the command line
parameters to create an executable from more than one object file!

The challenge of the multifile project is too daunting to contemplate
at this moment, but I did cut and past Pete's code and prototypes to my
code and put everything in a single file for now and, yes indeed, it
improved my desired results greatly.

The exp() routine does start to lose digits as the argument gets higher
and the output has very large magnitude. I don't know if this due to
propagated rounding error in Pete's code, or if this is just the nature
of the beast when dealing with the long double type on my platform. I
understand that there are computational strategies that one can use to
improve the results of exponential functions as the absolute value of
the argument increases, and I should research this further.

Thank you so much for sharing this and supporting my meandering.

Les

p.s. I know I have been advised not to worry, and I won't, but the
curiosity is killing me--if their is no float.h header file, where the
heck is all that good information, and why doesn't the compiler rebel
when I include an *.h file that really doesn't seem to exist? There is
some zen wisdom in all of this I am sure....