Hmmm ... Maybe I'm doing something wrong, but my installation of GPC (DJGPP) doesn't know what "Log1P" is, nor can I find it in the source code of the GPC unit, etc. Anyway, that's my problem, so I'll look into that further when I can.
As for my comment about "tweaked" coefficients (for platforms which don't already provide the function), the general rule is that if you truncate an infinite series, you need to "tweak" the coefficient values to provide optimum accuracy for the number of terms you are using. Simply truncating a series may produce inaccurate results, regardless of how many terms you apply. That is the general rule, and I remember it was certainly true with trigonometric functions. However, perhaps this effect is not so prevalent with the "log1p" coefficients? If so, that's good - it means we don't have to go to the trouble of finding "tweaked/optimized" coefficients (or find out how to derive them). Unfortunately, I don't remember the name of the book in which I found the "tweaked" coefficients for trigonometric functions, nor whether this also had similar coefficients for other functions (it was quite a few years ago now ;-).
Joe.
-----Original Message----- From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] Sent: Saturday, February 09, 2002 12:47 AM To: gpc@gnu.de Subject: Re: Complex_Arctan
Joe wrote:
BTW, AFAICT, DJGPP does not have a Log1P function, so I wasn't able to verify the above anyway. I didn't know the values to use for the coefficients "Log1PCoefs[]", so I couldn't try the suggested Log1P code, either. Although
??? You don't have to guess the coefficients, they are initialized in the code. (If you don't want to rebuild the RTS just to test it (I didn't either), move the "to begin do" stuff into your main program, and import the GPC unit (because of the SplitReal function).)
Emil presented a power series which I assume represents these coefficients, I don't know whether these need to be "tweaked". If I understand this correctly, it is an infinite series. However, AFAIK, when you apply a finite version of such a series, you need to "tweak" it's values to get good accuracy (no, I don't know how to do this).
I remember experimenting years ago with the power series for "sine" or "cosine" and being amazed at how a small finite "tweaked" power series was much more accurate than the "conventional" power series, no matter how many terms I calculated for it (after a while, the more terms you calculate, the greater the error you introduce). I presume the "log1p" power series also has this effect(?)
Not sure what sort of "tweaking" you have in mind. Anyway, my Log1P implementation tries to use as few terms of the power series as possible, and it sums them in reverse order. I don't know if this is optimal, but it appeared to work with acceptable precision: the largest error I found was 2 ulp's when compiled without optimization, and 1 ulp when compiled with -O2, both on an i586.
Emil
da Silva, Joe wrote:
Hmmm ... Maybe I'm doing something wrong, but my installation of GPC (DJGPP) doesn't know what "Log1P" is, nor can I find it in the source code of the GPC unit, etc. Anyway, that's my problem, so I'll look into that further when I can.
Emil posted a Pascal version of Log1P recently. And if the libc contains it (DJGPP's does), an external declaration will do. However, the next GPC release (there hasn't been any since this discussion) will have it built-in (though I'm calling it Ln1Plus rather than Log1P then because this seems more consistent with the Pascal naming, i.e. Ln for the natural logarithm vs. log() in C).
As for my comment about "tweaked" coefficients (for platforms which don't already provide the function), the general rule is that if you truncate an infinite series, you need to "tweak" the coefficient values to provide optimum accuracy for the number of terms you are using. Simply truncating a series may produce inaccurate results, regardless of how many terms you apply. That is the general rule, and I remember it was certainly true with trigonometric functions. However, perhaps this effect is not so prevalent with the "log1p" coefficients?
As I said, I'm no expert on numerics, and I haven't heard of this tweaking.
Anyway, what Emil's code does is to truncate the series at the point where the addends are smaller than 1 ulp so adding them would not change the result any further. The number of terms necessary for this is computed in advance and rounded up to be on the safe side.
While it is possible that several small addends could sum up to be more than 1 ulp, I don't think they could exceed 2 ulp (as Emil noted) for a quadratically convergent series. AFAIUI, that's why he used ln (1+x) = ln ((1+u) / (1-u)) = ln (1+u) - ln (1-u) with u = x/(2+x) because in the series of ln (1+u) and ln (1-u) the even terms cancel out, so the subseries consisting of the odd terms is quadratically convergent because |u| < 1/3.
Frank