Estimate logarithm quickly

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

(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

[trojdor]

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

[Albert Chan]

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

[EdS2]

(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

[Albert Chan]

(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

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

[Csaba Tizedes]

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

[Albert Chan]

(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

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

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

[Albert Chan]

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))
≈ 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 + ...

[Albert Chan]

(10-06-2021 10:58 PM)Albert Chan Wrote:  Thanks for the video. I finally "get" Doerfler's formula

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

[Albert Chan]

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

[Albert Chan]

(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 can proof this instead, which leads to the same thing.
To proof: ln(13/8)/ln(16/13) > (3.7-3)/(4-3.7) = 7/3 ≈ 2.33333

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)

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) * (3-y2^2)/(3-y1^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 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.

[Albert Chan]

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

---

We can use above ln(n) estimate for atanh(x) = 1/2 * ln((1+x)/(1-x))

atanh(x)/x ≥ 3/(1 + 2*sqrt(1-x*x))

Cas> series((atanh(x)/x) / (3/(1+2*sqrt(1-x*x))), x, 0, 9, polynom)

1 + 1/180*x^4 + 17/2520*x^6 + 139/20160*x^8

Since i^4=1, inequality hold even if x is purely imaginery! (x → x*i, and simplify)

atan(x)/x ≥ 3/(1 + 2*sqrt(1+x*x))

lua> y, x = 0.75, 0.5

lua> atanh(x)/x, 3/(1+2*sqrt(1-x*x))
1.0986122886681096      1.098076211353316
lua> atanh(y)/y, 3/(1+2*sqrt(1-y*y))
1.2972734327035422      1.2915026221291812

lua> atan(x)/x, 3/(1+2*sqrt(1+x*x))
0.9272952180016122      0.9270509831248422
lua> atan(y)/y, 3/(1+2*sqrt(1+y*y))
0.8580014783910458      0.8571428571428571

Convergence is from 1 side, with estimate not as good for bigger y.
We can compare ratio for y ≥ x > 0, using RHS numbers.

lua> (atanh(y)/y) / (atanh(x)/x) >= (1+2*sqrt(1-x*x)) / (1+2*sqrt(1-y*y))
true
lua> (atan(y)/y) / (atan(x)/x) >= (1+2*sqrt(1+x*x)) / (1+2*sqrt(1+y*y))
true

[Albert Chan]

(11-23-2021 09:01 PM)Albert Chan Wrote:  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)
...

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

To proof log2(13) > 3.7 with next convergent is equivalent to:

R = ln(13^3/2^11) / ln(16/13) > (3*3.7-11) / (4-3.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) * (1-t2²/3)/(1-t1²/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 1-sided 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)/(1-3/5*t1²)) / (1+(t2²/3)/(1-3/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

[klesl]

Sorry to hijack this topic, I found this estimation in czech written book (published in 1984)
log (x) = (z^2*a+b)*z
z = (x-1)/(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

[Albert Chan]

(10-22-2021 02:15 PM)Albert Chan Wrote:  ln(n) ≈ g - g/(2.7 + 24/g^2)       , where g = (n-1)/√n

Another way to derive above ln(n) approximation, via asinh
see Arc SOHCAHTOA method, for hyperbolics

log(n)
= 2*atanh((n-1) / (n+1))                // assumed n ≥ 1, for atanhq
= 2*atanhq((n-1)² / (n+1)²)           // TOA, O=(n-1)², A = (n+1)²
= 2*asinhq((n-1)² / (4n))               // SOH, H = A-O = 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/(g-p))                          → 27/10+24/g^2

With 3 asinh taylor terms to estimate log(n), formula always over-estimate, in absolute value.

lua> function fastlog(x) x=(x-1)/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)/(1-x))/2 -- atanh(x) upper-bound estimate
0.10033534884362541
lua> atanh(x)
0.10033534773107558

[vaklaff]

(04-26-2022 05:52 PM)klesl Wrote:  Sorry to hijack this topic, I found this estimation in czech written book (published in 1984)

Side question - what’s the name of the book?
Thanks in advance

[klesl]

Hry s kalkulátory (Playing with Calculators)