Estimate logarithm quickly

08212021, 03:39 PM
Post: #1




Estimate logarithm quickly
Estimate logorithm quickly, with tiny errors.
source: Dead Reckoning: Calculating without Instruments, p. 123 n*n = (n+x) * (nx) * (n*n/(n*nx*x)) 2*ln(n) = ln(n+x) + ln(nx) + ln(n*n/(n*nx*x)) ln(n*n/(n*nx*x)) = ln(1x*x/(n*nx*x)) = 2*atanh(x*x/(2*n*nx*x)) Let a = ln(n+x), b = ln(nx), c = n/x, and halved both side. ln(n) = (a+b)/2 + atanh(1/(2*c*c1)) atanh(1/(2*c*c1)) > atanh(1/c) / (2c) = (ab)/2 / (2c) → ln(n) > (a+b)/2 + (ab)/2 / (2c) If we divide both side by ln(base), inequality still holds We can thus generalize for log of any base. Example, n=3, x=1 a = log(3+1) = 2*log(2) b = log(31) = log(2) c = 3/1 = 3 log(3) > 3*log(2)/2 + (log(2)/2) / (2*3) = 19/12 * log(2) log(3)/log(2) > 19/12 Example, n=3, x=1.5 a = log(3+1.5) = 2*log(3)log(2) b = log(31.5) = log(3)log(2) c = 3/1.5 = 2 log(3) > (3*log(3)2*log(2))/2 + (log(3))/2 / (2*2) = log(2) + 13/8*log(3) log(2) > 5/8*log(3) log(3)/log(2) < 8/5 19/12 < log2(3) < 8/5 Both extreme happened to be convergents of log2(3) C:\> spigot C log(3)/log(2)  head 1/1 2/1 3/2 8/5 19/12 65/41 84/53 485/306 1054/665 24727/15601 

08212021, 03:58 PM
Post: #2




RE: Estimate logarithm quickly
(08212021 03:39 PM)Albert Chan Wrote: atanh(1/(2*c*c1)) > atanh(1/c) / (2c) = (ab)/2 / (2c) This is not yet proofed (nx > 0, c = n/x > 1) To proof it, it required this inequality (also not proofed, suggestions welcome): (1w)^(2w) * (1+w)^(2+w) < 1, where 0 < w < 1 XCAS> series((1w)^(2w) * (1+w)^(2+w), w, 0, 10) 1  w^4/3  4*w^6/15  10*w^8/63  4*w^10/45 + w^11*order_size(w) Taylor series suggested the expression peaked at w = 0, and plots support it. For now, assuming it is true ... (1w)^(2w) * (1+w)^(2+w) < 1 (1w*w)^2 < ((1w)/(1+w))^w (1/(1w*w))^2 > ((1+w)/(1w))^w Let c = 1/w, c > 1 (c*c/(c*c1))^2 > ((1+1/c)/(11/c)) ^ (1/c) ln(c*c/(c*c1)) * 2 > ln((1+1/c)/(11/c)) / c Simpilfy with identity ln(z) = 2*atanh((z1)/(z+1)) (2*atanh(1/(2*c*c1)) * 2 > (2*atanh(1/c)) / c atanh(1/(2*c*c1) > atanh(1/c) / (2c) 

08212021, 05:41 PM
Post: #3




RE: Estimate logarithm quickly
(08212021 03:58 PM)Albert Chan Wrote: To proof it, it required this inequality (also not proofed, suggestions welcome): My attempt to proof this, by showing there is no local extremum for 0 < w < 1 To get the extrema of above is messy. Instead, we solve extrema of ln(expression) XCAS> f := (2w)*ln(1w) + (2+w)*ln(1+w) XCAS> diff(f,w) ln(1w)  (2w)/(1w) + ln(1+w) + (2+w)/(1+w) XCAS> numer(ans()) // local extrema if df/dw = 0 ln(w+1) + ln(w+1) + w^2*ln(w+1)  w^2*ln(w+1) + 2*w Divide all by 2, then replace with identity atanh(x) = ln((x+1)/(1x))/2 (w^21) * atanh(w) + w = (w^21) * (w + w^3/3 + w^5/5 + ...) + w All coefficients of w's is positive, 1/(2k1)  1/(2k+1) = 2/(4k^21) > 0 XCAS> series((w^21)*atanh(w) + w, w,0,10) 2*w^3/3 + 2*w^5/15 + 2*w^7/35 + 2*w^9/63 + w^11*order_size(w) Within 0 < w < 1, there is no local extremum, f is bounded by edges. XCAS> limit(f, w=0, +1) → 0 XCAS> limif(f, w=1, 1) → infinity f = ln((1w)^(2w) * (1+w)^(2+w)) < 0 → (1w)^(2w) * (1+w)^(2+w) < 1 

08212021, 06:31 PM
Post: #4




RE: Estimate logarithm quickly
Albert, thanks for the book reference.
(Surprisingly, that's one I don't have...yet.) ENTER > = 

08212021, 11:42 PM
Post: #5




RE: Estimate logarithm quickly
This is an amazing formula, for estimating ln(n) to 5 sig. digits, n = [1/√2, √2]
(again, from "Dead Reckoning", p. 132) ln(n) ≈ 6*(n1) / (4*√n + n+1) use 1/ln(10) = 0.4342945 ... to scale back for log10. Example, for log10(3): 3^2 = 0.9e1 log10(0.9) ≈ 0.4343 * 6 * .1 / (4*√.9 + 1.9) ≈ 0.04576 log10(3) ≈ (0.04576 + 1) / 2 = 0.47712 Or, scale by 2^k factor, log10(2) ≈ 0.30103 3 * 4 = 1.2e1 log10(1.2) = 0.4343 * 6*.2/(4*√1.2 + 2.2) ≈ 0.07918 log10(3) ≈ 0.07918 + 1  2*0.30103 = 0.47712 

08232021, 06:44 AM
Post: #6




RE: Estimate logarithm quickly
(08212021 03:39 PM)Albert Chan Wrote: Estimate logorithm quickly, with tiny errors. Just to note, this is Simon Tatham's arbitraryprecision spigot calculator, which deals in continued fractions as well as decimals: https://www.chiark.greenend.org.uk/~sgta...html#cfrac 

10062021, 10:58 PM
Post: #7




RE: Estimate logarithm quickly
(10062021 05:12 PM)Eric Rechlin Wrote: Computing Logarithms by Integration  Richard Schwartz (17:02) Thanks for the video. I finally "get" Doerfler's formula ln(n) = (n1) / ∫(n^t, t, 0, 1) ≈ (n1) / ((1 + 4*√n + n)/6) // Simpson's Rule. I also learned about Boole's Rule. It is really just Richardson extrapolation, from Simpson's Rule. CAS> S1 := [1,0,4,0,1] ./ 6 CAS> S2 := [1,4,2,4,1] ./ 12 CAS> B1 := S2 + (S2S1)/(4^21) → [7/90, 16/45, 2/15, 16/45, 7/90] CAS> horner(B1, e^0.25) / (e1) → 1.00000050019 // estimated 1/ln(e), Boole's Rule This should match extrapolation of 2 Simpson's. Note: we work with reciprocal log, to simplify Richardson's extrapolation. CAS> D(x) := horner([1,4,1], x^0.5)/(6*(x1)) // Doerfler's formula, for 1/ln(x) CAS> s1 := D(e) → 1.00033715273 CAS> s2 := D(e^0.5)/2 → 1.00002154097 CAS> horner([1,4,2,4,1], e^0.25)/(12*(e1)) → 1.00002154097 Note: s2 is equivalent to more messy [1,4,2,4,1] Simpson's. version. Thus, it is a good idea to reduce log arguments, close to 1 Extrapolate for Boole's Rule result: CAS> s2 + (s2s1)/(4^21) → 1.00000050019 

10072021, 06:25 PM
Post: #8




RE: Estimate logarithm quickly
I have a very easy to remember method if I want to sketch log/log or log/lin data with hand:
log2 > 0.3 log3 > 0.5 log5 > 0.7 These values approx equally spaced on the log scale and with log1 and log10 values there are 5 points on one decade. Csaba 

10182021, 01:02 PM
(This post was last modified: 07282022 04:32 PM by Albert Chan.)
Post: #9




RE: Estimate logarithm quickly
(10062021 10:58 PM)Albert Chan Wrote:(10062021 05:12 PM)Eric Rechlin Wrote: Computing Logarithms by Integration  Richard Schwartz (17:02) Both methods actually are one and the same. Doerfler's Borchardt's algorithm + Richardson extrapolation: d(k,k) ≈ ∫(n^t,t,0,1) = (n1)/ln(n) CAS> d00 := a0 := (1+n)/2 CAS> g0 := sqrt(n) CAS> d01 := factor(a1 := (a0+g0)/2) → (n + 2*n^(1/2) + 1)/4 CAS> d11 := factor(d01 + (d01d00)/3) → (n + 4*n^(1/2) + 1)/6 CAS> g1 := sqrt(a1*g0) → (sqrt(sqrt(n)*(n+1)+2*n))/2 This seemingly looking "AGM" is really to evaluate doubled points. (*) We can rewrite g1 = (√√n) * (1+√n)/2, expanded: CAS> g1 := (n^(1/4) + n^(3/4))/2 CAS> d02 := factor(a2 := (a1+g1)/2) → (n + 2*n^(3/4) + 2*n^(1/2) + 2*n^(1/4) + 1)/8 CAS> d12 := factor(d02 + (d02d01)/3) → (n + 4*n^(3/4) + 2*n^(1/2) + 4*n^(1/4) + 1)/12 CAS> d22 := factor(d12 + (d12d11)/15) → (7*n + 32*n^(3/4) + 12*n^(1/2) + 32*n^(1/4) + 7)/90 Or, we dot multiply [a0,a1,a2] by weight, directly for Boole's Rule CAS> simplify( [a0, a1, a2](n=q^4) * [1,20,64]/45 ) \( \Large \frac {(7\cdot q^{4}+32\cdot q^{3}+12\cdot q^{2}+32\cdot q+7)} {90} \) (*) We can "see" doubling of points by thinking weights as a "number" weight of a3 = 122222221 = 11111111*11 weight of g2 = 1010101 = 11111111/11 weight of g3 = √(a3*g2) = 11111111 

10222021, 02:15 PM
(This post was last modified: 03102022 06:14 PM by Albert Chan.)
Post: #10




RE: Estimate logarithm quickly
Instead of Romberg integration, we can also do Gaussian Quadrature, to estimate ln(n)
Shifting points and weights to limit of 0 to 1, for 3 points, we have: Code: point weight ln(n) ≈ (n1) * 18/(((n^p + n^p)*5 + 8)*√n) , where p = √(0.15) Above formula give ln(e) ≈ 1.0000004796, slightly better than Booles' rule. However, it maybe expensive to calculate n^p, n^p The good news is (n^p + n^p) taylor series, half the terms get cancelled. Only even powers of p remains, no square roots anywhere. CAS> s := (1+x)^p + (1+x)^p CAS> series(s) 2 + p^2*x^2  p^2*x^3 + (1/12*p^4+11/12*p^2)*x^4 + (1/6*p^45/6*p^2)*x^5 + x^6*order_size(x) CAS> pade(s, x, 4, 3) (p=sqrt(0.15)) (3.5*x^224*x24)/(0.85*x^212*x12) Replacing s with above pade approximation, we have: ln(n) ≈ g  g/(2.7 + 24/g^2) , where g = (n1)/√n (*) With this, we have log(e) ≈ 1.00016029872 Not as good, but still twice as good as Doerfler's formula (log(e) ≈ 1/1.00033715273) If we limit n = [√0.5, √2], this Pade formula is at least 17 times better. lua> n = 1.2586 lua> (n1) * 6/(1 + 4*sqrt(n) + n)  Doerfler formula 0.22999976897818544 lua> g = (n1)/sqrt(n) lua> g  g / (2.7 + 24/g^2)  Pade formula 0.22999999788822215 lua> log(n) 0.2299999921106961 (*) If we consider whole expression (not just 3 points), we still get the same result. CAS> pade(ln(1+x)*sqrt(1+x)/x,x,4,3) → (17*x^2+240*x+240)/(27*x^2+240*x+240) → ln(1+x) ≈ x/sqrt(1+x) * (1  1/(27/10 + 24*(1+x)/x^2)) ≈ x  x^2/2 + x^3/3  x^4/4 + x^5/5  x^6/6 + x^7/6.992  x^8/7.962 + x^9/8.897  x^10/9.786 + ... 

11202021, 02:02 PM
Post: #11




RE: Estimate logarithm quickly
(10062021 10:58 PM)Albert Chan Wrote: Thanks for the video. I finally "get" Doerfler's formula ∫(n^t, t, 0, 1) = ∫(e^(t*ln(n)), t, 0, 1) = ∫(e^x, x, 0, ln(n)) / ln(n) = (n1) / ln(n) With exponential function, its derivative, (e^x)' = e^x, also grow exponentially. Any integrand polynomial fit that included endpoints will overestimate integral result. (Unless n=1, e^(t*ln(n)) = e^0 = 1, no longer exponential) Because integral is in the denominator, ln(n) is underestimated. (we compare absolute value, to cover cases when 0 < n < 1) ln(n) ≥ n1 / ((1 + 4*√n + n)/6) // Simpson's Rule Error = 0 when n=1, and increases when n is further away from 1. Example, ln(√3) will underestimate more than ln(√2) ln(3)/ln(2) = ln(√3)/ln(√2) > (6*(√31)/(4*√(√3) + √3+1)) / (6*(√21)/(4*√(√2) + √2+1)) ≈ 1.58492 ln(3)/ln(2) = 1.58496..., inequality holds, as expected 

11222021, 01:49 AM
Post: #12




RE: Estimate logarithm quickly
Here is an example, to prove 13^40 > 4^74, difference of only 1.2%
The video make use of (13^3 = 2197) ≈ (2^11 = 2048). Here, we don't use it. We compare ratio of 4th root, (difference is now tiny 0.3%, but it doesn't matter) FYI, 37/10 is a convergent of log2(13), thus below ratio is very close to 1. 13^10 / 2^37 = (13/8)^3 / (16/13)^7 RHS make the base close to 1, which is needed for Doerfler's formula to work well. Note: Constant 6 of Doerfler's formula is skipped. (it will be cancelled anyway) Comparing numerator/denominator ratio with its logarithm, an increasing function. (3/2*ln(169/64)) / (7/2*ln(256/169)) > (3/7) * ((169/641)/(4*(13/8) + 169/64+1)) / ((256/1691)/(4*(16/13) + 256/169+1)) = (3/7) * (105/649) / (29/419) = 18855/18821 ≈ 1.00181 > 1 // ⇒ 13^10 > 2^37, QED We give up a bit of accuracy by doing ln(x^2) instead of ln(x), to avoid square root mess. Had we use ln(x) instead, we would get a ratio of 1.00208 (true ratio ≈ 1.00210)  Actually, we can do better, (13/8 = 1.625) ≈ (256/169 ≈ 1.515) Both underestimated errors have similar size, thus mostly cancel each other. Apply Doerfler's formula for (3*ln(13/8)) / (7/2*ln(256/169)), we get a ratio of 1.00209 

11232021, 09:01 PM
(This post was last modified: 04262022 11:09 AM by Albert Chan.)
Post: #13




RE: Estimate logarithm quickly
(11222021 01:49 AM)Albert Chan Wrote: FYI, 37/10 is a convergent of log2(13), thus below ratio is very close to 1. Another way to proof above ratio > 1 is to proof log2(13) > 37/10 However, we have to get the right inequality (here, we wanted ">", not "<") We can proof this instead, which leads to the same thing. To proof: ln(13/8)/ln(16/13) > (3.73)/(43.7) = 7/3 ≈ 2.33333 We use a simple formula to test, ln(x) = 2*atanh(y), where y = (x1)/(x+1) Y(x) := (x1)/(x+1); // if x is fractions, Y(n/d) = (nd)/(n+d) y1>y2>0: atanh(y1)/atanh(y2) > (y1/y2) ln(13/8)/ln(16/13) > Y(13/8)/Y(16/13) = (5/21)/(3/29) = 145/63 ≈ 2.30159 atanh(y) ≈ y not enough. We add pade correction: atanh(y)/y > 1/ (1  y^2/3) y1>y2>0: atanh(y1)/atanh(y2) > (y1/y2) * (3y2^2)/(3y1^2) (145/63) * (3(3/29)^2)/(3(5/21)^2) = 43995/18821 ≈ 2.33755 > 7/3 QED  8/13 < 1 < 16/13 → log2(13) = (3/1, 4/1) Try next semiconvergent: (3+4)/(1+1) = 7/2 (*) 2^7/13^2 = 128/169 < 1 → log2(13) = (7/2, 4/1) (169/128 ≈ 1.32031) > (16/13 ≈1.23077) We can keep the same denominator, to ensure we are solving for tighter z lower bound. ln(169/128)/ln(16/13) > Y(169/128)/Y(16/13) = (41/297)/(3/29) = 1189/891 solve((2z7)/(4z) > 1189/891, z), we get z = [3.70010, 4) → log2(13) > 37/10 QED (*) We do not want log2(13) next convergent, (3+4*2)/(1+1*2) = 11/3, because it guaranteed to be in the denominator, solving the tighter upper bound instead. 

11262021, 06:24 PM
Post: #14




RE: Estimate logarithm quickly
(11202021 02:02 PM)Albert Chan Wrote: ln(n) ≥ n1 / ((1 + 4*√n + n)/6) // Simpson's Rule A simpler proof of inequality is convert it to atanh(y) Assume x>1, then y = (x1)/(x+1) > 0 To avoid square mess, we apply Doerfler's formula with squared argument. ln(x) = atanh(y = (x1)/(x+1))*2 → atanh(y) = ln(x = (1+y)/(1y))/2 XCAS> D2(x) := 3*(x*x1)/(1 + 4*x + x*x) XCas> factor(D2((1+y)/(1y)) /2) → 3*y/(3+y^2) This is just pade(atanh(y),y,4,2), which expands to: y/(1y^2/3) = y + y^3/3 + y^5/3² + y^7/3³ + ... For y>0, atanh(y) = y + y^3/3 + y^5/5 + y^7/7 + ... is bigger. For x>1, ln(x), which atanh(y) were derived from, is biggger than D2(x) Because of symmetry, For 0<x<1, D2(x) = D2(1/x), same as ln(x). Thus the proof can be extended from x > 1, to x > 0 XCAS> D2(2), D2(1/2) → (9/13 , 9/13) 

04262022, 05:01 PM
Post: #15




RE: Estimate logarithm quickly
(11232021 09:01 PM)Albert Chan Wrote: 8/13 < 1 < 16/13 → log2(13) = (3/1, 4/1) To proof log2(13) > 3.7 with next convergent is equivalent to: R = ln(13^3/2^11) / ln(16/13) > (3*3.711) / (43.7) = 0.1/0.3 R = atanh(t1 = 149/4245) / atanh(t2 = 3/29) > 1/3 Since t1 < t2, if we estimate atanh(x) ≈ x, we get R upper bound instead, not what we wanted. R < (t1/t2) ≈ 0.339301 Tighter atanh bound is better, but still the wrong bound. R < (t1/t2) * (1t2²/3)/(1t1²/3) ≈ 0.338230 (t1/t2) matched 2 digits. We may claim better R good for 0.338 > 1/3, but there is surer way. A trick is to convert atanh to asinh, with alternate sign series expansion see Arc SOHCAHTOA method, for hyperbolics R = atanhq(149²/4245²) / atanhq(3²/29²) // TOA = asinhq(s1 = 149²/(4245²149²)) / asinhq(s2 = 3²/(29²3²)) // SOH R = sqrt(s1/s2) * [1, (1+s2/6)/(1+s1/6)] = [0.337689, 0.338228] > 1/3 QED Instead of atanh 1sided estimate, asinh allowed ratio estimate be bracketed. (here, we don't need upper bound. R ≥ sqrt(s1/s2) > 1/3 is enough) Bonus: asinh converge faster than atanh. Below R bounds, last factor (1+ε) → 1/(1ε), for better pade estimate. R < (t1/t2) * (1+(t1²/3)/(13/5*t1²)) / (1+(t2²/3)/(13/5*t2²)) = 0.338226 271757 // error = 9602 ulp R > sqrt(s1/s2) * (1(s1/6)/(1+9/20*s1)) / (1(s2/6)/(1+9/20*s2)) = 0.338226 257560 // error = +4595 ulp 

04262022, 05:52 PM
Post: #16




RE: Estimate logarithm quickly
Sorry to hijack this topic, I found this estimation in czech written book (published in 1984)
log (x) = (z^2*a+b)*z z = (x1)/(x+1) and limits for x are: for decadic logarithm the x must be between 0,316 and 3,162 for natural logharithm the x must be between 0,6065 and 1,649 the constants a and b: for log(x): a = 0,36415, b = 0,86304 for ln(x): a= 0,70225, b = 1,99938 

04272022, 10:02 PM
Post: #17




RE: Estimate logarithm quickly
(10222021 02:15 PM)Albert Chan Wrote: ln(n) ≈ g  g/(2.7 + 24/g^2) , where g = (n1)/√n Another way to derive above ln(n) approximation, via asinh see Arc SOHCAHTOA method, for hyperbolics log(n) = 2*atanh((n1) / (n+1)) // assumed n ≥ 1, for atanhq = 2*atanhq((n1)² / (n+1)²) // TOA, O=(n1)², A = (n+1)² = 2*asinhq((n1)² / (4n)) // SOH, H = AO = 4n = 2*asinh(g/2) // remove n ≥ 1 assumption. = g * (1  g²/24 * (1  9/80*g² * (1  25/168*g² * ( ... ≈ g * (1  g²/24 / (1 + 9/80*g²)) = g  g/(2.7 + 24/g²) Or, we could simply use pade command. CAS> p := pade(2*asinh(g/2),g,5,4) → (17*g^3+240*g)/(27*g^2+240) CAS> partfrac(g/(gp)) → 27/10+24/g^2 With 3 asinh taylor terms to estimate log(n), formula always overestimate, in absolute value. lua> function fastlog(x) x=(x1)/sqrt(x); return x  x/(2.7+24/(x*x)) end lua> function ratio(x) return fastlog(x) / log(x) end lua> ratio(sqrt(2.0)) 1.0000002929130414 lua> ratio(sqrt(0.5)) 1.0000002929130412 lua> x = 0.1 lua> fastlog((1+x)/(1x))/2  atanh(x) upperbound estimate 0.10033534884362541 lua> atanh(x) 0.10033534773107558 

04282022, 04:33 PM
Post: #18




RE: Estimate logarithm quickly  
04282022, 05:02 PM
Post: #19




RE: Estimate logarithm quickly
Hry s kalkulátory (Playing with Calculators)


« Next Oldest  Next Newest »

User(s) browsing this thread: