+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: More special function confusion/annoyance (/thread-21050.html)
More special function confusion/annoyance - deSitter - 12-25-2023 06:32 AM
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
RE: More special function confusion/annoyance - rprosperi - 12-25-2023 02:19 PM
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.
RE: More special function confusion/annoyance - Albert Chan - 12-25-2023 04:00 PM
Complex number support for Si(z), Ci(z) had been in XCas long time ago.
https://github.com/geogebra/giac/blob/master/src/giac/cpp/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;
#endCas> si(3+4i) → 6.7479950814-3.49866372113*i
Cas> ci(3+4i) → -3.49575703398-5.17590521518*i
RE: More special function confusion/annoyance - deSitter - 12-25-2023 06:19 PM
(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/master/src/giac/cpp/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
RE: More special function confusion/annoyance - Albert Chan - 12-26-2023 02:11 PM
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;