The article's impl of its own series [1] centers on the octave [2/3, 4/3] but the error is actually more balanced across [1/sqrt(2), sqrt(2)]. (A series expansion of the error probably explains this - some of the experts here like pavpanchekha | pclmulqdq might be able to explain pithily.)
So, here is one way to evaluate the discussed series evaluated across the centered interval and I think it only needs 9 terms for "near IEEE double" NOT the author's 15 terms:
import sys, math # arg = number of terms of series
n = int(sys.argv[1]) if len(sys.argv) > 1 else 3
def lnATH(x): # Ln via ArcTanH form/series
u = (x - 1)/(x + 1) # Can expand loop to expr
return sum(2/(2*j+1)*u**(2*j+1) for j in range(n))
x = math.sqrt(0.5) # sqrt(1/2) centers error
x1 = x*2 # IEEE FP exponent gets within 1 octave
m = 3840 # 4K monitor; adjust to whatever
dx = (x1 - x)/m
for i in range(m): # Plot error in favorite tool
print(x, lnATH(x) - math.log(x))
x += dx
If you run it with 8 you begin to see round-off effects. With 9 terms, it becomes dominated by such effects.
Anyway 15/9 = 1.66X speed cost which seemed enough to be noteworthy. I mean, he does call it "log_fast.py" after all. (4 seems adequate for IEEE single precision - though you might be 1-bit shy of 24-bits of relative precision at the very edges of octaves if that matters to you.)
I'm just guessing here. What stands out is that the extrema of that octave have the same value. I would guess that this means that things are probably especially well-behaved in that region.
If you look at the Chebyshev series of log in that range, you may see some very small coefficients.
Re: centers on the octave [2/3, 4/3] but the error is actually more balanced across [1/sqrt(2), sqrt(2)]
For the computation of the natural logarithm using binary floating-point arithmetic there is no practical difference, i.e. the overall accuracy that can be achieved is basically the same. Here is a worked example based on minimax core approximations extracted from my personal code collection (design and implementation are my own; no guarantee of fitness for any particular purpose).
/* natural logarithm */
float my_logf (float a)
{
float m, r, s, t, i, f;
int32_t e;
#if LOG_VARIANT == 1
// max ulp err = 0.86280
const float cutoff_f = 0.707106781f;
#elif LOG_VARIANT == 3
// max ulp err = 0.85089
const float cutoff_f = 0.666666667f;
#endif // LOG_VARIANT
if ((a > 0.0f) && (a <= 0x1.fffffep+127f)) { // 3.40282347e+38f
m = frexpf (a, &e);
if (m < cutoff_f) {
m = m + m;
e = e - 1;
}
i = (float)e;
f = m - 1.0f;
s = f * f;
#if LOG_VARIANT == 1
/* Compute log1p(m) for m in [sqrt(0.5)-1, sqrt(2.0)-1] */
r = -0x1.384000p-4f; // -0.076232910
t = 0x1.084000p-3f; // 0.129028320
r = fmaf (r, s, -0x1.0f218cp-3f); // -0.132388204
t = fmaf (t, s, 0x1.226b04p-3f); // 0.141805679
r = fmaf (r, s, -0x1.5427a4p-3f); // -0.166091233
t = fmaf (t, s, 0x1.99a35ap-3f); // 0.200018600
r = fmaf (r, s, -0x1.000416p-2f); // -0.250015587
r = fmaf (t, f, r);
r = fmaf (r, f, 0x1.55555cp-2f); // 0.333333433
r = fmaf (r, f, -0x1.fffffap-2f); // -0.499999911
#elif LOG_VARIANT == 3
/* Compute log1p(f) for f in [-1/3, 1/3] */
r = -0x1.0ae000p-3f; // -0.130310059
t = 0x1.208000p-3f; // 0.140869141
r = fmaf (r, s, -0x1.f1988ap-4f); // -0.121483363
t = fmaf (t, s, 0x1.1e5740p-3f); // 0.139814854
r = fmaf (r, s, -0x1.55b36ep-3f); // -0.166846141
t = fmaf (t, s, 0x1.99d8b2p-3f); // 0.200120345
r = fmaf (r, s, -0x1.fffe02p-3f); // -0.249996200
r = fmaf (t, f, r);
r = fmaf (r, f, 0x1.5554fap-2f); // 0.333331972
r = fmaf (r, f, -0x1.000000p-1f); // -0.500000000
#endif // LOG_VARIANT
r = fmaf (r, s, f);
r = fmaf (i, 0x1.62e430p-01f, r); // 0.693147182 // log(2)
} else {
r = a + a; // silence NaNs if necessary
if (a < 0.0f) r = 0.0f/0.0f; // QNaN INDEFINITE
if (a == 0.0f) r = -INFINITY; // -INF
}
return r;
}
Thanks for your code post, but your fused-multiply-add Horner coeffs seem like they are not from the same power series I was discussing / the article was promoting.
The effects I mention are pretty clear if you run my code with both `x = 2/3.0` and with `x = math.sqrt(0.5)`, even with a much smaller `m` / coarser granularity.
For example { for the python &| copy-paste impaired &| lazy :-) }:
Note in the 2/3..4/3 case the exponent changes from e-06 to e-07 - a whole order of magnitude difference in |error|, while in the rt1/2..rt2 case the nearness of |error|. This same pattern manifests for higher numbers of terms like 8 or 9 or higher resolution like m=3840.
So, here is one way to evaluate the discussed series evaluated across the centered interval and I think it only needs 9 terms for "near IEEE double" NOT the author's 15 terms:
If you run it with 8 you begin to see round-off effects. With 9 terms, it becomes dominated by such effects.Anyway 15/9 = 1.66X speed cost which seemed enough to be noteworthy. I mean, he does call it "log_fast.py" after all. (4 seems adequate for IEEE single precision - though you might be 1-bit shy of 24-bits of relative precision at the very edges of octaves if that matters to you.)
[1] https://github.com/zachartrand/SoME-3-Living/blob/main/scrip... (and yeah, the code says "16 terms" & article "15 terms", but any way you slice it, I think it's a lot of extra terms).