Author Topic: graph of the gamma function (Euler-Mascheroni)  (Read 1580 times)

0 Members and 1 Guest are viewing this topic.

BSpinoza

• Newbie
• Posts: 80
graph of the gamma function (Euler-Mascheroni)
« on: August 20, 2021, 01:39:19 AM »
This Programm calculates and draws the graph of the gamma function (Euler-Mascheroni)
by using the Stirling's approximation (or Stirling's formula).
Code: (qb64) [Select]
`' ------------------------' PROGRAM: GAMMA_GRAPH.BAS' ------------------------' This Programm calculates and draws the graph' of the gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula).'     created: April 04, 2017 by Bruno Sch„fer, Losheim am See, Germany' last review: August 19, 2021' ----------------------------------------------------------------------' Some interesting values at points:''      â(1/2) = sqrt(ã)'      â(-1/2) = -2 * sqrt(ã)'      â(-1) = â(-2) = â(-3) = infinity ì' ----------------------------------------------------------------------SCREEN _NEWIMAGE(600, 600, 256)_TITLE "Graph of the GAMMA function"CLS , 15COLOR 0, 15WINDOW (-600, -600)-(600, 600)' ----------------------' help linesFOR y = -500 TO 500 STEP 50    LINE (-500, y)-(500, y), 7 'help lines  YNEXT yFOR x = -500 TO 500 STEP 50    LINE (x, -500)-(x, 500), 7 'help lines  XNEXT x' --------------'   axis ticksFOR y = -500 TO 500 STEP 50    LINE (-6, y)-(8, y), ' ticks y axisNEXT yFOR x = -500 TO 500 STEP 50 'ticks x axis    LINE (x, -8)-(x, 6) ', 7NEXT x' -----------------' axesLINE (0, -500)-(0, 500) 'y-AchseLINE (-500, 0)-(500, 0) 'x-AchseLINE (-500, 500)-(500, 500) 'upper border lineLINE (-500, -500)-(500, -500) 'lower border lineLINE (-500, -500)-(-500, 500) 'left border lineLINE (500, -500)-(500, 500) 'right border line' --------------------' labels y-axisLOCATE 2, 32: PRINT CHR\$(226) + "  function"; 'title_PRINTSTRING (278, 43), " 5"_PRINTSTRING (278, 93), " 4"_PRINTSTRING (278, 143), " 3"_PRINTSTRING (278, 193), " 2"_PRINTSTRING (278, 243), " 1"_PRINTSTRING (270, 343), " -1"_PRINTSTRING (270, 393), " -2"_PRINTSTRING (270, 443), " -3"_PRINTSTRING (270, 493), " -4"_PRINTSTRING (270, 543), " -5"' --------------------' labels x-axisLOCATE 20, 5: PRINT " -5"LOCATE 20, 13: PRINT "-4"LOCATE 20, 19: PRINT "-3"LOCATE 20, 25: PRINT "-2"LOCATE 20, 31: PRINT "-1"LOCATE 20, 37: PRINT "0"LOCATE 20, 44: PRINT "1"LOCATE 20, 50: PRINT "2"LOCATE 20, 57: PRINT "3"LOCATE 20, 63: PRINT "4"LOCATE 20, 69: PRINT "5"' -----------------' Drawing the curveFOR t = -5 TO 5 STEP 0.0005    y = gamma(t) * 100    x = t * 100    IF ABS(y) < 497 THEN        PSET (x, y), 12        DRAW "U2 R2 D2 L2"    END IFNEXT t' -------------- gamma function --------------------------' returns the value of the Gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula)' example x = 1.25'       Keisan (CASIO online calculator):'               gamma (X) = 0.90640247705547707798267128896691800074879192072'       QB64 with _FLOAT - Type precision:'            gamma##(x##) = 0.906402477055477FUNCTION gamma## (x##)    PI = _PI(1)    n = 550    a## = (x## + n)    b## = x##    FOR I = 1 TO (n - 1)        b## = b## * (x## + I)    NEXT I    gamma## = SQR((2 * PI) / a##) * a## ^ a## * EXP(a## * (-1) + (1 / (12 * a##)) - (1 / (360 * a## ^ 3))) / b##END FUNCTION`
"Simplicity is the result of maturity."  Friedrich Schiller

_vince

• Seasoned Forum Regular
• Posts: 348
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #1 on: August 21, 2021, 10:40:35 PM »
nice, what is a gamma function and what is it for?

BSpinoza

• Newbie
• Posts: 80
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #2 on: August 21, 2021, 11:03:25 PM »
"Simplicity is the result of maturity."  Friedrich Schiller

Ashish

• Forum Resident
• Posts: 632
• Never Give Up!
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #3 on: August 22, 2021, 12:16:03 AM »
Nice!!
I remember only one thing about gamma, that is, its relation with factorial.

Gamma(N + 1) = N factorial
if (Me.success) {Me.improve()} else {Me.tryAgain()}

My Projects - https://github.com/AshishKingdom?tab=repositories
OpenGL tutorials - https://ashishkingdom.github.io/OpenGL-Tutorials

david_uwi

• Newbie
• Posts: 75
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #4 on: August 22, 2021, 03:51:20 AM »
Gamma functions are really useful for evaluating integrals.
Integral(exp(a*x^2)) has no (known) algebraic solution and it turns up with regularity in physics and chemistry*.
But the same integral from zero to infinity (between limits) is given by a gamma function.

*for instance in diffusion and in calculation of the entropy of gases.

SMcNeill

• QB64 Developer
• Forum Resident
• Posts: 3657
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #5 on: August 22, 2021, 12:31:05 PM »
Gamma functions are really useful for evaluating integrals.
Integral(exp(a*x^2)) has no (known) algebraic solution and it turns up with regularity in physics and chemistry*.
But the same integral from zero to infinity (between limits) is given by a gamma function.

*for instance in diffusion and in calculation of the entropy of gases.

You guys are talking over my head again with all this integral stuff.  I just thought Gamma was another option on your video screens to brighten and darken the shadows on your display!  Who the heck knew it had anything to do with Buzz Lightyear??

“To Infinity and Beyond!!”
https://github.com/SteveMcNeill/Steve64 — A github collection of all things Steve!

Qwerkey

• Forum Resident
• Posts: 712
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #6 on: August 22, 2021, 12:53:28 PM »
Bruno (@BSpinoza) suggested that this program should be placed in the Library (Maths & Geometry).  From the Wikipedia pages, I was rather sceptical of this abstruse function being of any general use anywhere or of any interest to our members, so it's interesting to hear David's(@david_uwi) comments on its practical applicability.

bplus

• Forum Resident
• Posts: 7446
• b = b + ...
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #7 on: August 22, 2021, 01:30:55 PM »
Bruno (@BSpinoza) suggested that this program should be placed in the Library (Maths & Geometry).  From the Wikipedia pages, I was rather sceptical of this abstruse function being of any general use anywhere or of any interest to our members, so it's interesting to hear David's(@david_uwi) comments on its practical applicability.

A general plotting application would be so much more useful.

bplus

• Forum Resident
• Posts: 7446
• b = b + ...
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #8 on: August 22, 2021, 03:01:45 PM »
The n in line 92 seems arbitrary, why was 550 picked?

Code: (qb64) [Select]
`Print " 99! = 9.3326215443944152681699238856267e+155 " 'according to MS Windows calculatorPrint gamma##(100, 10)Print gamma##(100, 100)Print gamma##(100, 550)' -------------- gamma function --------------------------' returns the value of the Gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula)' example x = 1.25'       Keisan (CASIO online calculator):'               gamma (X) = 0.90640247705547707798267128896691800074879192072'       QB64 with _FLOAT - Type precision:'            gamma##(x##) = 0.906402477055477Function gamma## (x##, n)    PI = _Pi(1)    a## = (x## + n)    b## = x##    For I = 1 To (n - 1)        b## = b## * (x## + I)    Next I    gamma## = Sqr((2 * PI) / a##) * a## ^ a## * Exp(a## * (-1) + (1 / (12 * a##)) - (1 / (360 * a## ^ 3))) / b##End Function`
« Last Edit: August 22, 2021, 03:11:07 PM by bplus »

BSpinoza

• Newbie
• Posts: 80
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #9 on: August 23, 2021, 02:49:22 AM »
Code: (qb64) [Select]
`PRINT " 99! = 9.3326215443944152681699238856267e+155 " 'according to MS Windows calculatorPRINT gamma##(100, 10)PRINT gamma##(100, 100)PRINT gamma##(100, 550)' -------------- gamma function --------------------------' returns the value of the Gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula)' example x = 1.25'       Keisan (CASIO online calculator):'               gamma (X) = 0.90640247705547707798267128896691800074879192072'       QB64 with _FLOAT - Type precision:'            gamma##(x##) = 0.906402477055477FUNCTION gamma## (x##, n) 'Gamma function (Euler-Mascheroni), Stirling's    'n = 550    a## = (x## + n)    b## = x##    FOR I = 1 TO (n - 1)        b## = b## * (x## + I)    NEXT I    gamma## = SQR((2 * _PI(1)) / a##) * a## ^ a## * EXP(a## * (-1) + (1 / (12 * a##)) - (1 / (360 * a## ^ 3))) / b##END FUNCTION`
This delivers the result:

99! = 9.3326215443944152681699238856267e+155
9.332621544393954D+155
9.332621544394392D+155
9.332621544394415D+155

I increased the n (arbitrarly) until the result gives the best agreement with the digital places of the exact data. So n = 550 was ok for me.
"Simplicity is the result of maturity."  Friedrich Schiller

BSpinoza

• Newbie
• Posts: 80
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #10 on: August 23, 2021, 05:05:30 AM »
@bplus:
If you change my code, please do it correct, otherwise will not get the precision, which explains the use of that value "n = 550" (yes you can also use n= 499... it' a little bit arbitrary):

If you won't use _PI(1) directly in the formula, please exchange it with the correct data type.

you should use:

Code: (qb64) [Select]
` PI## = _PI(1)`instead of:
Code: (qb64) [Select]
` PI = _PI(1)`
So your correct version is:

Code: (qb64) [Select]
`PRINT " 99! = 9.3326215443944152681699238856267e+155 " 'according to MS Windows calculatorPRINT gamma##(100, 10)PRINT gamma##(100, 100)PRINT gamma##(100, 550)' -------------- gamma function --------------------------' returns the value of the Gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula)' example x = 1.25'       Keisan (CASIO online calculator):'               gamma (X) = 0.90640247705547707798267128896691800074879192072'       QB64 with _FLOAT - Type precision:'            gamma##(x##) = 0.906402477055477FUNCTION gamma## (x##, n) 'Gamma function (Euler-Mascheroni), Stirling's    PI## = _PI(1)    a## = (x## + n)    b## = x##    FOR I = 1 TO (n - 1)        b## = b## * (x## + I)    NEXT I    gamma## = SQR((2 * PI##) / a##) * a## ^ a## * EXP(a## * (-1) + (1 / (12 * a##)) - (1 / (360 * a## ^ 3))) / b##END FUNCTION`
"Simplicity is the result of maturity."  Friedrich Schiller

bplus

• Forum Resident
• Posts: 7446
• b = b + ...
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #11 on: August 23, 2021, 08:46:12 AM »
@bplus:
If you change my code, please do it correct, otherwise will not get the precision, which explains the use of that value "n = 550" (yes you can also use n= 499... it' a little bit arbitrary):

If you won't use _PI(1) directly in the formula, please exchange it with the correct data type.

you should use:

Code: (qb64) [Select]
` PI## = _PI(1)`instead of:
Code: (qb64) [Select]
` PI = _PI(1)`
So your correct version is:

Code: (qb64) [Select]
`PRINT " 99! = 9.3326215443944152681699238856267e+155 " 'according to MS Windows calculatorPRINT gamma##(100, 10)PRINT gamma##(100, 100)PRINT gamma##(100, 550)' -------------- gamma function --------------------------' returns the value of the Gamma function (Euler-Mascheroni)' by using the Stirling's approximation (or Stirling's formula)' example x = 1.25'       Keisan (CASIO online calculator):'               gamma (X) = 0.90640247705547707798267128896691800074879192072'       QB64 with _FLOAT - Type precision:'            gamma##(x##) = 0.906402477055477FUNCTION gamma## (x##, n) 'Gamma function (Euler-Mascheroni), Stirling's    PI## = _PI(1)    a## = (x## + n)    b## = x##    FOR I = 1 TO (n - 1)        b## = b## * (x## + I)    NEXT I    gamma## = SQR((2 * PI##) / a##) * a## ^ a## * EXP(a## * (-1) + (1 / (12 * a##)) - (1 / (360 * a## ^ 3))) / b##END FUNCTION`

I used PI exactly the same as from your original code, see your line 91. For the function, I just took out the line n = 550 and put n as argument to test different values for it, that was all that I changed.

So I guess you fixed your own code :)
« Last Edit: August 23, 2021, 08:51:43 AM by bplus »

jack

• Seasoned Forum Regular
• Posts: 367
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #12 on: August 23, 2021, 08:55:25 AM »
_Pi will only give Pi to double precision, if you want Pi to _Float precision then don't use _Pi(1)
here's an implementation for small arguments that uses the power series of 1/Gamma(x) and then simply uses recursion
this works well for small arguments

Code: (qb64) [Select]
`Function gamma## (x As _Float)    Dim As Long i    Static As _Float A(28), pi    Dim As _Float x1, ix, fx, y    If x = 1 Or x = 2 Then        gamma## = 1##        Exit Function    End If    pi = 3.1415926535897932384##    A(0) = 1.000000000000000000000F0    A(1) = 5.772156649015328606065F-1    A(2) = -6.558780715202538810770F-1    A(3) = -4.200263503409523552900F-2    A(4) = 1.665386113822914895017F-1    A(5) = -4.219773455554433674820F-2    A(6) = -9.621971527876973562114F-3    A(7) = 7.218943246663099542395F-3    A(8) = -1.165167591859065112113F-3    A(9) = -2.152416741149509728157F-4    A(10) = 1.280502823881161861531F-4    A(11) = -2.013485478078823865568F-5    A(12) = -1.250493482142670657345F-6    A(13) = 1.133027231981695882374F-6    A(14) = -2.056338416977607103450F-7    A(15) = 6.116095104481415817862F-9    A(16) = 5.002007644469222930055F-9    A(17) = -1.181274570487020144588F-9    A(18) = 1.043426711691100510491F-10    A(19) = 7.782263439905071254049F-12    A(20) = -3.696805618642205708187F-12    A(21) = 5.100370287454475979015F-13    A(22) = -2.058326053566506783222F-14    A(23) = -5.348122539423017982370F-15    A(24) = 1.226778628238260790158F-15    A(25) = -1.181259301697458769513F-16    A(26) = 1.186692254751600332579F-18    A(27) = 1.412380655318031781555F-18    A(28) = -2.298745684435370206592F-19    x1 = x - 1##    ix = Int(Abs(x1))    If x1 < 0 Then ix = -ix    fx = x1 - ix    If fx = 0## And x1 < 0## Then        Print "Error, (in GAMMA) singularity encountered"        Exit Function    End If    If fx <> 0 Then        y = A(28)        For j = 27 To 0 Step -1            y = A(j) + y * fx        Next        y = 1## / y        If ix > 0 Then            For i = 1## To ix                y = y * (fx + i)            Next        End If        If ix < 0 Then            x1 = -x1            y = pi * x1 / (gamma##(x1 + 1##) * Sin(pi * x1))        End If    ElseIf ix > 0 Then        y = 1##        For i = 1## To ix            y = y * i        Next    End If    gamma## = yEnd Function`

BSpinoza

• Newbie
• Posts: 80
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #13 on: August 24, 2021, 09:03:44 AM »
@bplus:
Quote
I used PI exactly the same as from your original code, see your line 91. For the function, I just took out the line n = 550 and put n as argument to test different values for it, that was all that I changed.

So I guess you fixed your own code :)

Dear bplus, sorry for that dissent.

If you use your code (posted August 22) and you use my code (the function from the original first post or your version, i've changed in August 23):

You will get different results with about 500 iteration steps:

PRINT gamma##(100, 550) :

results
my Version:         9.332621544394415D +155      -> all 16 decimal places are correct
your version:       9.332621674246341D +155       -> only the first 7 decimal places are correct

The only explanation I have for this different results is the use of

##PI = _PI(1)

PI = _PI(1)

Please check it by yourself!

"Simplicity is the result of maturity."  Friedrich Schiller

bplus

• Forum Resident
• Posts: 7446
• b = b + ...
Re: graph of the gamma function (Euler-Mascheroni)
« Reply #14 on: August 24, 2021, 11:24:41 AM »
Hi @BSpinoza

OK I ran your fix (Aug 23) of "my" code that was copied from your original post and the n = 155 line removed and made a parameter to the function. Fix is namely PI## = ...

The results are better on my system (Widows 10 laptop) and QB64 version as shown in bottom right corner but still not the full 16 digits as you claim.