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