# HP Forums

Full Version: Gamma function, SinhIntegral, CoshIntegral
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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.
The Prime answer is the absolute value of the other platform’s answer…
(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?
I do not know…
(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
(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```
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```
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :