Post Reply 
Estimate logarithm quickly
08-21-2021, 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) * (n-x) * (n*n/(n*n-x*x))

2*ln(n) = ln(n+x) + ln(n-x) + ln(n*n/(n*n-x*x))

ln(n*n/(n*n-x*x)) = ln(1-x*x/(n*n-x*x)) = 2*atanh(x*x/(2*n*n-x*x))

Let a = ln(n+x), b = ln(n-x), c = n/x, and halved both side.

ln(n) = (a+b)/2 + atanh(1/(2*c*c-1))

atanh(1/(2*c*c-1)) > atanh(1/c) / (2c) = (a-b)/2 / (2c)

ln(n) > (a+b)/2 + (a-b)/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(3-1) = 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(3-1.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
Find all posts by this user
Quote this message in a reply
08-21-2021, 03:58 PM
Post: #2
RE: Estimate logarithm quickly
(08-21-2021 03:39 PM)Albert Chan Wrote:  atanh(1/(2*c*c-1)) > atanh(1/c) / (2c) = (a-b)/2 / (2c)

This is not yet proofed (n-x > 0, c = n/x > 1)
To proof it, it required this inequality (also not proofed, suggestions welcome):

(1-w)^(2-w) * (1+w)^(2+w) < 1,       where 0 < w < 1

XCAS> series((1-w)^(2-w) * (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 ...

(1-w)^(2-w) * (1+w)^(2+w) < 1
(1-w*w)^2 < ((1-w)/(1+w))^w
(1/(1-w*w))^2 > ((1+w)/(1-w))^w

Let c = 1/w, c > 1

(c*c/(c*c-1))^2 > ((1+1/c)/(1-1/c)) ^ (1/c)

ln(c*c/(c*c-1)) * 2 > ln((1+1/c)/(1-1/c)) / c

Simpilfy with identity ln(z) = 2*atanh((z-1)/(z+1))

(2*atanh(1/(2*c*c-1)) * 2 > (2*atanh(1/c)) / c

atanh(1/(2*c*c-1) > atanh(1/c) / (2c)
Find all posts by this user
Quote this message in a reply
08-21-2021, 05:41 PM
Post: #3
RE: Estimate logarithm quickly
(08-21-2021 03:58 PM)Albert Chan Wrote:  To proof it, it required this inequality (also not proofed, suggestions welcome):

(1-w)^(2-w) * (1+w)^(2+w) < 1,       where 0 < w < 1

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 := (2-w)*ln(1-w) + (2+w)*ln(1+w)
XCAS> diff(f,w)

-ln(1-w) - (2-w)/(1-w) + 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)/(1-x))/2

(w^2-1) * atanh(w) + w = (w^2-1) * (w + w^3/3 + w^5/5 + ...) + w

All coefficients of w's is positive, 1/(2k-1) - 1/(2k+1) = 2/(4k^2-1) > 0

XCAS> series((w^2-1)*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((1-w)^(2-w) * (1+w)^(2+w)) < 0       → (1-w)^(2-w) * (1+w)^(2+w) < 1
Find all posts by this user
Quote this message in a reply
08-21-2021, 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 > =
Find all posts by this user
Quote this message in a reply
08-21-2021, 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*(n-1) / (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
Find all posts by this user
Quote this message in a reply
08-23-2021, 06:44 AM
Post: #6
RE: Estimate logarithm quickly
(08-21-2021 03:39 PM)Albert Chan Wrote:  Estimate logorithm quickly, with tiny errors.
...
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
...

Just to note, this is Simon Tatham's arbitrary-precision spigot calculator, which deals in continued fractions as well as decimals:
https://www.chiark.greenend.org.uk/~sgta...html#cfrac
Find all posts by this user
Quote this message in a reply
10-06-2021, 10:58 PM
Post: #7
RE: Estimate logarithm quickly
(10-06-2021 05:12 PM)Eric Rechlin Wrote:  Computing Logarithms by Integration - Richard Schwartz (17:02)

Thanks for the video. I finally "get" Doerfler's formula Smile

ln(n) = (n-1) / ∫(n^t, t, 0, 1) ≈ (n-1) / ((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 + (S2-S1)/(4^2-1)     → [7/90, 16/45, 2/15, 16/45, 7/90]
CAS> horner(B1, e^0.25) / (e-1)       → 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*(x-1)) // 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*(e-1))   → 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 + (s2-s1)/(4^2-1)              → 1.00000050019
Find all posts by this user
Quote this message in a reply
10-07-2021, 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
Find all posts by this user
Quote this message in a reply
10-18-2021, 01:02 PM (This post was last modified: 10-18-2021 04:30 PM by Albert Chan.)
Post: #9
RE: Estimate logarithm quickly
(10-06-2021 10:58 PM)Albert Chan Wrote:  
(10-06-2021 05:12 PM)Eric Rechlin Wrote:  Computing Logarithms by Integration - Richard Schwartz (17:02)

Thanks for the video. I finally "get" Doerfler's formula Smile

Both methods actually are one and the same.
Doerfler's Borchardt's algorithm + Richardson extrapolation: d(k,k) ≈ ∫(n^t,t,0,1) = (n-1)/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 + (d01-d00)/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 + (d02-d01)/3)    → (n + 4*n^(3/4) + 2*n^(1/2) + 4*n^(1/4) + 1)/12
CAS> d22 := factor(d12 + (d12-d11)/15)   → (7*n + 32*n^(3/4) + 12*n^(1/2) + 32*n^(1/4) + 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
Find all posts by this user
Quote this message in a reply
10-22-2021, 02:15 PM (This post was last modified: 10-22-2021 05:40 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
0.5            8/18
0.5 ± √(0.15)  5/18

ln(n) ≈ (n-1) * 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^4-5/6*p^2)*x^5 + x^6*order_size(x)

CAS> pade(s, x, 4, 3) (p=sqrt(0.15))

(-3.5*x^2-24*x-24)/(-0.85*x^2-12*x-12)

Replacing s with above pade approximation, we have:

ln(n) ≈ g - g/(2.7 + 24/g^2)       , where g = (n-1)/√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> (n-1) * 6/(1 + 4*sqrt(n) + n) -- Doerfler formula
0.22999976897818544
lua> g = (n-1)/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))
Find all posts by this user
Quote this message in a reply
11-20-2021, 02:02 PM
Post: #11
RE: Estimate logarithm quickly
(10-06-2021 10:58 PM)Albert Chan Wrote:  Thanks for the video. I finally "get" Doerfler's formula Smile

ln(n) = (n-1) / ∫(n^t, t, 0, 1) ≈ (n-1) / ((1 + 4*√n + n)/6) // Simpson's Rule.

∫(n^t, t, 0, 1) = ∫(e^(t*ln(n)), t, 0, 1) = ∫(e^x, x, 0, ln(n)) / ln(n) = (n-1) / ln(n)

With exponential function, its derivative, (e^x)' = e^x, also grow exponentially.
Any integrand polynomial fit that included end-points will over-estimate integral result.
(Unless n=1, e^(t*ln(n)) = e^0 = 1, no longer exponential)

Because integral is in the denominator, |ln(n)| is under-estimated.
(we compare absolute value, to cover cases when 0 < n < 1)

|ln(n)| |n-1| / ((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*(√3-1)/(4*√(√3) + √3+1)) / (6*(√2-1)/(4*√(√2) + √2+1))
≈ 1.58492      

ln(3)/ln(2) = 1.58496..., inequality holds, as expected
Find all posts by this user
Quote this message in a reply
11-22-2021, 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/64-1)/(4*(13/8) + 169/64+1)) / ((256/169-1)/(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 under-estimated 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
Find all posts by this user
Quote this message in a reply
11-23-2021, 09:01 PM (This post was last modified: 11-29-2021 10:13 PM by Albert Chan.)
Post: #13
RE: Estimate logarithm quickly
(11-22-2021 01:49 AM)Albert Chan Wrote:  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

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 use a simple formula to test, ln(x) = 2*atanh(y), where y = (x-1)/(x+1)

Y(x) := (x-1)/(x+1);       // if x is fractions, Y(n/d) = (n-d)/(n+d)

For y>0, atanh'(y) = 1/(1-y^2) is increasing function ⇒ 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

Let z = log2(3),
solve((z-3)/(4-z) > 145/63, z), we get z = [3.69712, 4), not good enough.

Note that solved z is always bracketed by a range. (if LHS go negative, we have no solution)
With above setup, denominator must be one side of the bound.
This is how we can be sure we are solving for right bound for z.

We could fix by a pade correction, or start over with another ratio of logs.
|atanh(y)| ≥ |y / (1 - y^2/3)|

y1>y2>0: atanh(y1)/atanh(y2) > y1/y2 * (3-y2^2)/(3-y1^2)

(5/21)/(3/29) = 145/63 ≈ 2.30159
145/63 * (3-(3/29)^2)/(3-(5/21)^2) = 43995/18821 ≈ 2.33755

ln(13/8)/ln(16/13) ≈ 2.33823, correction not bad.

solve((z-3)/(4-z) > 43995/18821, z), we get z = [3.70038, 4)    → log2(13) > 37/10 QED

---

8/13 < 1 < 16/13                  → log2(13) = (3/1, 4/1)

Try next semi-convergent: (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((2z-7)/(4-z) > 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.
Find all posts by this user
Quote this message in a reply
11-26-2021, 06:24 PM
Post: #14
RE: Estimate logarithm quickly
(11-20-2021 02:02 PM)Albert Chan Wrote:  |ln(n)| |n-1| / ((1 + 4*√n + n)/6) // Simpson's Rule

A simpler proof of inequality is convert it to atanh(y)
Assume x>1, then y = (x-1)/(x+1) > 0

To avoid square mess, we apply Doerfler's formula with squared argument.
ln(x) = atanh(y = (x-1)/(x+1))*2       → atanh(y) = ln(x = (1+y)/(1-y))/2

XCAS> D2(x) := 3*(x*x-1)/(1 + 4*x + x*x)
XCas> factor(D2((1+y)/(1-y)) /2)       → -3*y/(-3+y^2)

This is just pade(atanh(y),y,4,2), which expands to:

y/(1-y^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)
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)