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

0 Members and 1 Guest are viewing this topic.

Offline 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 , 15
COLOR 0, 15
WINDOW (-600, -600)-(600, 600)
' ----------------------
' help lines
FOR y = -500 TO 500 STEP 50
    LINE (-500, y)-(500, y), 7 'help lines  Y
NEXT y
FOR x = -500 TO 500 STEP 50
    LINE (x, -500)-(x, 500), 7 'help lines  X
NEXT x
' --------------
'   axis ticks
FOR y = -500 TO 500 STEP 50
    LINE (-6, y)-(8, y), ' ticks y axis
NEXT y
FOR x = -500 TO 500 STEP 50 'ticks x axis
    LINE (x, -8)-(x, 6) ', 7
NEXT x
' -----------------
' axes
LINE (0, -500)-(0, 500) 'y-Achse
LINE (-500, 0)-(500, 0) 'x-Achse
LINE (-500, 500)-(500, 500) 'upper border line
LINE (-500, -500)-(500, -500) 'lower border line
LINE (-500, -500)-(-500, 500) 'left border line
LINE (500, -500)-(500, 500) 'right border line
' --------------------
' labels y-axis
LOCATE 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-axis
LOCATE 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 curve
FOR 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 IF
NEXT 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.906402477055477
FUNCTION 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

Offline _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?

Offline BSpinoza

  • Newbie
  • Posts: 80
"Simplicity is the result of maturity."  Friedrich Schiller

Offline 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

Offline 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.

Offline SMcNeill

  • QB64 Developer
  • Forum Resident
  • Posts: 3657
    • Steve’s QB64 Archive Forum
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!

Offline 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.

Offline 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.

Offline 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 calculator
Print 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.906402477055477
Function 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 »

Offline 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 calculator
PRINT 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.906402477055477

FUNCTION 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

Offline 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 calculator
PRINT 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.906402477055477


FUNCTION 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

Offline 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 calculator
PRINT 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.906402477055477


FUNCTION 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 »

Offline 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## = y
End Function

Offline 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) 

instead of

  PI = _PI(1)

Please check it by yourself!








"Simplicity is the result of maturity."  Friedrich Schiller

Offline 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.

Bspinoza.PNG