Rebol3 Code Examplex


Gamma function

Evaluate or approximate the Gamma function, a continuous extension of the factorial to non-integer values.

Rebol [
    title: "Rosetta code: Gamma function"
    file:  %Gamma_function.r3
    url:   https://rosettacode.org/wiki/Gamma_function
]

gamma: function/with [
    "Compute gamma function via Spouge/Taylor series for 1/Gamma(x)"
    x [number!]
][
    y: x - 1.0                 ;; series is expanded around (x - 1)
    sm: first approx
    foreach n next approx [
        sm: sm * y + n         ;; Horner's method, highest-degree first
    ]
    1.0 / sm                   ;; series approximates 1/Gamma(x)
][
    ;; Coefficients in natural (lowest-degree first) order for readability;
    ;; `reverse` reorders once at load time so the foreach loop above works.
    approx: reverse [
         1.00000000000000000000   0.57721566490153286061  -0.65587807152025388108
        -0.04200263503409523553   0.16653861138229148950  -0.04219773455554433675
        -0.00962197152787697356   0.00721894324666309954  -0.00116516759185906511
        -0.00021524167411495097   0.00012805028238811619  -0.00002013485478078824
        -0.00000125049348214267   0.00000113302723198170  -0.00000020563384169776
         0.00000000611609510448   0.00000000500200764447  -0.00000000118127457049
         0.00000000010434267117   0.00000000000778226344  -0.00000000000369680562
         0.00000000000051003703  -0.00000000000002058326  -0.00000000000000534812
         0.00000000000000122678  -0.00000000000000011813   0.00000000000000000119
         0.00000000000000000141  -0.00000000000000000023   0.00000000000000000002
    ]
]

repeat i 10 [
    printf ["gamma" -3 "/3 = "] [i gamma i / 3]
]

Output:

gamma  1/3 = 2.67893853470775
gamma  2/3 = 1.3541179394264
gamma  3/3 = 1.0
gamma  4/3 = 0.892979511569249
gamma  5/3 = 0.902745292950934
gamma  6/3 = 1.0
gamma  7/3 = 1.190639348759
gamma  8/3 = 1.50457548825154
gamma  9/3 = 1.99999999999397
gamma 10/3 = 2.77815847933857