473,320 Members | 1,865 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,320 software developers and data experts.

Portable log1p()

Many C platforms lack the C99 log1p(), which computes
log(1+x) accurately even when x is near 0.

I'm looking for a portable alternative.
Somewhere I found code which can be re-expressed as:

double rep_log1p(double x)
{
return ((1+x)-1) ? log(1+x)*(x/((1+x)-1)) : x;
}

but I do not understand the math, and I am worried of
a "dead spot" somewhere.
I have derived

double my_log1p(double x)
{
return log(1+x)-(((1+x)-1)-x)/(1+x);
}

[the idea here is that if we define
y = 1+x after rounding
e = 1+x-y in the mathematical sense
the derivative of the function log(1+y) is 1/(1+y)
thus log(1+y+e) is close to log(1+y) + e/(1+y) ]
Both functions seem to work, though(my_log1p) seems
somewhat more accurae, and faster.
Any comment on the accuracy or portability ?

In particular, is it guaranteed that a conformant
implementation will NOT optimize things, defeating
the attempt to get accurate results ?

If no, is there a way to insure this ?
François Grieu
Nov 14 '05 #1
3 2935
"Francois Grieu" <fg****@francenet.fr> wrote in message
news:fg**************************@individual.net.. .
Many C platforms lack the C99 log1p(), which computes
log(1+x) accurately even when x is near 0.

I'm looking for a portable alternative.
Somewhere I found code which can be re-expressed as:

double rep_log1p(double x)
{
return ((1+x)-1) ? log(1+x)*(x/((1+x)-1)) : x;
}

but I do not understand the math, and I am worried of
a "dead spot" somewhere.
I ran our log1p tests against this and got:

[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED
log1p: Summary for 1-argument, 53-bit double function:

log1p: Testing 2040 finite function values, result was
log1p: smaller 426 times
log1p: equal 1161 times
log1p: larger 453 times.

log1p: Testing 7 special function values gave 1 mismatches.

log1p: RMS relative bit error was 0.89.
log1p: Max relative bit error was 11.86 for [1001] log1p(+0.000110055)
log1p: Max relative ulp error was 3359 ulp for [1001] log1p(+0.000110055)

64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
5 0 10 21 46 139 234 1161 245 106 37 16 9 6 5
I have derived

double my_log1p(double x)
{
return log(1+x)-(((1+x)-1)-x)/(1+x);
}

[the idea here is that if we define
y = 1+x after rounding
e = 1+x-y in the mathematical sense
the derivative of the function log(1+y) is 1/(1+y)
thus log(1+y+e) is close to log(1+y) + e/(1+y) ]
Both functions seem to work, though(my_log1p) seems
somewhat more accurae, and faster.
For this version, our tests gave:

[2043] log1p(-1): got = -nan, want = -inf -- UNEXPECTED
[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED
log1p: Summary for 1-argument, 53-bit double function:

log1p: Testing 2040 finite function values, result was
log1p: smaller 1 times
log1p: equal 2033 times
log1p: larger 6 times.

log1p: Testing 7 special function values gave 2 mismatches.

log1p: RMS relative bit error was -9.00.
log1p: Max relative bit error was 1.00 for [2001] log1p(-2.22045e-16)
log1p: Max relative ulp error was 2 ulp for [2001] log1p(-2.22045e-16)

64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 2 2033 4 1 0 0 0 0 0
Any comment on the accuracy or portability ?


Our 53-bit version of log1p gave:

log1p: Summary for 1-argument, 53-bit double function:

log1p: Testing 2040 finite function values, result was
log1p: smaller 15 times
log1p: equal 1998 times
log1p: larger 27 times.

log1p: Testing 7 special function values gave 0 mismatches.

log1p: RMS relative bit error was -7.86.
log1p: Max relative bit error was 1.00 for [2001] log1p(-2.22045e-16)
log1p: Max relative ulp error was 2 ulp for [2001] log1p(-2.22045e-16)

64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 16 1998 25 1 0 0 0 0 0

So your function is actually slightly more accurate than ours,
not quite as robust, and probably not quite as fast. (We use
a polynomial approximation over the difficult region.)

BTW, for all tests I used mingw with our log, which is pretty good.
In brief:

log: RMS relative bit error was -9.79.
log: Max relative bit error was 0.90 for [1706] log(+1.70613)
log: Max relative ulp error was 1 ulp for [1048] log(+1.0471)

64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 0 2048 2 0 0 0 0 0 0

HTH,

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com
Nov 14 '05 #2
In article <5L******************@nwrddc02.gnilink.net>,
"P.J. Plauger" <pj*@dinkumware.com> wrote:
ran our log1p tests
Testing agaisnt professional, name brand validation
suite; now THAT'S the power of usenet! Much appreciated!
double my_log1p(double x)
{
return log(1+x)-(((1+x)-1)-x)/(1+x);
}


[2043] log1p(-1): got = -nan, want = -inf -- UNEXPECTED
[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED


Indeed, this should not be.

So your function is actually slightly more accurate than ours,
not quite as robust, and probably not quite as fast. (We use
a polynomial approximation over the difficult region.)


Speed is very dependent on if log() is in hardware or not.
When it is, my_log1p() is hard to beat in the 0.002 to 0.5
and and -0.25 to -0.002 ranges, I guess.
When log() is in SW, polynomial approximation seems the way
to go.
Meanwhile, a private email pointed out that the above math
is close to the one used by GSL (GNU Scientific Library),

#include <GNU_GPL_v2.h>
double gsl_log1p (const double x)
{
volatile double y;
y = 1 + x;
return log(y)-((y-1)-x)/y ; /* cancels errors with IEEE arithmetic */
}

which however has an interesting extra twist: the use of
"volatile" prevents - I believe/hope somewhat reliably -
the compiler from "simplifying" ((1+x)-1)-x, or performing
similar damage.

I'm still trying to find what the C standards says about how
an implementation is allowed to reshuffle expressions like
log(1+x)-(((1+x)-1)-x)/(1+x)

In also worried of implementations where computations
are carried on e.g 80 bits and rounded to 64 bits when stored.
François Grieu
Nov 14 '05 #3
On Sat, 08 May 2004 12:03:24 +0200, Francois Grieu <fg****@francenet.fr>
wrote:
I'm still trying to find what the C standards says about how
an implementation is allowed to reshuffle expressions like
log(1+x)-(((1+x)-1)-x)/(1+x)


If nothing else, the "as-if" rule applies: it can do anything it wants as
long as the results are identical - as if it made no change.
--
#include <standard.disclaimer>
_
Kevin D Quitt USA 91387-4454 96.37% of all statistics are made up
Per the FCA, this address may not be added to any commercial mail list
Nov 14 '05 #4

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

13
by: James Harris | last post by:
Hi, Can someone recommend a book that will teach me how to approach C programming so that code is modularised, will compile for different environments (such as variations of Unix and Windows),...
2
by: Francois Grieu | last post by:
Many C platforms lack the C99 log1p(), which computes log(1+x) accurately even when x is near 0. I'm looking for a portable (C) alternative. Somewhere I found code which can be re-expressed...
10
by: Jason Curl | last post by:
Dear C group, I'm very interested in writing portable C, but I only have GNU, Sparc and Cygwin to compile on. What I find is the biggest problem to writing portable C is what headers to...
8
by: suresh_C# | last post by:
Dear All, What is difference between Portable Executable (PE) file and a Assembly? Thanks, Mahesh
2
by: Tull Clancey | last post by:
Hi all. I'm nearing completion of a host app that needs to send data to, and receive data back from a portable, an HP device. The application running on the portable will be a .net application....
131
by: pemo | last post by:
Is C really portable? And, apologies, but this is possibly a little OT? In c.l.c we often see 'not portable' comments, but I wonder just how portable C apps really are. I don't write...
409
by: jacob navia | last post by:
I am trying to compile as much code in 64 bit mode as possible to test the 64 bit version of lcc-win. The problem appears now that size_t is now 64 bits. Fine. It has to be since there are...
10
by: bramnizzle | last post by:
I don't know if this is the right thread or not, but... In my endless pursuit to hold on to my dinosaur laptop (Dell 1100 Inspiron - circa 2004)...I want to keep as much space free for my...
23
by: asit | last post by:
what is the difference between portable C, posix C and windows C ???
0
by: DolphinDB | last post by:
Tired of spending countless mintues downsampling your data? Look no further! In this article, you’ll learn how to efficiently downsample 6.48 billion high-frequency records to 61 million...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
1
by: PapaRatzi | last post by:
Hello, I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...
0
by: CloudSolutions | last post by:
Introduction: For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
0
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.