Gamma function, SinhIntegral, CoshIntegral
11-06-2022, 06:27 PM
Post: #1
 robmio Member Posts: 170 Joined: Jan 2020
Gamma function, SinhIntegral, CoshIntegral
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.
11-06-2022, 07:03 PM
Post: #2
 lrdheat Senior Member Posts: 860 Joined: Feb 2014
RE: Gamma function, SinhIntegral, CoshIntegral
The Prime answer is the absolute value of the other platform’s answer…
11-06-2022, 07:24 PM
Post: #3
 robmio Member Posts: 170 Joined: Jan 2020
RE: Gamma function, SinhIntegral, CoshIntegral
(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?
11-06-2022, 08:45 PM
Post: #4
 lrdheat Senior Member Posts: 860 Joined: Feb 2014
RE: Gamma function, SinhIntegral, CoshIntegral
I do not know…
11-07-2022, 02:31 PM
Post: #5
 Albert Chan Senior Member Posts: 2,687 Joined: Jul 2018
RE: Gamma function, SinhIntegral, CoshIntegral
(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, 03:34 PM
Post: #6
 robmio Member Posts: 170 Joined: Jan 2020
RE: Gamma function, SinhIntegral, CoshIntegral
(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
11-08-2022, 11:49 AM
Post: #7
 robmio Member Posts: 170 Joined: Jan 2020
RE: Gamma function, SinhIntegral, CoshIntegral
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)