Post Reply 
(12C) log1p function
01-31-2019, 01:10 AM (This post was last modified: 02-01-2019 11:47 AM by Albert Chan.)
Post: #1
(12C) log1p function
This trick is from John D Cook's Tricks for computing log(1+X)

I picked the formula from the comment section instead.
The result is slightly more accurate, and no worry about divide by zero.

log1p(X) = ln(1 + X) - (X + 1 - 1 - X) / (X + 1)

Example: log1p(X = 1.23456789e-6)

1.23456789 [EEX] 6 [CHS] [Enter] ; X value

[Enter] 1 + 1 - [X<>Y] - [CHS]     ; epsilon => -4.3211e-10
[Lst-X] 1 + ÷        ; correction to epsilon => -4.321094e-10
[Lst-X] [LN]          ; LN(1+X) => 1.234999237e-6
+                         ; log1p(X) => 1.234567128e-6 (all digits correct)
Find all posts by this user
Quote this message in a reply
01-31-2019, 04:08 AM
Post: #2
RE: (12C) log1p function
From HP-15C Advanced Functions Handbook
Level 1: Correctly Rounded, or Nearly So

[Image: attachment.php?aid=6851]

This trick is used e.g. in Accurate TVM for HP 35s and was adapted for the HP-15C and the HP-34C.

Code:
01 ENTER
02 ENTER
03 EEX
04 +
05 LN
06 X<>Y
07 LSTx
08 EEX
09 -
10 ÷
11 *

Example:

1.23456789 EEX 6 CHS
R/S
0.000001235

CLEAR PREFIX
1234567127

Cheers
Thomas


Attached File(s) Thumbnail(s)
   
Find all posts by this user
Quote this message in a reply
01-31-2019, 05:48 AM
Post: #3
RE: (12C) log1p function
This is how the WP 34S does this computation. One of Kahan's clever algorithms I believe.

- Pauli
Find all posts by this user
Quote this message in a reply
01-31-2019, 08:45 AM (This post was last modified: 01-31-2019 09:03 AM by Dieter.)
Post: #4
RE: (12C) log1p function
(01-31-2019 01:10 AM)Albert Chan Wrote:  This trick is from John D Cook's Tricks for computing log(1+X)

I picked the formula from the comment section instead.

What about this one:
http://www.hpmuseum.org/forum/thread-1072.html

But the key sequence you posted is not correct.
There must be a second [ENTER] at the beginning and the end it's [LSTx] [LN] (without the "1").

Dieter
Find all posts by this user
Quote this message in a reply
01-31-2019, 09:42 AM
Post: #5
RE: (12C) log1p function
Hi, Dieter

Typo fixed. Thanks.

A (quite) accurate ln1+x function, or "how close can you get" part II thread had the same formula. Smile
Was there a part I ?

My rpn.exe confirmed +correction way much better, and avoided divide-by-0 test.

Quote:>rpn
1.23456789e-3 s
1+ log s1 ? # log1p value
0.001233806437707756434884160067124232513609
r 1+ =-10 s2 ? log s3 ? # 1+X, LN(1+X) value
1.001234568
0.001233806547572121411586143127131058379588

- ?10 # correction wanted
-1.09864365E-10

r2 1- r $ - r2 / ?10 # correction = -(X+1-1-X)/(X+1)
-1.09864365E-10

r3 r x r2 1- / r3 - ?10 # correction = log1p(X) * X / (1+X-1) - log1p(X)
-1.099321546E-10
Find all posts by this user
Quote this message in a reply
01-31-2019, 09:43 PM (This post was last modified: 01-31-2019 09:46 PM by Dieter.)
Post: #6
RE: (12C) log1p function
(01-31-2019 09:42 AM)Albert Chan Wrote:  Typo fixed. Thanks.

The second "Enter" at the beginning is still missing. Without it the code will not work.

(01-31-2019 09:42 AM)Albert Chan Wrote:  A (quite) accurate ln1+x function, or "how close can you get" part II thread had the same formula. Smile
Was there a part I ?

Yes, but that was about Bernoulli numbers:
http://www.hpmuseum.org/forum/thread-870.html

But I also posted a similar method for the e^x–1 function:
http://www.hpmuseum.org/forum/thread-5508.html

(01-31-2019 09:42 AM)Albert Chan Wrote:  My rpn.exe confirmed +correction way much better, and avoided divide-by-0 test.

Quote:>rpn
1.23456789e-3 s
1+ log s1 ? # log1p value
0.001233806437707756434884160067124232513609
r 1+ =-10 s2 ? log s3 ? # 1+X, LN(1+X) value
1.001234568
0.001233806547572121411586143127131058379588

- ?10 # correction wanted
-1.09864365E-10

r2 1- r $ - r2 / ?10 # correction = -(X+1-1-X)/(X+1)
-1.09864365E-10

r3 r x r2 1- / r3 - ?10 # correction = log1p(X) * X / (1+X-1) - log1p(X)
-1.099321546E-10

Please excuse me, Albert, but this is one more case where I see a lot of numbers and cryptic code, but I don't understand anything. Who may understand the meaning of "r 1+ =-10 s2 ? log s3 ? # 1+X, LN(1+X) value" ?

Dieter
Find all posts by this user
Quote this message in a reply
01-31-2019, 10:55 PM (This post was last modified: 02-01-2019 01:02 PM by Albert Chan.)
Post: #7
RE: (12C) log1p function
Hi, Dieter

2 Enters are there, but on 2 lines. First one after value entered.

Regarding the cryptic commands, it was for my own toy rpn program.
s for save memory, default to 0
r for recall memory, default to 0
= is round to digits, minus sign meant only for this number, otherwise globally
# is starting comment, and ignored until end-of-line
? is display top of stack, default round to globally set digits (default 40)

Code:
RPN.EXE command             HP-12C keys (approx.)

1.23456789e-3 s             1.23456789 EEX 3 CHS STO 0
1+ log s1 ?                 1 + LN STO 1 <display>
r 1+ =-10 s2 ?              RCL 0 1 + <round-to-10-digits> STO 2 <display>
log s3 ?                    LN STO 3 <display>
- ?10                       - <display 10 digits>

r2 1- r $ - r2 / ?10        RCL 2 1 - RCL 0 [X<>Y] -  RCL 2 / <display 10 digits>
r3 r x r2 1- / r3 - ?10     RCL 3 RCL 0 * RCL 2 1 - / RCL 3 - <display 10 digits>
Find all posts by this user
Quote this message in a reply
01-31-2019, 11:13 PM
Post: #8
RE: (12C) log1p function
.
Hi, Dieter:

(01-31-2019 09:43 PM)Dieter Wrote:  [...]I don't understand anything. Who may understand the meaning of "r 1+ =-10 s2 ? log s3 ? # 1+X, LN(1+X) value" ?

Probably the same people who understands "UNROT OVER SWAP * DUP OBJ\-> DROP2 / 4. ROLL" ... Smile (actual code, not made up)

Regards.
V.
.

  
All My Articles & other Materials here:  Valentin Albillo's HP Collection
 
Visit this user's website Find all posts by this user
Quote this message in a reply
02-01-2019, 07:07 AM
Post: #9
RE: (12C) log1p function
(01-31-2019 01:10 AM)Albert Chan Wrote:  log1p(X) = ln(1 + X) + (X + 1 - 1 - X) / (X + 1)

That should be

ln(1 + X) - (X + 1 - 1 - X) / (X + 1)

as the [CHS] in your instructions indicate, and then it's the same formula as Dieter's.

Cheers, Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
02-01-2019, 12:30 PM (This post was last modified: 02-03-2019 02:28 AM by Albert Chan.)
Post: #10
RE: (12C) log1p function
Hi Werner

Typo fixed. Thanks

Y = 1+X, approximately
eps = -(Y - 1 - X)
Z = eps / Y ; correction to eps

correction = LN(1+X, exactly) - LN(1+X, approximately)
= LN(Y + eps) - LN(Y)
= LN(1 + eps/Y)
= LN(1 + Z)
= Z - Z²/2 + Z³/3 ...

-> log1p(X) ~ LN(Y) - (Y-1-X)/Y
Find all posts by this user
Quote this message in a reply
02-02-2019, 12:28 PM
Post: #11
RE: (12C) log1p function
(01-31-2019 11:13 PM)Valentin Albillo Wrote:  Probably the same people who understands "UNROT OVER SWAP * DUP OBJ\-> DROP2 / 4. ROLL" ... Smile (actual code, not made up)

That's from DavidM's solution in your recent thread [VA] SRC#003- New Year 2019 Special:
(01-17-2019 06:04 PM)DavidM Wrote:  A similar approach using Thomas' first optimization:
Code:
\<<
   [['\pi' 2019. 2019.][1. '\pi' 2019.][1. 1. '\pi']] \->NUM
   0. DUP 1. \->V3
   0.
   DO
      UNROT
      OVER SWAP *
      DUP OBJ\-> DROP2
      /
      4. ROLL
   UNTIL
      OVER SAME
   END
   NIP NIP
\>>
Result: 12.6389823194
Execution time: about 8.5s on a real 50g
Size: 149 bytes (with embedded matrix build), 65 bytes (initial matrix as argument)

As always the context makes it easier to understand what a program does.
Using a local variable M for the constant matrix might help as well:
Code:
« → M
  « [ 0 0 1 ] 0
    DO
      SWAP
      M SWAP *
      DUP OBJ→ DROP2 /
    UNTIL
      ROT OVER SAME
    END
  »
»

But nah, I didn't understand Albert Chan's cryptic commands until he explained them.

Cheers
Thomas
Find all posts by this user
Quote this message in a reply
02-10-2019, 09:08 AM (This post was last modified: 04-19-2022 01:27 AM by Albert Chan.)
Post: #12
RE: (12C) log1p function
Dieter's asinh formula look very nice: log1p(x) = asinh ( (x²+2x)/(2x+2) ) = asinh(x - x/(2+2/x))

We can create "half angle" version, using symmetrical relation: atanh(sin(z)) = asinh(tan(z))

log1p(x) = 2 atanh( x/(x+2) ) = 2 asinh( x/ √((x+2)² - x²) )

log1p(x) = 2 asinh ( x/√(4x+4) )

Comment: it may be easier to derive relation with Arc SOHCAHTOA method, for hyperbolics

atanhq(x²/(x+2)²)      // TOA, O=x², A=(x+2)²
= asinhq(x²/(4x+4))   // SOH, H = A-O = 4x+4
Find all posts by this user
Quote this message in a reply
02-10-2019, 02:19 PM (This post was last modified: 02-11-2019 02:48 AM by Albert Chan.)
Post: #13
RE: (12C) log1p function
Using relation again: atanh(sin(z)) = asinh(tan(z))
Assume x positive, but close to 0

atanh( 1-x ) = asinh( (1-x) / √(2x-x²) )
atanh( x-1 ) = -atanh( 1-x )

Example, x = 1.2345678e-9, using Casio FX-115MS, showing internal digits:

atanh( 1-x ) = atanh(0.999999998766) = 10.6030760457

asinh( (1-x) / √(2x-x²) ) = asinh(20124.612604) = 10.6028460338 (all digits correct)

Edit: another way, also very good: atanh( 1-x ) = ½ ln(2/x - 1)
Find all posts by this user
Quote this message in a reply
Post Reply 




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