Ln(x) using repeated square root extraction
03-21-2016, 12:03 AM (This post was last modified: 04-11-2020 12:00 AM by Gerson W. Barbosa.)
Post: #1
 Gerson W. Barbosa Senior Member Posts: 1,465 Joined: Dec 2013
Ln(x) using repeated square root extraction
Algorithm:

1. Enter x;
2. Take the square root 6 times;
3. Multiply by 2;
4. Subtract 1;
5. Take the square root;
6. Subtract 1;
7. Multiply by 64.

For x = 2, on an 8-digit calculator like the Canon LC-37, you should get ln(2) ~ 0.693152. For 10-digit calculators, in step 2 take the square root 9 times and multiply by 512 in the last step. Thus, on the HP-12C, which has two guard digits, we get ln(2) ~ 0.693147648. Ordinary 10-digit calculator might give less approximate results.

A 4-step less accurate algorithm is available in this old thread:

The accuracy range should be tested for your particular calculator. Use these at your own risk.

Gerson.

Method:

$\ln (x)=\int_{1}^{x}\frac{du}{u}=\int_{1}^{x}u^{-1}du=\lim_{v\rightarrow 0}\left | \frac{u^{v}}{v} \right |_{1}^{x}=\lim_{v\rightarrow 0} \left ( \frac{x^{v}}{v}-\frac{1^{v}}{v} \right )=\lim_{v\rightarrow 0} \left ( \frac{x^{v}-1}{v} \right )$

By testing the limit with a few values for v close to zero, an empirical second term has been able to be added to the expression:

$\ln (x)=\lim_{v\rightarrow 0} \left ( \frac{x^{v}-1}{v}-\frac{v\cdot \ln^{2}(x)}{2} \right )$

After turning the limit into an equality and solving the resulting quadratic equation for ln(x), we get

$\ln (x)=\lim_{v\rightarrow 0} \left ( \frac{\sqrt{2\cdot x^{v}-1}-1}{v} \right )$

Example:

Let x = 2 and v = 0.001

Then

$\ln (2)\approx \frac{\sqrt{2\cdot 2^{0.001}-1}-1}{0.001}\approx 0.69314723$

which is good to 6 decimal places ( ln(2) = 0.69314718 ).

P.S.: Perhaps 8 square root extractions and final multiplication by 256 (2^8) is a better compromise for the range 2..100 on the Canon LC-37. Contrary to what I thought its SILVA-CELL 189 battery is still working after all these years, so I was able to make some more tests. :-)

-----------------------------------

Good Friday 2020 Update

The following expression allows us to go further:

$\ln (x)=\lim_{v\rightarrow 0} \left ( \frac{x^{v}-1}{v}-\frac{v\cdot \ln^{2}(x)}{2!} -\frac{v^{2}\cdot \ln^{3}(x)}{3!} -\frac{v^{3}\cdot \ln^{4}(x)}{4!} - \cdots \right )$

By going up to the third power of ln(x) and solving the corrisponding cubic equation, we get

$\ln (x)=\lim_{v\rightarrow 0}\frac{1}{v}\left ( \sqrt[3]{\sqrt{9x^{2v}-6x^{v}+2}+3x^{v}-1} - \frac{1}{\sqrt[3]{\sqrt{9x^{2v}-6x^{v}+2}+3x^{v}-1}} - 1 \right )$

For example, if x = 2 and v = 0.01 then ln(2) ~ 0.69314719

This approximation is rather complicated for manual calculations. Also, it requires one cubic root extraction. It can be used, however, as an alternative method for implementing the ln(x) function to 7 significant digits on calculators that lack it, like the HP-16C. The e^x function has also been implemented. Since it's based on the first formula, obtained by solving a quadratic equation, and uses a constant number of loops, it gives less significant digits. The cubic root algorithm is optimized for the narrow range required by the approximation formula and was based one provided by Albert Chan here.

Code:
 001- LBL A; LN 002- 1 003- EEX 004- CHS 005- 2 006- CF 1 007- 1 008- + 009- STO I 010- CLx 011- LASTx 012- x⇆y 013- x≤y 014- GTO 2 015- LBL 0 016- √x 017- RCL I 018- x>y 019- GSB 1 020- R↓ 021- x⇆y 022- ENTER  023- + 024- x⇆y 025- GTO 0 026- RTN 027- LBL 1 028- R↓ 029- STO 0 030- R↓ 031- 2 032- × 033- STO 1 034- RCL 0 035- ENTER  036- × 037- 9 038- × 039- 6 040- RCL 0 041- × 042- - 043- 2 044- + 045- √x 046- 3 047- RCL 0 048- × 049- + 050- 1 051- - 052- GSB 3 053- ENTER  054- 1/x 055- - 056- 1 057- - 058- RCL 1 059- × 060- F?1 061- CHS 062- RTN 063- LBL 2 064- SF1 065- 1/x 066- RTN  067- LBL 3 068- 2 069- STO I 070- ENTER  071- + 072- x⇆y 073- × 074- STO 0 075- LASTx 076- 2 077- - 078- 6 079- / 080- √x 081- 1 082- + 083- LBL 4 084- ENTER  085- ENTER  086- RCL 0 087- x⇆y 088- / 089- LASTx 090- ENTER  091- × 092- - 093- 3 094- / 095- √x 096- + 097- 2 098- / 099- DSZ 100- GTO 4 101- RTN 102- LBL B; eᵡ 103- 1 104- 3 105- STO I 106- R↓ 107- 8 108- 1 109- 9 110- 2 111- + 112- LASTx 113- / 114- ENTER 115- × 116- 1 117- + 118- 2 119- / 120- LBL 5 121- ENTER  122- × 123- DSZ 124- GTO 5 125- RTN

Examples:

2 GSB A -> 0.6931471(360)
GSB B -> 2.0000(13165)

12345 GSB A -> 9.421006(848)
GSB B -> 12345.0(1957)

6.789 EEX 79 GSB A -> 183.8195(343) GSB B -> 6.(686887E79)

0.12345 GSB A -> -2.091919(104)
GSB B -> 0.123450(119)

230 GSB B -> 7.496895E99
GSB A -> 229.97041(15)

1 GSB A -> 0.000000000
GSB B -> 1.000000000
GSB B -> 2.7182(92170)
GSB B -> 15.1544(4181)
GSB A -> 2.71829(4016)
GSB A -> 1.000004(736)
GSB A -> 0.000004736

0.693147181 CHS GSB B -> 0.500000(16)

10 GSB A -> 2.302585(344) STO 2

2 GSB A RCL 2 / -> 0.3010299(435)

0 GSB A -> Error 0

2 CHS GSB A -> Error 0
03-21-2016, 05:09 AM (This post was last modified: 03-21-2016 05:12 AM by Gerson W. Barbosa.)
Post: #2
 Gerson W. Barbosa Senior Member Posts: 1,465 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
For the sake of completeness here is the recursive limit definition of ln(x), even though this is not quite necessary for our purpose:

$\ln (x)=\lim_{n\rightarrow \infty} \left [ n\cdot \left ( x^{\frac{1}{n}}-1 \right )-\sum_{k=2}^{\infty }\frac{\ln ^{k}(x)}{k!\cdot n^{k-1}} \right ]$

If the limit is removed and n is set to 1 then we'll have the following recursive series representation:

$\ln (x)= \ x-1 -\sum_{k=2}^{\infty }\frac{\ln ^{k}(x)}{k!}$
03-21-2016, 07:11 AM
Post: #3
 Paul Dale Senior Member Posts: 1,752 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
I've got a feeling I tried something along these lines for the LN function of the 34S. The LN function has always been very slow on the 34S and I've tried several implementations to try to get both performance and accuracy. The current implementation is much faster than the original but still far slower than I'd like. I ended up with:

Code:
/* Natural logarithm.  *  * Take advantage of the fact that we store our numbers in the form: m * 10^e  * so log(m * 10^e) = log(m) + e * log(10)  * do this so that m is always in the range 0.1 <= m < 2.  However if the number  * is already in the range 0.5 .. 1.5, this step is skipped.  *  * Then use the fact that ln(x^2) = 2 * ln(x) to range reduce the mantissa  * into 1/sqrt(2) <= m < 2.  *  * Finally, apply the series expansion:  *   ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)  * which converges quickly for an argument near unity.  */

- Pauli
03-03-2022, 11:20 PM
Post: #4
 Thomas Klemm Senior Member Posts: 1,687 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
Formula

We can use the following identity:

$$\log(x) = n \log\left(x^{\frac{1}{n}}\right) = n \log\left(\sqrt[n]{x}\right)$$

For say $$n = 2^{10} = 1024$$ and $$1 \leqslant x \leqslant 100$$ the value of $$\sqrt[n]{x}$$ is close to $$1$$.

Thus we can use the Taylor series to calculate the logarithm:

$$\log(1 + \varepsilon) = \varepsilon - \frac{\varepsilon^2}{2} + \frac{\varepsilon^3}{3} - \frac{\varepsilon^4}{4} + \frac{\varepsilon^5}{5} + \mathcal{O}(\varepsilon^6)$$

Program

Here's a program for the HP-42S that calculates both the logarithm and its approximation:
Code:
LN  LASTX SQRT  SQRT  SQRT  SQRT  SQRT SQRT  SQRT  SQRT  SQRT  SQRT 1  - 4  1/X  RCL× ST Y 3  1/X  X<>Y  -  RCL× ST Y 2  1/X  X<>Y  -  RCL× ST Y 1  X<>Y  -  × 1024  ×

It's easy to extend if you want to use more terms.

Example

x = 2

0.69314718056
0.69314718056

Comparison

We can compare this to your solution:

$$\sqrt{2(1 + \varepsilon) - 1} - 1$$

The Taylor series agrees for the first two terms:

$$\varepsilon - \frac{\varepsilon^2}{2} + \frac{\varepsilon^3}{2} - \frac{5 \varepsilon^4}{8} + \frac{7 \varepsilon^5}{8} + \mathcal{O}(\varepsilon^6)$$
03-05-2022, 05:32 PM
Post: #5
 Dan C Member Posts: 98 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
I just tried this on my CASIO fx-4000P, and the result was 0.693147648
I did the square root 9 times, and multiplied with 512 in the last step.
03-05-2022, 05:53 PM
Post: #6
 Dan C Member Posts: 98 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
(03-05-2022 05:32 PM)Dan C Wrote:  I just tried this on my CASIO fx-4000P, and the result was 0.693147648
I did the square root 9 times, and multiplied with 512 in the last step.

I just had to try this on my CASIO fx-4500P also, and the result was the same,
0.693147648.

Now i just have to try this on my FX-602P also
03-05-2022, 08:02 PM
Post: #7
 Dan C Member Posts: 98 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
On my FX-602P this gives 0.693147371
hmm, different from my fx-4000p and fx-4500p.
What does this tell of the FX-602P?
Less accurate?
03-05-2022, 08:29 PM
Post: #8
 Namir Senior Member Posts: 823 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
The extended limit for ln(x) does very well with values of v = 0.001 or less. Very impressive!

Namir
03-05-2022, 08:39 PM
Post: #9
 Thomas Klemm Senior Member Posts: 1,687 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
We experience cancellation when calculating $$\varepsilon = \sqrt[n]{x} - 1$$.
The program was executed on free42 with 34 decimal digits of precision.
Thus this effect is not noticed.

It is a trade-off between coming close to $$1$$ and how many terms to use in the Taylor series.

For a calculator like the HP-15C with 10 digits using $$n = 2^6$$ might be the best choice.

For $$x = 2$$ we end up with the following sequence of iterated square roots:

2.000000000
1.414213562
1.189207115
1.090507733
1.044273783
1.021897149
1.010889286

This leaves us with $$\varepsilon = 0.010889286$$ and thus 5 terms of the Taylor series should be enough.
We notice the cancellation since we have now only 8 significant digits left.

Here's the program for the HP-15C or similar calculators:
Code:
   001 {          11 } √x̅    002 {          11 } √x̅    003 {          11 } √x̅    004 {          11 } √x̅    005 {          11 } √x̅    006 {          11 } √x̅    007 {           1 } 1    008 {          30 } −    009 {          36 } ENTER    010 {          36 } ENTER    011 {          36 } ENTER    012 {           5 } 5    013 {          10 } ÷    014 {           4 } 4    015 {          15 } 1/x    016 {          34 } x↔y    017 {          30 } −    018 {          20 } ×    019 {           3 } 3    020 {          15 } 1/x    021 {          34 } x↔y    022 {          30 } −    023 {          20 } ×    024 {           2 } 2    025 {          15 } 1/x    026 {          34 } x↔y    027 {          30 } −    028 {          20 } ×    029 {           1 } 1    030 {          34 } x↔y    031 {          30 } −    032 {          20 } ×    033 {           6 } 6    034 {           4 } 4    035 {          20 } ×

For $$x = 2$$ we get:

0.6931471776

While for $$\ln(2)$$ we get:

0.6931471806

This leaves us with a difference of:

0.0000000030
03-06-2022, 12:04 AM
Post: #10
 Albert Chan Senior Member Posts: 1,896 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
(03-03-2022 11:20 PM)Thomas Klemm Wrote:  We can compare this to your solution:

$$\sqrt{2(1 + \varepsilon) - 1} - 1$$

The Taylor series agrees for the first two terms:

$$\varepsilon - \frac{\varepsilon^2}{2} + \frac{\varepsilon^3}{2} - \frac{5 \varepsilon^4}{8} + \frac{7 \varepsilon^5}{8} + \mathcal{O}(\varepsilon^6)$$

log1p(ε) = 2*atanh(y) = 2*(y+y^3/3+y^5/5+...) ≈ 2y, where y=ε/(2+ε)

2ε/(2+ε) is simpler than √(1+2ε)-1, and twice as accurate, when ε is tiny.

CAS> series(2*ε/(2+ε),ε) // estimate for log1p(ε)

ε - 1/2*ε^2 + 1/4*ε^3 - 1/8*ε^4 + 1/16*ε^5 + ε^6*order_size(ε)

CAS> (1/3-1/2)/(1/3-1/4)

-2.

Is is funny sqrt approximation formula give better estimate for log1p, than sqrt itself.

2*x/(x+2)
03-06-2022, 12:50 AM
Post: #11
 Albert Chan Senior Member Posts: 1,896 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
(03-05-2022 08:39 PM)Thomas Klemm Wrote:  For $$x = 2$$ we end up with the following sequence of iterated square roots:

2.000000000
1.414213562
1.189207115
1.090507733
1.044273783
1.021897149
1.010889286

This leaves us with $$\varepsilon = 0.010889286$$ and thus 5 terms of the Taylor series should be enough.
We notice the cancellation since we have now only 8 significant digits left.

It is more work, but we could improve accuracy by "pulling out" 1.

CAS> f(x) := x/(sqrt(1+x)+1) // == sqrt(1+x)-1
CAS> 1. // x-1 = 2-1 = 1
CAS> f(Ans)
0.414213562373
0.189207115003
9.05077326653e−2
4.42737824274e−2
2.18971486541e−2
1.08892860517e−2

The recursive formula, log1p(x) = 2*log1p(x/(sqrt(1+x)+1)), is similar to atan formula
(05-31-2021 09:51 PM)Albert Chan Wrote:  $$\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{1+x^2}+1} \right)$$
03-06-2022, 09:33 AM
Post: #12
 Thomas Klemm Senior Member Posts: 1,687 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
I was curious to try Pauli's idea and use:

$$\log\left(\frac{1+\varepsilon}{1-\varepsilon}\right) = 2 \left( \varepsilon + \frac{\varepsilon^3}{3} + \frac{\varepsilon^5}{5} + \mathcal{O}(\varepsilon^7) \right)$$

We set:

$$\frac{1+\varepsilon}{1 - \varepsilon} = x$$

$$\varepsilon = \frac{x - 1}{x + 1}$$

Here's the corresponding program for the HP-15C:
Code:
   001 {          11 } √x̅    002 {          11 } √x̅    003 {          11 } √x̅    004 {          11 } √x̅    005 {          11 } √x̅    006 {          11 } √x̅    007 {          36 } ENTER    008 {          36 } ENTER    009 {           1 } 1    010 {          30 } −    011 {          34 } x↔y    012 {           1 } 1    013 {          40 } +    014 {          10 } ÷    015 {          36 } ENTER    016 {          36 } ENTER    017 {          36 } ENTER    018 {           5 } 5    019 {          10 } ÷    020 {          20 } ×    021 {           3 } 3    022 {          15 } 1/x    023 {          40 } +    024 {          20 } ×    025 {          20 } ×    026 {           1 } 1    027 {          40 } +    028 {          20 } ×    029 {           1 } 1    030 {           2 } 2    031 {           8 } 8    032 {          20 } ×

The computation of $$\varepsilon$$ is now a bit more complicated, but we have fewer terms to compute.
In the end, the program is shorter.

The result for $$x = 2$$ is:

0.6931471773

Compared to the previous result, it is only slightly off in the last digit.
03-09-2022, 07:47 PM
Post: #13
 Albert Chan Senior Member Posts: 1,896 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
Instead of reducing argument, √...√x, we can blow it up, and get log from AGM

From Series[EllipticK[1-(4x)^2], {x,0,4}], and ignored imaginery parts:

K(m=1-(4/x)^2) = log(x) + (log(x)-1)*(4/x^2) + O(1/x^4)

lua> x = 2^14
lua> a, b = 1, 4/x
lua> c = b/x
lua> for i=1,6 do a,b = (a+b)/2, sqrt(a*b); print(a,b) end
Code:
0.5001220703125         0.015625 0.25787353515625        0.08839913658307309 0.17313633586966154     0.15098277337311447 0.162059554621388       0.16168056210089243 0.1618700583611402      0.16186994744240293 0.16187000290177156     0.16187000290176207

With 6 sqrt, (a,b) already close enough so that (a+b)/2 ≈ AGM

lua> k = pi/(a+b)
lua> k / 14
0.6931471898242749
lua> (k+c)*(1-c) / 14
0.6931471805599455
lua> log(x) / 14 -- = log(2)
0.6931471805599453

Note: if c=4/x^2 small enough, below machine epsilon, log(x) = k, no need for correction.
03-10-2022, 05:18 AM
Post: #14
 Thomas Klemm Senior Member Posts: 1,687 Joined: Dec 2013
RE: Ln(x) using repeated square root extraction
In the initial thread An old logarithm algorithm the "pages 33 and 34 of the manual for the Texas Instruments SR-10" are quoted.

However in this manual for the Texas Instruments electronic slide rule calculator SR-10 we find on page 30:
Quote:Logarithmic and Exponential Function

$$\ln a = \left[ \left( - \frac{3}{5} b^2 + 1 \right)^{-1} \times 5 + 4 \right] \frac{2b}{9}$$

$$0.7 < a < 1.6$$

where $$b = \frac{a - 1}{a + 1} = \left( a + 1 \right)^{-1} \times (-2) + 1$$

This expression yields values with an error less than 0.0003% over the range of a from 0.7 to 1.6.

The first 3 terms of the Taylor series of this expression agree with those of $$\log\left(\frac{1+\varepsilon}{1-\varepsilon}\right)$$:

$$\frac{2\varepsilon}{9} \left( 4 + \frac{5}{1 - \frac{3\varepsilon^2}{5}} \right) = 2 \varepsilon + \frac{2\varepsilon^3}{3} + \frac{2\varepsilon^5}{5} + \frac{6\varepsilon^7}{25} + \frac{18\varepsilon^9}{125} + \mathcal{O}(\varepsilon^{11})$$

Again a program for the HP-15C:
Code:
   001 {          11 } √x̅    002 {          11 } √x̅    003 {          11 } √x̅    004 {          11 } √x̅    005 {          11 } √x̅    006 {          11 } √x̅    007 {           1 } 1    008 {          40 } +    009 {           2 } 2    010 {          34 } x↔y    011 {          10 } ÷    012 {           1 } 1    013 {          34 } x↔y    014 {          30 } −    015 {          36 } ENTER    016 {       43 11 } g x²    017 {           3 } 3    018 {          20 } ×    019 {           5 } 5    020 {          10 } ÷    021 {           1 } 1    022 {          34 } x↔y    023 {          30 } −    024 {           5 } 5    025 {          34 } x↔y    026 {          10 } ÷    027 {           4 } 4    028 {          40 } +    029 {          20 } ×    030 {           1 } 1    031 {           2 } 2    032 {           8 } 8    033 {          20 } ×    034 {           9 } 9    035 {          10 } ÷

The result for $$x = 2$$ is:

0.6931471786
03-10-2022, 03:17 PM (This post was last modified: 03-12-2022 03:49 PM by Albert Chan.)
Post: #15
 Albert Chan Senior Member Posts: 1,896 Joined: Jul 2018
RE: Ln(x) using repeated square root extraction
(03-10-2022 05:18 AM)Thomas Klemm Wrote:  The first 3 terms of the Taylor series of this expression agree with those of $$\log\left(\frac{1+\varepsilon}{1-\varepsilon}\right)$$:

$$\frac{2\varepsilon}{9} \left( 4 + \frac{5}{1 - \frac{3\varepsilon^2}{5}} \right) = 2 \varepsilon + \frac{2\varepsilon^3}{3} + \frac{2\varepsilon^5}{5} + \frac{6\varepsilon^7}{25} + \frac{18\varepsilon^9}{125} + \mathcal{O}(\varepsilon^{11})$$
FYI, approximation formula based from atanh(ε) pade approximation:

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

Let d = 1 - 3/5*ε^2, we have ε^2/3 = 5/9*(1-d)

ε + (ε^3/3)/d = ε*(1 + 5/9*(1-d)/d) = ε/9*(9 + 5*(1/d-1)) = ε/9*(4 + 5/d)

ln((1+ε)/(1-ε)) = 2*atanh(ε) ≥ 2ε/9*(4 + 5/(1-3/5*ε^2))

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

We doubled accuracy if we use above formula instead (Bonus, less code steps)
Again, taylor series of ln((1+ε)/(1-ε)), using above formula for ln:

XCAS> g := (x-1)/sqrt(x)
XCAS> f := g - g/(27/10+24/g^2)
XCAS> series(f(x=(1+ε)/(1-ε)),ε,0,10)

2*ε + 2/3*ε^3 + 2/5*ε^5 + 123/400*ε^7 + 3217/12000*ε^9 + O(ε^11)

XCAS> (c - 6/25) / (c - 123/400) | c=2/7.

-2.09836065574
 « Next Oldest | Next Newest »

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