HP Forums
Estimate logarithm quickly - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Estimate logarithm quickly (/thread-17376.html)



Estimate logarithm quickly - Albert Chan - 08-21-2021 03:39 PM

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


RE: Estimate logarithm quickly - Albert Chan - 08-21-2021 03:58 PM

(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)


RE: Estimate logarithm quickly - Albert Chan - 08-21-2021 05:41 PM

(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


RE: Estimate logarithm quickly - trojdor - 08-21-2021 06:31 PM

Albert, thanks for the book reference.
(Surprisingly, that's one I don't have...yet.)


RE: Estimate logarithm quickly - Albert Chan - 08-21-2021 11:42 PM

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


RE: Estimate logarithm quickly - EdS2 - 08-23-2021 06:44 AM

(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/~sgtatham/spigot/spigot.html#cfrac


RE: Estimate logarithm quickly - Albert Chan - 10-06-2021 10:58 PM

(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


RE: Estimate logarithm quickly - Csaba Tizedes - 10-07-2021 06:25 PM

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


RE: Estimate logarithm quickly - Albert Chan - 10-18-2021 01:02 PM

(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


RE: Estimate logarithm quickly - Albert Chan - 10-22-2021 02:15 PM

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))


RE: Estimate logarithm quickly - Albert Chan - 11-20-2021 02:02 PM

(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


RE: Estimate logarithm quickly - Albert Chan - 11-22-2021 01:49 AM

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


RE: Estimate logarithm quickly - Albert Chan - 11-23-2021 09:01 PM

(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.


RE: Estimate logarithm quickly - Albert Chan - 11-26-2021 06:24 PM

(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)