27 ms·

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 expans

by cb321 2y ago

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.) [1] https://github.com/zachartrand/SoME-3-Living/blob/main/scripts/log_fast.py 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).

pclmulqdq 2y ago
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.
njuffa 2y ago
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; }
- cb321 2y ago
  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 :-) }: $ python3 j.py | headtail 3 3 # [2/3..4/3], m=300 0.6666666666666666 3.7747748309890916e-06 0.6688888888888889 3.5668329603222304e-06 0.6711111111111111 3.3693446750415745e-06 --- 1.3266666666666544 -3.120628879194065e-07 1.3288888888888766 -3.2510972441945896e-07 1.3311111111110987 -3.3859621345611046e-07 $ python3 k.py | headtail 3 3 # [rt1/2..rt2], m=300 0.7071067811865476 1.279804798859363e-06 0.7094638037905028 1.19727672415193e-06 0.711820826394458 1.1195756742798935e-06 --- 1.4071424945612245 -1.1572535666970474e-06 1.4094995171651796 -1.1970089550117358e-06 1.4118565397691347 -1.2378520898859158e-06 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.

by cb321 2y ago

pclmulqdq 2y ago

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.

njuffa 2y ago

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; }

cb321 2y ago
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 :-) }: $ python3 j.py | headtail 3 3 # [2/3..4/3], m=300 0.6666666666666666 3.7747748309890916e-06 0.6688888888888889 3.5668329603222304e-06 0.6711111111111111 3.3693446750415745e-06 --- 1.3266666666666544 -3.120628879194065e-07 1.3288888888888766 -3.2510972441945896e-07 1.3311111111110987 -3.3859621345611046e-07 $ python3 k.py | headtail 3 3 # [rt1/2..rt2], m=300 0.7071067811865476 1.279804798859363e-06 0.7094638037905028 1.19727672415193e-06 0.711820826394458 1.1195756742798935e-06 --- 1.4071424945612245 -1.1572535666970474e-06 1.4094995171651796 -1.1970089550117358e-06 1.4118565397691347 -1.2378520898859158e-06 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.