More special function confusion/annoyance

[Albert Chan]

As noted in expintegrals.py si_generic(), si(ε) has cancellation issues.
With arbitrary precision, it switched to hypergeometric series instead.

For 12-digits precision, we only need a few terms.

sin(x)/x = 1 - x^2/3! + x^4/5! - x^6/7! + ...

Integrate above, from 0 to ε, we have:

si(ε) = ε * (1 - ε^2/(3*3!) + ε^4/(5*5!) - ε^6/(7*7!) + ...)
si(ε) = ε * (1 - ε^2/(3*3!) / (1+0.03*ε^2) + 191/8820000*ε^6 - ...)

> solve(191/8820000*ε^6 = 1e−13, ε = 0.01)      → 4.08074261156e−2

Siimplify with identity (x/k)/(1+x) = (1 - 1/(1+x))/k , we have:

si(ε, |ε|<0.04) ≈ ε * (1 - (1 - 1/(1+0.03*ε^2))/0.54)

Updated si(z) code

Code:

si(z) :=
BEGIN
IF abs(z)<.04 THEN RETURN z - (z-z/(1+.03*z*z))/.54 END
RETURN -.5*((Ei(i*z)-Ei(-i*z))*i + sign(re(z))*pi);
END;