XProfan

Deutsch

p.specht

Jakob Bernoulli's Vorschlag einer kontinuierlichen Funktion statt der nur für ganzzahlige positive Werte definierten Fakultät wurde von Leonard Euler 1729 aufgegriffen und gipfelte in einer sehr komplizierten, im Negativen komplexwertigen Funktion, der 'Gamma-Funktion'. Diese kann bis dato nur mit erheblichem Aufwand mehr oder weniger genau numerisch angenähert berechnet werden.

Daß diese Funktion aber nicht die einzige Möglichkeit darstellt, vernünftige Zwischenwerte zwischen natürlichzahligen Fakultäten zu berechnen, wurde einige Jahrzehnte später klar: Mathematiker wie etwa Hadamard machten Vorschläge, um die störenden Unendlichkeitsstellen und Diskontinuitäten auf der negativen Seite der Zahlenachse zu beseitigen.

Regelmäßig, und das bis in unserer Tage, werden immer wieder Versuche unternommen, die sperrige Gammafunktion endlich durch elegantere Varianten zu ersetzen. Ein Vorschlag von Herrn (p.t.) Peter Luschny aus 2006 beispielsweise schlägt eine besonders elegante und glatte Funktionskurve mit vielen interessanten mathematischen Eigenschaften vor. Wie stets gefordert, hat auch diese 'Luschny-L'-Funktion die Eigenschaft, an ganzzahligen Eingabewerten die entsprechenden, auch schulüblich ausmultiplizierbaren Fakultätswerte anzunehmen. Die Berechnung erfolgt unter Anwendung weiterer, jedoch gut erforschter zahlentheoretischer Funktionen. Anbei eine XProfan-Demoversion des entsprechenden, auf seiner Homepage veröffentlichten, C-Sharp-Programmes.

Hinweis: Das Urheberrecht liegt vollständig bei Herrn P. Luschny! Ein Hinweis auf seine Homepage findet sich im Programmtext.

Kompilieren Markieren Separieren

WindowTitle "Die Luschny-L-Funktion - glatter und eleganter als das Bernoulli/Euler-Gamma"
' C#-Quelle: Appx.4 of https://www.luschny.de/math/factorial/hadamard/HadamardsGammaFunction.html
' Achtung: Die vollen Urheberrechte liegen allein bei Herrn (p.t.) Peter Luschny!
'(DT) Demo-Translation to XProfan11.2a in 2014-10 by P.Specht, Wien (Austria); Ohne Gewähr!
WindowStyle 24:Window 0,0-%maxx,%maxy-40:font 2:set("decimals",14)
var xh&=width(%hwnd)/2:var yh&=height(%hwnd)/2
line 0,yh&-2*xh&,yh&:line xh&,0-xh&,2*yh&
usepen 0,2,255
moveto 0,yh&-LFactorial(&Loop/10)

whileloop -100,100

    print " ";format$("%g",&Loop/10),tab(10);format$("%g",LFactorial(&Loop/10));mkstr$(" ",42-%pos)
    case %csrlin>42:locate 1,1
    lineto xh&+5*&Loop,yh&-10*LFactorial(&Loop/10)

endwhile

beep
waitinput
end

proc Psi :parameters x!

    declare psi!,piCotpix!,result!,C!,r!
    case (x!<=0) and (x!=Round(x!,0)):return -9999999

    if x!<0' reflection formula

        psi! = Psi(1-x!)
        piCotpix! = -Pi()/Tan(-Pi()*x!)
        return psi! - piCotpix!

    endif

    case x!<=val("1e-6"):return -0.57721566490153286 + 1.6449340668482264365 * x! - 1/x!
    result! = 0
    C! = 12

    while x!<C!

        result!=result!-1/x!
        x!=x!+1

    endwhile

    r! = 1/x!
    result!=result!+ Ln(x!) - 0.5 * r!
    r! = r!*r!
    r!=r!*(1/12-(r!*(1/120-(r!*(1/252-(r!*(1/240-(r!*1/132))))))) )
    return result!-r!

endproc'Psi

proc LnFactorial :parameters z!

    '//const double a0 = 1.0 / 12.0;
    '//const double a1 = 1.0 / 30.0;
    '//const double a2 = 53.0 / 210.0;
    '//const double a3 = 195.0 / 371.0;
    '//const double a4 = 22999.0 / 22737.0;
    '//const double a5 = 29944523.0 / 19733142.0;
    '//const double a6 = 109535241009.0 / 48264275462.0;
    var a0! = 0.0833333333333333333333333;
    var a1! = 0.0333333333333333333333333;
    var a2! = 0.252380952380952380952381;
    var a3! = 0.525606469002695417789757;
    var a4! = 1.01152306812684171174737;
    var a5! = 1.51747364915328739842849;
    var a6! = 2.26948897420495996090915;
    var sqrt2Pi! = 0.91893853320467274;
    z! = z! + 1
    return sqrt2Pi! + (Z! - 0.5) * Ln(Z!) - Z! +  \
    a0!/(Z!+a1!/(Z!+a2!/(Z!+a3!/(Z!+a4!/(Z!+a5!/(Z!+a6!/Z!))))))

endproc

proc Factorial :parameters x!

    case (x!<0) and (x!=Round(x!,0)):return -999999999
    case x!=0:return 1
    var y! = x!
    var p! = 1

    while y!<8:p!=p!*y!:y!=y!+1:endwhile

        var r! = Exp(LnFactorial(y!))

        if x!<8:r!=(r!*x!)/(p!*y!):endif

            return r!

        endproc

        proc LFactorial :parameters x!

            case x!=0:return 0.5
            var y! = if(x!<0,-x!*0.5,x!*0.5)
            var G! = y! * (Psi(y!+0.5) - Psi(y!)) - 0.5
            case x!<0:return G!/Factorial(-x!)
            y! = Pi()*x!
            var S! = Sin(y!)/y!
            return (1 - S! * G!) * Factorial(x!)

        endproc