Gamma function, SinhIntegral, CoshIntegral - Printable Version +- 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: Gamma function, SinhIntegral, CoshIntegral (/thread-19088.html) Gamma function, SinhIntegral, CoshIntegral - robmio - 11-06-2022 06:27 PM 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 00) 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. RE: Gamma function, SinhIntegral, CoshIntegral - lrdheat - 11-06-2022 07:03 PM The Prime answer is the absolute value of the other platform’s answer… RE: Gamma function, SinhIntegral, CoshIntegral - robmio - 11-06-2022 07:24 PM (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? RE: Gamma function, SinhIntegral, CoshIntegral - lrdheat - 11-06-2022 08:45 PM I do not know… RE: Gamma function, SinhIntegral, CoshIntegral - Albert Chan - 11-07-2022 02:31 PM (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 RE: Gamma function, SinhIntegral, CoshIntegral - robmio - 11-07-2022 03:34 PM (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``` RE: Gamma function, SinhIntegral, CoshIntegral - robmio - 11-08-2022 11:49 AM 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 00) 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```