 (12C) log1p function - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (12C) log1p function (/thread-12308.html) (12C) log1p function - Albert Chan - 01-31-2019 01:10 AM 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) RE: (12C) log1p function - Thomas Klemm - 01-31-2019 04:08 AM From HP-15C Advanced Functions Handbook Level 1: Correctly Rounded, or Nearly So 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 RE: (12C) log1p function - Paul Dale - 01-31-2019 05:48 AM This is how the WP 34S does this computation. One of Kahan's clever algorithms I believe. - Pauli RE: (12C) log1p function - Dieter - 01-31-2019 08:45 AM (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 RE: (12C) log1p function - Albert Chan - 01-31-2019 09:42 AM Hi, Dieter Typo fixed. Thanks. A (quite) accurate ln1+x function, or "how close can you get" part II thread had the same formula. 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 RE: (12C) log1p function - Dieter - 01-31-2019 09:43 PM (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. 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 RE: (12C) log1p function - Albert Chan - 01-31-2019 10:55 PM 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  r 1+ =-10 s2 ?              RCL 0 1 +  STO 2  log s3 ?                    LN STO 3  - ?10                       -  r2 1- r \$ - r2 / ?10        RCL 2 1 - RCL 0 [X<>Y] -  RCL 2 /  r3 r x r2 1- / r3 - ?10     RCL 3 RCL 0 * RCL 2 1 - / RCL 3 -  RE: (12C) log1p function - Valentin Albillo - 01-31-2019 11:13 PM . 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" ... (actual code, not made up) Regards. V. . RE: (12C) log1p function - Werner - 02-01-2019 07:07 AM (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 RE: (12C) log1p function - Albert Chan - 02-01-2019 12:30 PM 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 RE: (12C) log1p function - Thomas Klemm - 02-02-2019 12:28 PM (01-31-2019 11:13 PM)Valentin Albillo Wrote:  Probably the same people who understands "UNROT OVER SWAP * DUP OBJ\-> DROP2 / 4. ROLL" ... (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 RE: (12C) log1p function - Albert Chan - 02-10-2019 09:08 AM Dieter's asinh formula look very nice: ﻿﻿log1p(x) = asinh ( (x²+2x)/(2x+2) ) = asinh(x - x/(2+2/x)) We can create another version, using this nicely symmetrical relation: atanh(sin(z)) = asinh(tan(z)) log1p(x) = 2 atan( x/(x+2) ) = 2 asinh( x/ √((x+2)² - x²) ) log1p(x) = 2 asinh ( x/√(4x+4) ) RE: (12C) log1p function - Albert Chan - 02-10-2019 02:19 PM 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)