Accurate x - log(1+x) - 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: Accurate x - log(1+x) (/thread-18330.html) Accurate x - log(1+x) - Albert Chan - 05-05-2022 07:52 PM (04-27-2020 07:54 PM)Albert Chan Wrote:  I tried to re-use my formula, for log(probability of no repetition) ln_P(n,s) := -s/(12*n*(n-s)) - (s + (n-s+0.5)*ln(1-s/n)) 30 people birthday "collision" ≈ -expm1(ln_P(365,30)) = 0.706316242724 Using exact formula: ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿1 - perm(365,30)/365^30 = 0.706316242719 For big n, small x = s/n, we can drop the first term. Second term has a problem, even if we have accurate log1p (s + n*ln(1 - s/n)) / n = x + log(1-x) ≈ x - x = 0 There might be a better way. For now, I use: x + log(1-x) = x - (x + x^2/2 + x^3/3 + x^4/4 + ...) = -(x^2/2 + x^3/3 + x^4/4 + ...) We have correct formula for ln_P(n,s), but numerically very hard to evaluate with big n. It would be nice if we can accurately calculate x - log1p(x) For small x, we can sum taylor/pade series. log1p(x) = 2*atanh(x/(x+2)) Code: function tinyx_sub_log1p(x) --   x - log1p(x)     local y = x/(x+2)       -- = x - 2*atanh(y)     local z = y*y           -- = y*(x+2) - 2*(y + y^3/3 + y^5/5 + ...)     z = z*(1/3 + z*(1/5 + z*(1/7 + z/9 / (1-9/11*z))))     return y * (x-2*z)      -- relerr O(x^11/805376) end y = x/(x+2) ≈ x/(0+2) = x/2 abs_err(z) ≈ (1/13 - 9/121) * y^12 ≈ (4/1573) * (x/2)^12 = x^12/1610752 (x-2z) ≈ x - 2/3*y² ≈ x - 2/3*(x²/4) = x - x²/6 ≈ x relerr = abs_err(2z)/(x-2z) ≈ x^11/805376 relerr(x=0.1000) ≈ 1.24e-17, below Free42 Binary machine epsilon relerr(x=0.0025) ≈ 2.96e-35, below Free42 Decimal machine epsilon For big x, we use "half-angle" formula for log1p Keep reducing log1p argument, until we can use tinyx_sub_log1p log1p(x) = 2 * log1p(√(1+x) - 1) = 2 * log1p(x / (√(1+x) + 1)) Let y=sqrt(1+x), z = (y-1) = x/(y+1) Let f(x) = x - log1p(x) f(x) = z*(y+1) - 2*log1p(z) = z*(y-1) + 2*f(z) = z^2 + 2*f(z) Both RHS terms are non-negative, without cancellation errors. Code: function x_sub_log1p(x)     if abs(x) < 0.1 then return tinyx_sub_log1p(x) end     x = x / (sqrt(1+x)+1)   -- "half-angle" formula     return x_sub_log1p(x)*2 + x*x end RE: Accurate x - log(1+x) - Albert Chan - 05-05-2022 08:18 PM (04-27-2020 07:54 PM)Albert Chan Wrote:  I tried to re-use my formula, for log(probability of no repetition) ln_P(n,s) := -s/(12*n*(n-s)) - (s + (n-s+0.5)*ln(1-s/n)) 30 people birthday "collision" ≈ -expm1(ln_P(365,30)) = 0.706316242724 Using exact formula: ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿1 - perm(365,30)/365^30 = 0.706316242719 Let's try our new x_sub_log1p function. (note: y = x - log1p(x) ≥ 0) Code: function ln_nr(n,s)         -- log(probability of no repetition)     local x = -s/n     local y = x_sub_log1p(x)     return x/(12*(n-s)) + n*y + (s-0.5)*(x-y) end function P_collision(n,s) return -expm1(ln_nr(n,s)) end lua> P_collision(365, 30) 0.7063162427241915 lua> P_collision(2^128, 5e12) 3.6734198463188465e-014 lua> P_collision(2^128, 2.6153e18) 0.009999839905579181 Let's see what happen if we replace with naive implementation ... lua> function x_sub_log1p(x) return x - log1p(x) end lua> P_collision(365, 30) -- still OK, n is not huge 0.7063162427241918 lua> P_collision(2^128, 5e12) -- bad 7.346839692638293e-014 lua> P_collision(2^128, 2.6153e18) -- bad 0.019899683013021148 RE: Accurate x - log(1+x) - Albert Chan - 05-05-2022 08:57 PM From [VA] SRC #010 - Pi Day 2022 Special (03-14-2022 08:03 PM)Valentin Albillo Wrote:  unexpected infinite product for $$\pi$$ announced in the title, the awesome expression:      which beautifully relates $$\pi$$ and e. (03-22-2022 01:01 AM)Albert Chan Wrote:   (03-22-2022 12:14 AM)Valentin Albillo Wrote:  To settle down the question, if someone with access to Mathematica or some other arbitrary-precision software can compute the product for N=100,000 using 100 digits, say, or as many as necessary to ensure full 34 correct digits or more, and post here the resulting value I'd appreciate it. Thanks in advance. >>> from mpmath import * >>> mp.dps = 100 >>> pn = lambda n: exp(nsum(lambda x: 1+log1p(-1/(x*x))*x*x,[2,n]) + 1.5) >>> n = mpf(100000) >>> N = 2*n+1 >>> x = pn(n) >>> print x 3.141608361513791562872866895754895278060325823725833279147116393910631517290786764227775828378244404 lua> n = 100000 lua> s = 0 lua> for x=n,2,-1 do y=x*x; s = s + (1 + log1p(-1/y)*y) end lua> exp(s + 1.5) -- error = -440 ULP 3.141608361513987 lua> s = 0 lua> for x=n,2,-1 do y=x*x; s = s - x_sub_log1p(-1/y)*y end lua> exp(s + 1.5) -- error = 0 ULP 3.1416083615137915 Using x_sub_log1p(x), we get true PN, without ULP error ! RE: Accurate x - log(1+x) - Albert Chan - 05-06-2022 02:03 PM We can use x_sub_log1p(), to get accurate (e^x - 1 - x), for small x y = x - log(1+x) e^y = e^x / (1+x) (e^y - 1) = (e^x - 1 - x) / (1+x) (e^x - 1 - x) = (e^y - 1) * (1+x) Code: function expm1_sub_x(x)     local z = expm1(x_sub_log1p(x))     return z + z*x end lua> x = 1/1024 lua> expm1(x) - x      -- error = 967 ULP 4.769924165351404e-007 lua> expm1_sub_x(x) -- error = -1 ULP 4.769924165352429e-007 Comment: next post expm1_sub(x) take care of big x too. It is also more consistent throughout, func_sub(x) = func(x) - x RE: Accurate x - log(1+x) - Albert Chan - 05-09-2022 12:41 AM Let f(x) = x - log(1+x) f(x) - f(-x) = (x - log(1+x)) - (-x - log(1-x)) = 2*(x - atanh(x)) LHS f(±x) never negative, thus possible massive cancellations. We can go from accurate (atanh(x)-x) to accurate (log(1+x)-x), but not in reverse. Code for RHS should be factored out. Note that I flipped the order, to be consistent throughout, func_sub(x) = func(x) - x I also use better pade approximation for atanh_sub_tiny(x), making it more useful by itself. Code: function atanh_sub_tiny(x)  -- abserr ≈ 256/2760615 * x^15     local y = x*x           -- pade(atanh(x)-x,x,14,8)     return x*y * (5005 - y*(5082 - y*969))                / (15015 - y*(24255 - y*(11025 - y*1225))) end function log1p_sub_tiny(x)  -- abserr ≈ -1/176679360 * (x/(1+x/2))^15     local y = x/(x+2)     return atanh_sub_tiny(y)*2 - x*y    -- approx -x*y ≤ 0 end function log1p_sub(x)       -- Free42 Decimal, 0.2 -> 0.01     if abs(x) < 0.2 then return log1p_sub_tiny(x) end     if x < -.56 or x > .89 then return log1p(x)-x end     x = x / (sqrt(1+x)+1)   -- "half-angle" formula     return log1p_sub(x)*2 - x*x end function expm1_sub(x)     if x > .86 then return expm1(x) - x end     if x < -.7 then return exp(x)-(1+x) end     local z = expm1(-log1p_sub(x))     return z + z*x end function ln_nr(n,s)         -- log(probability of no repetition)     local x = -s/n     local y = log1p_sub(x)     return x/(12*(n-s)) - n*y + (s-0.5)*(x+y) end Example, expression ((1+r)^-n - 1)/r + n lua> r, n = 1e-6, 1e3 lua> expm1(-n*log1p(r))/r + n --> error = 788 ULP 0.5003328749086222 lua> z = log1p_sub(r) lua> (expm1_sub(-n*(z+r)) - n*z)/r --> error = -1 ULP 0.5003328749087098 lua> mapm = require'mapm' lua> r = mapm.new(r) -- exact binary to decimal conversion lua> (((1+r)^-n - 1)/r + n):tofixed(20) 0.50033287490870967825 RE: Accurate x - log(1+x) - Albert Chan - 04-04-2023 11:05 PM I added HP71B implementation of log1p_sub(x), called FNL(X) https://www.hpmuseum.org/forum/thread-19689-post-170860.html#pid170860 Update 2/02/24: FNL recursion part is not necessary see https://www.hpmuseum.org/forum/thread-21235-post-183567.html#pid183567 (02-02-2024 04:27 PM)Albert Chan Wrote:  FNL(X) recursion replaced with direct fast-path: 200 DEF FNL(X) ! = ln(1+X) - X, but more accurate 210 IF X<-.3832 OR X>.5163 THEN FNL=LOGP1(X)-X @ GOTO 250 220 X2=X/(X+2) @ X4=X2*X2 230 X4=X4*(5005-X4*(5082-X4*969))/(15015-X4*(24255-X4*(11025-X4*1225))) 240 FNL=X2*(X4+X4-X) 250 END DEF