HP71B IBOUND fooled

[Albert Chan]

(05-02-2022 01:42 AM)Albert Chan Wrote:  50 Y=EXP(-X*X/4) @ Y=Y/2/COS(ASIN(Y)/3+PI/6)
60 DEF FNU(U)=U*(1-U)/SQR(-LN(U*U*(3-2*U)))
70 T=TIME @ DISP INTEGRAL(0,Y,P,FNU(IVAR))*3/SQR(PI),TIME-T

What if x is big ? (upper integral limit, u0 is tiny. Or, 1-u ≈ 1)

\(\displaystyle Q(x)
= {3 \over \sqrt{\pi}} \int_0^{u_0} \frac{u\;(1-u)\;du}{\sqrt{-\log(u^2\;(3-2u))}}
≈ {3 \over \sqrt{2\;\pi}} \int_0^{u_0} \frac{u\;du}{\sqrt{-\log(u)}}
= 3 \; Q\left(\sqrt{-4\;\log(u_0)}\;\right)
\)

XCas> y := exp(-x*x/4)
XCas> u0 := y/2 / cos(asin(y)/3 + pi/6)
XCas> series(2*sqrt(-log(u0)), x=inf, polynom)      → \(\displaystyle x + \frac{\ln\left(3\right)}{x}\)

For big x, we have Q(x) / Q(x + ln(3)/x) ≈ 3
Perhaps we can replace constant 3 by other multiplier M ?

\(\displaystyle M =
\frac{Q(x)}{Q \left(x+\frac{k}{x} \right)}
= \frac{erfc \left({x \over \sqrt2}\right)/2} {erfc \left({x \over \sqrt2}
+ {k \over \sqrt2 \; x} \right)/2}
= e^k \left[ 1
+ \frac{k\;(k+2)}{2 x^2}
+ \frac{k\;(k^3+4k^2-16)}{8x^4} \;+\; ... \right]\)

Example: find y such that Q(y) = Q(x = 5) / 2

lua> x, m = 5, log(2)
lua> k = m
lua> x + k/x -- rough y
5.138629436111989
lua> k = m - log1p(k*(k+2)/(2*x^2))
lua> k = m - log1p(k*(k+2)/(2*x^2))
lua> x + k/x
5.131772468987268
lua> k = m - log1p(k*(k+2)/(2*x^2) + k*(k^3+4*k^2-16)/(8*x^4))
lua> k = m - log1p(k*(k+2)/(2*x^2) + k*(k^3+4*k^2-16)/(8*x^4))
lua> x + k/x
5.132077991427134

lua> -icdf(cdf(-x)/2) -- actual y
5.132018332044298

Comment: For k correction, we do not need accurate log1p()
Even simple log1p(ε) ≈ ε may suffice.

Or, we can revert series for k

lua> t1 = - m*(m+2)/(2*x^2)
lua> t2 = m*(m^2+4*m+6)/(2*x^4)
lua> t3 = - m*(15*m^3+92*m^2+252*m+360)/(24*x^6)
lua> k = t3 + t2 + t1 + k
lua> x + k/x
5.131972797544282

To improve y estimate further, apply Aitken Delta^2 process

lua> k = k - t3*t3 / (t3 - t2)
lua> x + k/x
5.132010307411825