Log-Arcsine Algorithm - Printable Version

Log-Arcsine Algorithm - 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: Log-Arcsine Algorithm (/thread-19068.html)

Log-Arcsine Algorithm - Thomas Klemm - 11-02-2022 10:22 PM

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

RE: Log-Arcsine Algorithm - Nihotte(lma) - 11-03-2022 06:32 PM

(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

RE: Log-Arcsine Algorithm - Albert Chan - 11-04-2022 03:29 PM

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

RE: Log-Arcsine Algorithm - Albert Chan - 11-04-2022 04:58 PM

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

RE: Log-Arcsine Algorithm - Thomas Klemm - 11-05-2022 01:27 AM

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

RE: Log-Arcsine Algorithm - Albert Chan - 11-05-2022 10:47 AM

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

RE: Log-Arcsine Algorithm - Thomas Klemm - 11-05-2022 11:18 PM

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

RE: Log-Arcsine Algorithm - Albert Chan - 11-06-2022 12:09 PM

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)

RE: Log-Arcsine Algorithm - Thomas Klemm - 11-06-2022 02:43 PM

Nice. I like where this is heading.