HP71B IBOUND fooled

[Albert Chan]

From Kahan's https://people.eecs.berkeley.edu/~wkahan/MathSand.pdf, page 23

\(\displaystyle Q(x) = \frac{1}{\sqrt{2\pi}} \int_x^\infty \exp\left(-\frac{u^2}{2} \right)\;du \)

To make integral limits finite, we let t = exp(-u*u/2) ⇒ dt = (t) (-u) du

\(\displaystyle Q(x) = \frac{1}{\sqrt{2\pi}}
\int_0^{t_0} \frac{dt}{\sqrt{-2 \ln(t)}} \qquad\) // t₀ = exp(-x²/2) ≤ 1

For HP71B, this integral have problems: t=0, infinite slope; t=1, infinite value.
INTEGRAL internal u-transform is enough to make end-points finite, but not zero.

Kahan's method "B": Let t = sin(s²) ⇒ dt/ds = cos(s²) (2s)
For both t=0 (s=0) or t=1 (s=√(pi/2)), dt/ds = 0

We could also work in reverse, with dt/du = (1-u) (6u) ⇒ t = u² (3-2u)
Unlike Kahan's transformation, getting u end-point involved solving cubics.
Luckily, with 0 ≤ x ≤ 1, only 1 u solution satisfied.

\(\displaystyle \large \int_0^x f(t)\;dt \qquad\) // 0 ≤ x ≤ 1

\(\displaystyle = \int_0^{s_0} f(\; \sin(s^2)\; )\;\cos(s^2)\;(2s)\;ds\qquad\quad\) // s0 = sqrt(asin(x))
\(\displaystyle = \int_0^{u_0} f(3u^2-2u^3)\;(1-u)\;(6u)\;du\qquad\,\) // u0 = (sqrt(x)/2) / cos(asin(sqrt(x))/3 + pi/6))

We could also scale integral to limits 0 to 1, then do u-transform. (u0 = 1)
But, this left scaling factor inside integrand, not as efficient.
Getting u0 involve only 1 time cost.

10 INPUT "x,eps = ? ";X,P
20 Y=SQR(ASIN(EXP(-X*X/2)))
30 DEF FNS(T) @ S=SIN(T*T) @ FNS=T*SQR((S+1)*(S-1)/LN(S)) @ END DEF
40 T=TIME @ DISP INTEGRAL(0,Y,P,FNS(IVAR))/SQR(PI),TIME-T

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
80 GOTO 10

Code:

>RUN
x,eps = ? 0, 1E-3
 .500003367483        .44
 .500002347538        .37
x,eps = ? 0, 1E-6
 .500000000524        1.52
 .500000057737        .7
x,eps = ? 0, 1E-9
 .499999999995        5.79
 .500000004166        1.88

x,eps = ? 10, 1E-3
 7.61991738243E-24    .21
 7.61991738248E-24    .23
x,eps = ? 10, 1E-6
 7.61985318691E-24    .93
 7.61985318696E-24    .67
x,eps = ? 10, 1E-9
 7.61985302414E-24    5.83
 7.61985302419E-24    4.87

x,eps = ? 20, 1E-3
 2.75363048487E-89    .19
 2.75363048486E-89    .21
x,eps = ? 20, 1E-6
 2.75362414276E-89    .62
 2.75362414274E-89    .7
x,eps = ? 20, 1E-9
 2.75362411868E-89    2.83
 2.75362411866E-89    2.25

X,EPS = ? 30, 1E-3
 4.90671903956E-198   .22
 4.90671903951E-198   .24
X,EPS = ? 30, 1E-6
 4.90671417978E-198   .52
 4.90671417974E-198   .42
X,EPS = ? 30, 1E-9
 4.90671392717E-198   2.85
 4.90671392713E-198   2.4

https://en.wikipedia.org/wiki/Q-function