More special function confusion/annoyance

[deSitter]

The Si(z) function (integral sine) doesn't take complex argument, but the Ci(z) and Ei(z) functions do.

Sigh.

Workarounds. So annoying. It's not as if function theory was invented yesterday. Why has no one put up a sane library of special functions? And they should all take complex arguments when that makes mathematical sense. Who writes this stuff? Do they need my help? Not kidding, I'm retired and ready for it.

-drl

[rprosperi]

Seems like you're the right guy to do the job, so why not just jump in and do it.

If so, other folks that need it will benefit too.

And if not, then you better understand why no one else has done it yet.

[Albert Chan]

Complex number support for Si(z), Ci(z) had been in XCas long time ago.
https://github.com/geogebra/giac/blob/ma...p/usual.cc

Even old XCas 1.5.0-37 (win32) had them.

XCas> Si(3+4i)      → 6.7479950814-3.49866372113*i
XCas> Ci(3+4i)      → -3.49575703398-5.17590521518*i

It is just not ported to HP Prime Cas yet.
In the meantime, this is a temporary workaround, based on mpmath expintegrals.py

Code:

#cas
si(z) :=
BEGIN
LOCAL v:=(Ei(i*z)-Ei(-i*z))*-.5i;
RETURN v - sign(re(z)) * pi/2;
END; 

ci(z) := 
BEGIN
LOCAL v:=(Ei(i*z)+Ei(-i*z))*.5;
LOCAL (x:=re(z)),(y:=im(z));
IF x>0 THEN RETURN v END;
y := x ? sign(y+!y) : sign(y)/2;
RETURN v + y*pi*i;
END; 
#end

Cas> si(3+4i)      → 6.7479950814-3.49866372113*i
Cas> ci(3+4i)      → -3.49575703398-5.17590521518*i

[deSitter]

(12-25-2023 04:00 PM)Albert Chan Wrote:  Complex number support for Si(z), Ci(z) had been in XCas long time ago.
https://github.com/geogebra/giac/blob/ma...p/usual.cc

Even old XCas 1.5.0-37 (win32) had them.

XCas> Si(3+4i)      → 6.7479950814-3.49866372113*i
XCas> Ci(3+4i)      → -3.49575703398-5.17590521518*i

It is just not ported to HP Prime Cas yet.
In the meantime, this is a temporary workaround, based on mpmath expintegrals.py

Code:

#cas
si(z) :=
BEGIN
LOCAL v:=(Ei(i*z)-Ei(-i*z))*-.5i;
RETURN v - sign(im(z)) * pi/2;
END; 

ci(z) := 
BEGIN
LOCAL v:=(Ei(i*z)+Ei(-i*z))*.5;
LOCAL (x:=re(z)),(y:=im(z));
IF x>0 THEN RETURN v END;
y := x ? sign(y+!y) : sign(y)/2;
RETURN v + y*pi*i;
END; 
#end

Cas> si(3+4i)      → 6.7479950814-3.49866372113*i
Cas> ci(3+4i)      → -3.49575703398-5.17590521518*i

Good job, thanks for that. We will get to the mountain top this way.

-drl

[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;