Gamma function, SinhIntegral, CoshIntegral

[robmio]

Hello, I found that the HP PRIME "Gamma" function does not provide the same results as Wolfram's MATHEMATICA. For instance:
HP PRIME --> Gamma (4/5, -6) --> 294.845140024
MATHEMATICA --> Gamma [4/5, -6] --> 238.757-172.621*i.
For this reason I have written a program to calculate the Incomplete Gamma function, extending it for numeric values that do not work on the HP PRIME. The code has been named “Upper_Inc_Γ”, which needs two other small programs to work well. One program calculates the value of the integral of the hyperbolic sine (Shi (x)), while the other program calculates the value of the integral of the hyperbolic cosine (Chi (x)).

Shi(x) code:

Code:


#cas
// Shi(z) Sinh Integral z ∈ ℂ
Shi(z):=
BEGIN
IF type(z)==DOM_IDENT OR type(z)==DOM_SYMBOLIC
THEN RETURN 1/2*(Ei(z)-Ei(−z))+G_0;
END;
CASE
IF ARG(z)==0 OR ARG(exact(z))==pi THEN 
1/2*(Ei(z)-Ei(−z)) END;
IF 0<ARG(z)<pi THEN 1/2*(Ei(z)-
Ei(−z))-sqrt(-1)*pi/2 END;
DEFAULT 1/2*(Ei(z)-Ei(−z))+sqrt(-1)*pi/2
END;
END;
#end

Chi(x) code:

Code:


#cas
// Chi(z) → Cosh Integral
Chi(z):=
BEGIN
IF type(z)==DOM_IDENT OR type(z)==DOM_SYMBOLIC
THEN RETURN (G_0+Ei(-z)+Ei(z))/2+G_1;
ELSE
IF type(approx(z))==DOM_INT OR
type(approx(z))==DOM_FLOAT AND 
SIGN(approx(z))==1 THEN 
RETURN 1/2*Ei(-z)+1/2*Ei(z);
ELSE
IF type(approx(z))==DOM_INT OR
type(approx(z))==DOM_FLOAT AND
SIGN(approx(z))==-1 THEN
RETURN sqrt(-1)*pi+1/2*Ei(-z)+1/2*Ei(z);
ELSE
IF type(approx(z))==DOM_COMPLEX AND
SIGN(im(approx(z)))==1 THEN
RETURN (sqrt(-1)*pi+Ei(-z)+Ei(z))/2;
ELSE
IF type(approx(z))==DOM_COMPLEX AND
SIGN(im(approx(z)))==-1 THEN
RETURN ((−sqrt(-1))*pi+Ei(-z)+Ei(z))/2;
END;
END;
END;
END;
END;
END;
#end

Upper_Inc_Γ(n,x):

Code:


#cas
Upper_Inc_Γ(a,z):=
BEGIN
LOCAL funz, t, fnz;
// fnz:=(1/(-a)!)*(Shi(z)-Chi(z));
CASE
IF a==0 AND z==0 THEN
RETURN "Error: a==0 & z==0" END;
IF (type(a)==DOM_INT AND a>0) OR 
type(a)==DOM_COMPLEX THEN
RETURN -z^a*integrate(t^(a-1)*e^(-t*z),t
,0,1)+Gamma(a) END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a>0 OR a<−1) THEN
RETURN (Gamma(a,z)-Gamma(a))*(−1)^(a
-1)+Gamma(a) END;
IF type(a)==DOM_INT AND a≤0 THEN
RETURN (−1)^(-a)*((sum((k-1)!/(
(-a)!*(-z)^k),k,1,-a))*exp(-z)+
(1/(-a)!)*(Shi(z)-Chi(z))) END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a<0 AND a>−1) THEN
RETURN (−z^(a)*exp(−z)+(
(Gamma(a+1,z)-Gamma(a+1))*(−1)^(a+1
-1)+Gamma(a+1)))/a END;
DEFAULT
RETURN Gamma(a,z);
END;
END;
#end

It would be nice if xCas implemented the Gamma function - what do you think?

Best regards, Roberto.

[lrdheat]

The Prime answer is the absolute value of the other platform’s answer…

[robmio]

(11-06-2022 07:03 PM)lrdheat Wrote:  The Prime answer is the absolute value of the other platform’s answer…

I hadn't really noticed: why does HP PRIME return the absolute value?

[lrdheat]

I do not know…

[Albert Chan]

(11-06-2022 06:27 PM)robmio Wrote:  HP PRIME --> Gamma (4/5, -6) --> 294.845140024
MATHEMATICA --> Gamma [4/5, -6] --> 238.757-172.621*i.

Mathematica Gamma, converted to HP Prime Gamma
Note: it is not Abs[Gamma[4/5,-6]] ≈ 294.624

Gamma[4/5] + Abs[Gamma[4/5] - Gamma[4/5,-6]]
= Gamma[4/5] + (Gamma[4/5] - Gamma[4/5,-6]) / (-1)^(4/5)
= 294.845140024 ...

HP Prime Gamma, back to Mathematica Gamma:

CAS> gamma(a,x) := when(x<0, [Gamma(a),Gamma(a,x)] * [1+(-1)^a,-(-1)^a], Gamma(a,x))

CAS> gamma(4/5,-6.)      → 238.757077078-172.62130796*i

Quote:I hadn't really noticed: why does HP PRIME return the absolute value?

It is just a guess, but some integral result is more elegant.

CAS> int(e^x^3)      → 1/3*(Gamma(1/3,-x^3) - Gamma(1/3))
CAS> Ans(x=6.)        → 5.96393809188e91

Mathematica:

∫(e^x^3) = -(x Γ(1/3, -x^3))/(3 (-x^3)^(1/3))

∫(e^x^3, x=0..6) ≈ (5.964E91 + 0.7733*i) - (-0.4465 + 0.7733*i) ≈ 5.964E91

[robmio]

(11-07-2022 02:31 PM)Albert Chan Wrote:  

(11-06-2022 06:27 PM)robmio Wrote:  HP PRIME --> Gamma (4/5, -6) --> 294.845140024
MATHEMATICA --> Gamma [4/5, -6] --> 238.757-172.621*i.

Mathematica Gamma, converted to HP Prime Gamma
Note: it is not Abs[Gamma[4/5,-6]] ≈ 294.624

Gamma[4/5] + Abs[Gamma[4/5] - Gamma[4/5,-6]]
= Gamma[4/5] + (Gamma[4/5] - Gamma[4/5,-6]) / (-1)^(4/5)
= 294.845140024 ...

HP Prime Gamma, back to Mathematica Gamma:

CAS> gamma(a,x) := when(x<0, [Gamma(a),Gamma(a,x)] * [1+(-1)^a,-(-1)^a], Gamma(a,x))

CAS> gamma(4/5,-6.)      → 238.757077078-172.62130796*i

Quote:I hadn't really noticed: why does HP PRIME return the absolute value?

It is just a guess, but some integral result is more elegant.

CAS> int(e^x^3)      → 1/3*(Gamma(1/3,-x^3) - Gamma(1/3))
CAS> Ans(x=6.)        → 5.96393809188e91

Mathematica:

∫(e^x^3) = -(x Γ(1/3, -x^3))/(3 (-x^3)^(1/3))

∫(e^x^3, x=0..6) ≈ (5.964E91 + 0.7733*i) - (-0.4465 + 0.7733*i) ≈ 5.964E91

Dear Albert, you are right: "xCas" does not return the absolute value of the "Gamma" function: I noticed it too, just before connecting to the forum.
I would like you to judge my program for calculating the gamma function. I have not yet implemented the program for calculating the Gamma function with lower bound of integration as a complex number.
Best regards, Roberto.

Code:


#cas
Upper_Inc_Γ(a,z):=
BEGIN
LOCAL funz, t,fnz;
// fnz:=(1/(-a)!)*(Shi(z)-Chi(z));
CASE
IF a==0 AND z==0 THEN
RETURN "Error: a==0 & z==0" END;
IF (type(a)==DOM_INT AND a>0) OR 
((type(a)==DOM_COMPLEX) AND 
SIGN(RE(a))==1) THEN
RETURN -z^a*integrate(t^(a-1)*e^(-t*z),t
,0,1)+Gamma(a) END;
IF (type(a)==DOM_COMPLEX AND SIGN(RE
(a))==-1) OR ((type(a)==DOM_RAT AND 
IM(a)≠0) AND SIGN(RE(a))==-1) OR 
((type(a)==DOM_RAT AND IM(a)≠0) AND 
SIGN(RE(a))==+1) THEN
RETURN simplify(integrate(t^(a-1)*
exp(-t),t,z,inf)) END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a>0 OR a<−1) THEN
RETURN (Gamma(a,z)-Gamma(a))*(−1)^(a
-1)+Gamma(a) END;
IF type(a)==DOM_INT AND a≤0 THEN
RETURN (−1)^(-a)*((sum((k-1)!/(
(-a)!*(-z)^k),k,1,-a))*exp(-z)+
(1/(-a)!)*(Shi(z)-Chi(z))) END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a<0 AND a>−1) THEN
RETURN (−z^(a)*exp(−z)+(
(Gamma(a+1,z)-Gamma(a+1))*(−1)^(a+1
-1)+Gamma(a+1)))/a END;
DEFAULT
RETURN Gamma(a,z);
END;
END;
#end

[robmio]

Good morning everyone,
since I have found that the latest version of my program about the Gamma function does not return a correct result when the arguments are made up of complex numbers, I made a change.
For safety, I also attach the program for the integral hyperbolic sine, for the integral hyperbolic cosine, and for the "Pochhammer symbol":

Shi(z):

Code:


#cas
// Shi(z) Sinh Integral z ∈ ℂ
Shi(z):=
BEGIN
IF type(z)==DOM_IDENT OR type(z)==DOM_SYMBOLIC
THEN RETURN 1/2*(Ei(z)-Ei(−z))+G_0;
END;
CASE
IF ARG(z)==0 OR ARG(exact(z))==pi THEN 
1/2*(Ei(z)-Ei(−z)) END;
IF 0<ARG(z)<pi THEN 1/2*(Ei(z)-
Ei(−z))-sqrt(-1)*pi/2 END;
DEFAULT 1/2*(Ei(z)-Ei(−z))+sqrt(-1)*pi/2
END;
END;
#end

Chi(z):

Code:


#cas
// Chi(z) → Cosh Integral
Chi(z):=
BEGIN
IF type(z)==DOM_IDENT OR type(z)==DOM_SYMBOLIC
THEN RETURN (G_0+Ei(-z)+Ei(z))/2+G_1;
ELSE
IF type(approx(z))==DOM_INT OR
type(approx(z))==DOM_FLOAT AND 
SIGN(approx(z))==1 THEN 
RETURN 1/2*Ei(-z)+1/2*Ei(z);
ELSE
IF type(approx(z))==DOM_INT OR
type(approx(z))==DOM_FLOAT AND
SIGN(approx(z))==-1 THEN
RETURN sqrt(-1)*pi+1/2*Ei(-z)+1/2*Ei(z);
ELSE
IF type(approx(z))==DOM_COMPLEX AND
SIGN(im(approx(z)))==1 THEN
RETURN (sqrt(-1)*pi+Ei(-z)+Ei(z))/2;
ELSE
IF type(approx(z))==DOM_COMPLEX AND
SIGN(im(approx(z)))==-1 THEN
RETURN ((−sqrt(-1))*pi+Ei(-z)+Ei(z))/2;
END;
END;
END;
END;
END;
END;
#end

Pochhammer(a,n):

Code:


#cas
Pochhammer(a,n):=
BEGIN
IF n==0 THEN
RETURN 1;
ELSE
RETURN product((a+j),j,0,n-1);
END;
END;
#end

Upper_Inc_Γ(a,z):

Code:


#cas
Upper_Inc_Γ(a,z):=
BEGIN
LOCAL funz, t,fnz;
// fnz:=(1/(-a)!)*(Shi(z)-Chi(z));
CASE
IF a==0 AND z==0 THEN
RETURN "Error: a==0 & z==0" END;
IF (type(a)==DOM_INT AND a>0) THEN
RETURN -z^a*integrate(t^(a-1)*e^(-t*z),t
,0,1)+Gamma(a) END;
IF (IM(a)≠0 AND type(a)==DOM_RAT)
 OR type(a)==DOM_COMPLEX THEN
RETURN subst(simplify(Gamma(a)-z^a*
e^(-z)*sum(z^k/(Pochhammer(a,k+1)),k,
0,n)),n,100.);
END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a>0 OR a<−1) THEN
RETURN (Gamma(a,z)-Gamma(a))*(−1)^(a
-1)+Gamma(a) END;
IF type(a)==DOM_INT AND a≤0 THEN
RETURN (−1)^(-a)*((sum((k-1)!/(
(-a)!*(-z)^k),k,1,-a))*exp(-z)+
(1/(-a)!)*(Shi(z)-Chi(z))) END;
IF (SIGN(z)==-1 AND type(a)==DOM_RAT)
 AND (a<0 AND a>−1) THEN
RETURN (−z^(a)*exp(−z)+(
(Gamma(a+1,z)-Gamma(a+1))*(−1)^(a+1
-1)+Gamma(a+1)))/a END;
DEFAULT
RETURN Gamma(a,z);
END;
END;
#end