Log-Arcsine Algorithm
11-02-2022, 10:22 PM
Post: #1
 Thomas Klemm Senior Member Posts: 2,078 Joined: Dec 2013
Log-Arcsine Algorithm
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
11-03-2022, 06:32 PM
Post: #2
 Nihotte(lma) Member Posts: 74 Joined: Mar 2020
RE: Log-Arcsine Algorithm
(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

Thanks Thomas ! That's fine !!

Laurent
11-04-2022, 03:29 PM
Post: #3
 Albert Chan Senior Member Posts: 2,621 Joined: Jul 2018
RE: Log-Arcsine Algorithm
(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
s0  = x = sin(θ)

The program then iterate for s1 = 2*sin(θ/2), s2 = 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
...
11-04-2022, 04:58 PM (This post was last modified: 11-05-2022 01:16 PM by Albert Chan.)
Post: #4
 Albert Chan Senior Member Posts: 2,621 Joined: Jul 2018
RE: Log-Arcsine Algorithm
(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
s0 = (x - 1/x)/2 = sinh(y)

then

s1 = 2*sinh(y/2)
s2 = 4*sinh(y/4) ... --> s = y = ln(x),

sin(θ) = sinh(θ*i)/i, this explained why Log-Arcsine Algorithm have same recurrence relation.
11-05-2022, 01:27 AM
Post: #5
 Thomas Klemm Senior Member Posts: 2,078 Joined: Dec 2013
RE: Log-Arcsine Algorithm
With a minor change we can also calculate the inverse hyperbolic sine function.

$$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}
11-05-2022, 10:47 AM
Post: #6
 Albert Chan Senior Member Posts: 2,621 Joined: Jul 2018
RE: Log-Arcsine Algorithm
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)))
11-05-2022, 11:18 PM
Post: #7
 Thomas Klemm Senior Member Posts: 2,078 Joined: Dec 2013
RE: Log-Arcsine Algorithm
(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.
11-06-2022, 12:09 PM
Post: #8
 Albert Chan Senior Member Posts: 2,621 Joined: Jul 2018
RE: Log-Arcsine Algorithm
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)
11-06-2022, 02:43 PM
Post: #9
 Thomas Klemm Senior Member Posts: 2,078 Joined: Dec 2013
RE: Log-Arcsine Algorithm
Nice. I like where this is heading.
 « Next Oldest | Next Newest »

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