Log-Arcsine Algorithm

[Thomas Klemm]

Recently I stumbled upon the Log-Arcsine Algorithm in Peter Henrici's book: Computational Analysis with the HP-25 Pocket Calculator. (pp. 207)

The following recurrence relation is used:

\(
\begin{align}
s_{n+1} = s_n \sqrt{\frac{2 s_n}{s_n + s_{n-1}}}
\end{align}
\)

This little program can be used with most HP calculators:

Code:

ENTER
ENTER
+
LASTX
R↑
+
÷
SQRT
×

Example

0 ENTER 1
R/S
R/S
R/S
…

1.41421356237
1.53073372946
1.56072257613
1.56827424527
1.57016557848
1.57063862547
1.57075690057
1.57078647018
1.57079386264
1.57079571076
1.57079617279
1.57079628829
1.57079631717
1.57079632439
1.57079632619
1.57079632664
1.57079632676
1.57079632679
1.57079632679
1.57079632679
…


Let's slap a loop around it and we have yet another program for the next \(\pi\)-day:

Code:

00 { 26-Byte Prgm }
01▸LBL "pi"
02 0
03 1
04▸LBL 00
05 ENTER
06 ENTER
07 +
08 LASTX
09 R↑
10 +
11 ÷
12 SQRT
13 ×
14 X≠Y?
15 GTO 00
16 2
17 ×
18 END

This allows us to calculate both \(\log\) and \(\arcsin\) on an HP-16C:

Code:

001 - 43,22, A  LBL A      
002 -       36  ENTER      
003 -       36  ENTER      
004 -    43 26  1/x        
005 -       40  +          
006 -        2  2          
007 -       10  /          
008 -       34  x<>y       
009 -       36  ENTER      
010 -    43 26  1/x        
011 -       30  -          
012 -        2  2          
013 -       10  /          
014 -       20  x          
015 -    43 36  LSTx       
016 -    22  0  GTO 0      
017 - 43,22, b  LBL B      
018 -       36  ENTER      
019 -       36  ENTER      
020 -       20  x          
021 -        1  1          
022 -       34  x<>y       
023 -       30  -          
024 -    43 25  SQRT       
025 -       34  x<>y       
026 -       20  x          
027 -    43 36  LSTx       
028 - 43,22, 0  LBL 0      
029 -       36  ENTER      
030 -       36  ENTER      
031 -       40  +          
032 -    43 36  LSTx       
033 -    43 33  R^         
034 -       40  +          
035 -       10  /          
036 -    43 25  SQRT       
037 -       20  x          
038 -    43 34  PSE        
039 -    43  0  x!=y       
040 -    22  0  GTO 0

Examples

10
GSB A

2.302585094

0.5
GSB B

0.523598775

Here is the same program for the HP-42S:

Code:

00 { 56-Byte Prgm }
01▸LBL "log"
02 ENTER
03 ENTER
04 1/X
05 +
06 2
07 ÷
08 X<>Y
09 ENTER
10 1/X
11 -
12 2
13 ÷
14 ×
15 LASTX
16 GTO 00
17▸LBL "asin"
18 ENTER
19 X↑2
20 1
21 X<>Y
22 -
23 SQRT
24 X<>Y
25 ×
26 LASTX
27▸LBL 00
28 ENTER
29 ENTER
30 +
31 LASTX
32 R↑
33 +
34 ÷
35 SQRT
36 ×
37 X≠Y?
38 GTO 00
39 END

References

• Computational analysis with the HP-25 pocket calculator (Peter Henrici).pdf
• HP16C Trigonometric Functions (no inverses so far)
• Boldly Going – HP-16C Scientific Functions Part 2
  

[Nihotte(lma)]

(11-02-2022 10:22 PM)Thomas Klemm Wrote:  Recently I stumbled upon the Log-Arcsine Algorithm in Peter Henrici's book: Computational Analysis with the HP-25 Pocket Calculator. (pp. 207)

...

This allows us to calculate both \(\log\) and \(\arcsin\) on an HP-16C:
...

References

• Computational analysis with the HP-25 pocket calculator (Peter Henrici).pdf
• HP16C Trigonometric Functions (no inverses so far)
• Boldly Going – HP-16C Scientific Functions Part 2
  

Thanks Thomas ! That's fine !!

Laurent

[Albert Chan]

(11-02-2022 10:22 PM)Thomas Klemm Wrote:  The following recurrence relation is used:

\(
\begin{align}
s_{n+1} = s_n \sqrt{\frac{2 s_n}{s_n + s_{n-1}}}
\end{align}
\)
...
Example

0 ENTER 1
R/S
R/S
R/S
…

1.41421356237
1.53073372946
1.56072257613
1.56827424527
1.57016557848
1.57063862547
1.57075690057
1.57078647018
1.57079386264
...

Above initial conditions for θ = asin(x):

s_-1 = x*√(1-x*x) = sin(θ)*cos(θ) = sin(2θ)/2
s₀  = x = sin(θ)

The program then iterate for s₁ = 2*sin(θ/2), s₂ = 4*sin(θ/4), ... --> s_∞ = θ = asin(x)

Sequence is identical to asin(x) half-angle formula. Note: asinq(x) = asin(√x)

(03-31-2022 11:07 PM)Albert Chan Wrote:  asinq(x) = 2 * asinq(x/2/(1+sqrt(1-x)))

>>> x = 1 # iterate for asin(√1)
>>> for k in range(1,20):
...            x /= 2*(1+sqrt(1-x))
...            print 2**k * sqrt(x)
...
1.41421356237
1.53073372946
1.56072257613
1.56827424527
1.57016557848
1.57063862547
1.57075690057
1.57078647018
1.57079386264
...

[Albert Chan]

(11-02-2022 10:22 PM)Thomas Klemm Wrote:  The following recurrence relation is used:

\(
\begin{align}
s_{n+1} = s_n \sqrt{\frac{2 s_n}{s_n + s_{n-1}}}
\end{align}
\)

If we let y = ln(x), and

s_-1 = (x² - 1/x²)/4 = sinh(2y)/2
s₀ = (x - 1/x)/2 = sinh(y)

then

s₁ = 2*sinh(y/2)
s₂ = 4*sinh(y/4) ... --> s_∞ = y = ln(x),

sin(θ) = sinh(θ*i)/i, this explained why Log-Arcsine Algorithm have same recurrence relation.

[Thomas Klemm]

With a minor change we can also calculate the inverse hyperbolic sine function.
We just have to start with:

\(
x \sqrt{1 + x^2}
\)

and

\(
x
\)

Example

2 SQRT 1

XEQ 00

0.88137358702

---

To find the recurrence relation we can define the following functions with somewhat suggestive names:

\(
\begin{align}
s(z) &= \frac{z - \frac{1}{z}}{2} \\
\\
c(z) &= \frac{z + \frac{1}{z}}{2} \\
\end{align}
\)

It's easy to show (see below) that:

\(
\begin{align}
2 s(z) c(z) &= s(z^2) \\
\\
2 c(z)^2 &= 1 + c(z^2) \\
\end{align}
\)

From the latter we conclude:

\(
\begin{align}
c(z) = \sqrt{\frac{1 + c(z^2)}{2}}
\end{align}
\)

Using the substitution \(z \to z^{\tfrac{1}{2}}\) we get:

\(
\begin{align}
c(z^{\frac{1}{2}}) &= \sqrt{\frac{1 + c(z)}{2}} \\
\\
s(z)
&= 2 \, s(z^{\frac{1}{2}}) \, c(z^{\frac{1}{2}}) \\
\\
&= 2 \, s(z^{\frac{1}{2}}) \, \sqrt{\frac{1 + c(z)}{2}} \\
\\
&= 2 \, s(z^{\frac{1}{2}}) \, \sqrt{\frac{s(z) + s(z)c(z)}{2s(z)}} \\
\\
&= 2 \, s(z^{\frac{1}{2}}) \, \sqrt{\frac{s(z) + \frac{1}{2}s(z^2)}{2s(z)}} \\
\end{align}
\)

Or then after rearranging it:

\(
\begin{align}
2 \, s(z^{\frac{1}{2}}) = s(z) \, \sqrt{\frac{2s(z)}{s(z) + \frac{1}{2}s(z^2)}}
\end{align}
\)

With:

\(
\begin{align}
\ell(h) = \frac{1}{2h}(x^h - x^{-h}) = \frac{s(x^h)}{h}
\end{align}
\)

Or then with \(z = x^h = x^{2^{-n}}\) we have:

\(
\begin{align}
\ell_n = \ell(2^{-n}) = 2^n s(x^{2^{-n}}) = 2^n s(z) = s_n
\end{align}
\)

This leads to the recurrence relation:

\(
\begin{align}
s_{n+1}
&= \ell(2^{-n-1}) \\
&= 2^{n+1} s(x^{2^{-n-1}}) \\
&= 2^n \cdot 2 s(z^{\tfrac{1}{2}}) \\
\\
&= 2^n s(z) \sqrt{\frac{2s(z)}{s(z) + \tfrac{1}{2}s(z^2)}} \\
\\
&= 2^n s(z) \sqrt{\frac{2 \cdot 2^n s(z)}{2^n s(z) + \tfrac{1}{2} \cdot 2^n s(z^2)}} \\
\\
&= s_n \sqrt{\frac{2 s_n}{s_n + 2^{n-1} s(z^2)}} \\
\\
&= s_n \sqrt{\frac{2 s_n}{s_n + s_{n-1}}} \\
\end{align}
\)

We can plug in \(z = e^x\) or \(z = e^{ix}\) to get related formulas for \(\sinh(x)\) and \(\sin(x)\):

\(
\begin{align}
2 \sinh(\frac{1}{2} z) &= \sinh(z) \sqrt{\frac{2\sinh(z)}{\sinh(z) + \frac{1}{2}\sinh(2z)}} \\
\\
2 \sin(\frac{1}{2} z) &= \sin(z) \sqrt{\frac{2\sin(z)}{\sin(z) + \frac{1}{2}\sin(2z)}} \\
\end{align}
\)

That's due to the fact that:

\(
\begin{align}
s(e^x) = \frac{e^x - e^{-x}}{2} = \sinh(x) \\
\\
s(e^{ix}) = \frac{e^{ix} - e^{-ix}}{2} = \sin(x) \\
\end{align}
\)

---

\(
\begin{align}
2 s(z) c(z)
&= 2 \frac{z - \frac{1}{z}}{2} \frac{z + \frac{1}{z}}{2} \\
&= \frac{\left(z - \frac{1}{z} \right) \left(z + \frac{1}{z} \right)}{2} \\
&= \frac{z^2 - \frac{1}{z^2}}{2} \\
&= s(z^2) \\
\\
2 c(z)^2
&= 2 \left( \frac{z + \frac{1}{z}}{2} \right)^2 \\
&= \frac{\left(z + \frac{1}{z} \right)^2}{2} \\
&= \frac{z^2 + 2 + \frac{1}{z^2}}{2} \\
&= 1 + \frac{z^2 + \frac{1}{z^2}}{2} \\
&= 1 + c(z^2) \\
\end{align}
\)

[Albert Chan]

It may be easier to get recurrence relation, directly from cosines.

Let θ = x/2^(n+1), s(n) = 2^n * sin(x/2^n)

cos(θ) = sin(2θ)/(2*sin(θ)) = s(n)/s(n+1)
cos(2θ) = sin(4θ)/(2*sin(2θ)) = s(n-1)/s(n)

cos(θ) = √((1+cos(2θ))/2)

s(n)/s(n+1) = √((1+s(n-1)/s(n))/2) --> s(n+1) = s(n) * √(2*s(n) / (s(n) + s(n-1)))

[Thomas Klemm]

(11-05-2022 10:47 AM)Albert Chan Wrote:  It may be easier

I find it interesting that the recurrence relation is not dependent on trigonometric formulas.

[Albert Chan]

Here is a table of initial conditions, for all inverse trigonometric / hyperbolic / log functions

s(n+1) = s(n) * √(2*s(n) / (s(n) + s(n-1)))      → s(∞) = func(x)

Code:

asin    x*√(1-x²)       x
asinh   x*√(1+x²)       x

acos    x*√(1-x²)       √(1-x²)
acosh   x*√(x²-1)       √(x²-1)

atan    x/(1+x²)        x/√(1+x²)
atanh   x/(1-x²)        x/√(1-x²)

log     (x-1/x)/2       (x-1)/√(x)
log1p   x*(1+x/2)/(1+x) x/√(1+x)

[Thomas Klemm]

Nice. I like where this is heading.