long double versions of functions in gcc under Cygwin

[Keith Thompson]

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]

I'm in the Knuth camp:
http://portal.acm.org/citation.cfm?id=356640&dl=&coll=GUIDE&CFID=15151515&CFTOKEN=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.)

 

[lcw1964]

#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_MAX);
x_min = fs_logl(LDBL_MIN);
}
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;
}

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]

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.

 

[lcw1964]

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]

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?

 

[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.

 

[lcw1964]

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]

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!

 

[lcw1964]

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....

 

[Keith Thompson]

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....

<OT>
You can ask gcc to tell you where it searchs for headers. Use the
"-v" option. Don't assume that the float.h that it uses is under
/usr/include.
</OT>

More generally, system headers aren't necessarily implemented as
files; they can be implemented magically in the compiler itself, or by
some sort of precompiled thingie that isn't a C source file. <OT>As
far as I know, gcc doesn't do this.</OT>

 

[Chris Torek]

More generally, system headers aren't necessarily implemented as
files; they can be implemented magically in the compiler itself, or by
some sort of precompiled thingie that isn't a C source file. <OT>As
far as I know, gcc doesn't do this.</OT>

Regarding the OT bit: no, but it does something just as sneaky.
The compiler pre-defines a bunch of underscore-prefixed names that
can then be used in the actual headers (which are often squirreled
away all over the file system) to decide how the compiler was
invoked. This can in turn change the way various types work,
including wchar_t. (Note that "wchar_t" has to match the compiler's
type for L"string"s, so this is actually quite a reasonable thing
to do in the first place.)

 

[ena8t8si]

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.

I don't agree with the basic premise, but for the
sake of discussion let's say I do. Do you really
think most gotos would be eliminated by C having
the features you named? In my experience most gotos
arise in one of three circumstances: bad coding,
branch to end of function to preserve single return
point, and machine generated C code (eg, YACC/LEX).
Breaking out of multiple loops seems like a distant
fourth (and usually remediable by putting the multiple
loops in their own function and using return); and
exception handling, it's hard to see how that's even
on the radar. It may be that C would benefit from
these features (personally I don't think it would,
but that's a separate issue), but even if it had
them it doesn't seem like goto usage would be affected
much.

My question above wasn't meant to be rhetorical;
I'm interested to hear your reactions.

 

[Keith Thompson]

I don't agree with the basic premise, but for the
sake of discussion let's say I do. Do you really
think most gotos would be eliminated by C having
the features you named? In my experience most gotos
arise in one of three circumstances: bad coding,
branch to end of function to preserve single return
point, and machine generated C code (eg, YACC/LEX).
Breaking out of multiple loops seems like a distant
fourth (and usually remediable by putting the multiple
loops in their own function and using return); and
exception handling, it's hard to see how that's even
on the radar. It may be that C would benefit from
these features (personally I don't think it would,
but that's a separate issue), but even if it had
them it doesn't seem like goto usage would be affected
much.

My question above wasn't meant to be rhetorical;
I'm interested to hear your reactions.

Bad coding obviously won't be affected. Yacc and Lex raise another
issue, which I'll come back to.

Branching to the end of a function (for example because you need to
execute some cleanup code rather than doing an immediate return) is
something that's usually done in, ahem, *exceptional* circumstances.
I'm not sure that C++-style exceptions would be a good fit for C. As
I recall, any C++ object can be thrown as an exception. Ada's
exception mechanism might be worth considering. In the 1983 version
of the language, it's simpler than C++'s mechanism; an exception is
just an exception, and a handler can handle either named exception or
all exceptions. Of course, any such mechanism would have to be very
lightweight.

The idea is to work with named entities that relate to the problem
domain rather than jumping to a named point in the code.

Getting back to Yacc and Lex, I haven't taken much of a look at the
code they generate. But a finite state machine is one of the few
places where gotos make sense. Ok, that was too strong; it's a case
where they make more sense than in most places.

Most control constructs correspond to something in the problem domain:
if/then/else makes a decision, a loop iterates over some set of
thingies, a function performs some part of the problem, etc. But a
goto just operates on the control flow of your program.

The one time when a goto statement matches something in the problem
domain is in a finite state machine, where a goto represents a
transition to a specified state (and the state is represented by a
chunk of code).

On the other hand, a finite state machine can also be implemented as a
case statement in a loop, with the state being represented by the
value of a variable.

Anyway, it's not terribly likely that any such changes are going to be
made.

 

[ena8t8si]

Branching to the end of a function (for example because you need to
execute some cleanup code rather than doing an immediate return) is
something that's usually done in, ahem, *exceptional* circumstances.
I'm not sure that C++-style exceptions would be a good fit for C. As
I recall, any C++ object can be thrown as an exception. Ada's
exception mechanism might be worth considering. In the 1983 version
of the language, it's simpler than C++'s mechanism; an exception is
just an exception, and a handler can handle either named exception or
all exceptions. Of course, any such mechanism would have to be very
lightweight.

The idea is to work with named entities that relate to the problem
domain rather than jumping to a named point in the code.

The pattern I was talking about looks something like
the following:

int typfun()
{
int rc;
...
x = get_resource_one();
if(x==NULL){ rc = -1; goto BAIL1; }

y = get_resource_two();
if(y==NULL){ rc = -2; goto BAIL2; }

z = get_resource_three();
if(z=NULL){ rc = -3; goto BAIL3; }

rc = 0;
release_resource(z);

BAIL3:
release_resource(y);

BAIL2:
release_resource(x);

BAIL3:
return rc;
}

Obviously there are variations, and I've written
the if()'s on one line to save space, but I
think you get the idea.

A pattern like this doesn't map very well onto
try/catch blocks, or other exception handling
mechanisms. I've used languages with exception
handling built in going back to the 1970's and
I think I'm pretty familiar with the landscape;
for expressing this pattern of control structure
(and this pattern is one of the common ones for
goto), just using goto allows a simpler and more
more clear function structure than try/catch
blocks. Or would you say otherwise?

 

[pete]

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.

http://groups.google.com/group/comp.lang.c/msg/e7cff29955d4971b

"We used to do successive multiplication, until we discovered
that accuracy degrades quickly with that approach.
The only really safe way to evaluate pow(x, y) accurately
is to compute exp(y * ln(x)) to extra internal precision.
(Took us a lot of rewrites, and a lot of careful testing,
to find that out.)" -- P.J. Plauger, Dinkumware, Ltd.

 

[jaysome]

The pattern I was talking about looks something like
the following:

int typfun()
{
int rc;
...
x = get_resource_one();
if(x==NULL){ rc = -1; goto BAIL1; }

y = get_resource_two();
if(y==NULL){ rc = -2; goto BAIL2; }

z = get_resource_three();
if(z=NULL){ rc = -3; goto BAIL3; }

rc = 0;
release_resource(z);

BAIL3:
release_resource(y);

BAIL2:
release_resource(x);

BAIL3:
return rc;
}

Obviously there are variations, and I've written
the if()'s on one line to save space, but I
think you get the idea.

A pattern like this doesn't map very well onto
try/catch blocks, or other exception handling
mechanisms.

This pattern maps very well onto try/catch blocks. First consider a
re-write of your example code (z=NULL corrected).

int typfun(void)
{
int rc;
...
x=y=z=NULL;
x = get_resource_one();
if(x==NULL){ rc = -1; goto CLEANUP; }
y = get_resource_two();
if(y==NULL){ rc = -2; goto CLEANUP; }
z = get_resource_three();
if(z==NULL){ rc = -3; goto CLEANUP; }

rc = 0;

CLEANUP:
release_resource(x);
release_resource(y);
release_resource(z);

return rc;
}

Now consider how it might be written in a language that provides
try/catch, and assumes release_resource() handles NULL, like it
should, like free() and delete and delete [] do:

int typfun()
{
int rc = 0;
...
x=y=z=NULL;
try
{
x = get_resource_one();
if(x==NULL){ rc = -1; throw rc; }
y = get_resource_two();
if(y==NULL){ rc = -2; throw rc; }
z = get_resource_three();
if(z==NULL){ rc = -3; throw rc; }
}
catch(int i)
{
rc = i;
}

release_resource(x);
release_resource(y);
release_resource(z);

return rc;
}

There's not much difference.

 

[Keith Thompson]

The pattern I was talking about looks something like
the following:

int typfun()
{
int rc;
...
x = get_resource_one();
if(x==NULL){ rc = -1; goto BAIL1; }

y = get_resource_two();
if(y==NULL){ rc = -2; goto BAIL2; }

z = get_resource_three();
if(z=NULL){ rc = -3; goto BAIL3; }

rc = 0;
release_resource(z);

BAIL3:
release_resource(y);

BAIL2:
release_resource(x);

BAIL3:
return rc;
}

Obviously there are variations, and I've written
the if()'s on one line to save space, but I
think you get the idea.

Sure. (The last label should be BAIL1, not BAIL3).

Just off the top of my head ...

With exception handlers, that would probably be written as three
nested blocks rather than straight-line code with gotos. As long as
you're creating nested blocks, I suppose you might as well just use
if-then-elses.

The code could be a bit simpler if you could assume that
release_resource() does nothing with a null pointer, as free() does.

Something that just occurred to me is that exceptions (in Ada or C++)
propagate across function calls. If we're (hypothetically) doing this
just to handle control flow within functions, that's not necessary.
Perhaps one answer is, rather than exception handling, to allow a
named break to exit from any scope, not just a loop or switch
statement.

So maybe something like this:

int typfun() /* This is not C. */
{
int rc;
...
USE_RESOURCE_ONE:
{
x = get_resource_one();
if(x==NULL){ rc = -1; break USE_RESOURCE_ONE };

USE_RESOURCE_TWO:
{
y = get_resource_two();
if(y==NULL){ rc = -2; break USE_RESOURCE_TWO; }

USE_RESOURCE_THREE:
{
z = get_resource_three();
if(z=NULL){ rc = -3; break USE_RESOURCE_THREE; }

rc = 0;
release_resource(z);
}

release_resource(y);
}

release_resource(x);
}

return rc;
}

Or, if release_resource() works correctly with a null pointer argument:

int typfun() /* This is not C. */
{
int rc;
...
USE_RESOURCES:
{
x = get_resource_one();
if(x==NULL){ rc = -1; break USE_RESOURCES; }

y = get_resource_two();
if(y==NULL){ rc = -2; break USE_RESOURCES; }

z = get_resource_three();
if(z=NULL){ rc = -3; break USE_RESOURCES; }

rc = 0;
} /* end of USE_RESOURCES block */

release_resource(z);
release_resource(y);
release_resource(x);

return rc;
}

 

[Richard Bos]

The pattern I was talking about looks something like
the following:

int typfun()
{
int rc;
...
x = get_resource_one();
if(x==NULL){ rc = -1; goto BAIL1; }

y = get_resource_two();
if(y==NULL){ rc = -2; goto BAIL2; }

z = get_resource_three();
if(z=NULL){ rc = -3; goto BAIL3; }

rc = 0;
release_resource(z);

BAIL3:
release_resource(y);

BAIL2:
release_resource(x);

BAIL3:
return rc;
}

Obviously there are variations, and I've written
the if()'s on one line to save space, but I
think you get the idea.

A pattern like this doesn't map very well onto
try/catch blocks, or other exception handling
mechanisms.

It works better if you manage to write release_resource() to do nothing
on a null pointer (like malloc(), and unfortunately unlike fclose()).
You can then write this, which IMO looks cleaner:

int typfun()
{
int rc;
...
x = get_resource_one();
if(x==NULL){ rc = -1; goto BAIL; }

y = get_resource_two();
if(y==NULL){ rc = -2; goto BAIL; }

z = get_resource_three();
if(z==NULL){ rc = -3; goto BAIL; }

rc = 0;

BAIL:
release_resource(z);
release_resource(y);
release_resource(x);
return rc;
}

That is probably more amenable to being converted to a try/catch
mechanism, too - though I still wouldn't choose that conversion myself.

Richard

 

[ena8t8si]

This pattern maps very well onto try/catch blocks. First consider a
re-write of your example code (z=NULL corrected).

Yes thank you for the correction.

int typfun(void)
{
int rc;
...
x=y=z=NULL;
x = get_resource_one();
if(x==NULL){ rc = -1; goto CLEANUP; }
y = get_resource_two();
if(y==NULL){ rc = -2; goto CLEANUP; }
z = get_resource_three();
if(z==NULL){ rc = -3; goto CLEANUP; }

rc = 0;

CLEANUP:
release_resource(x);
release_resource(y);
release_resource(z);

return rc;
}

Now consider how it might be written in a language that provides
try/catch, and assumes release_resource() handles NULL, like it
should, like free() and delete and delete [] do:

int typfun()
{
int rc = 0;
...
x=y=z=NULL;
try
{
x = get_resource_one();
if(x==NULL){ rc = -1; throw rc; }
y = get_resource_two();
if(y==NULL){ rc = -2; throw rc; }
z = get_resource_three();
if(z==NULL){ rc = -3; throw rc; }
}
catch(int i)
{
rc = i;
}

release_resource(x);
release_resource(y);
release_resource(z);

return rc;
}

There's not much difference.

Two comments. First, not all resource allocation/release can be
rewritten in the style you suggest using release(NULL). [*]
Second, why would someone use throw when what's it's doing
basically just a goto in disguise? If the release(NULL) technique
works, it seems better to dispense with try/catch altogether:

int typfun()
{
int rc = 0;
...
x=y=z=NULL;
do {
x = get_resource_one();
if(x==NULL){ rc = -1; break; }
y = get_resource_two();
if(y==NULL){ rc = -2; break; }
z = get_resource_three();
if(z==NULL){ rc = -3; break; }
} while(0);

release_resource(x);
release_resource(y);
release_resource(z);

return rc;
}

[*] For performance reasons, among others.