Is there a way to calculate tgammaf(50) on a device using cuda?

Question

I am using stirling's approximation to calculate this, but I am having trouble storing the huge answer I get from Stirling's approximation. Is there any good way to store such a big number on the device?

Answer 1

The CUDA standard math library supports both the double-precision function tgamma and the single-precision functiontgammaf. Since Γ(50) is on the order of 10⁶², the result overflows the representable range of the float type and tgammaf (50.0f) thus returns infinity. However, the computation can proceed without overflow by using the double-precision function tgamma and storing the result to a double variable, e.g.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

__global__ void kernel (int a)
{
    float r = tgammaf ((float)a);
    double rr = tgamma ((double)a);
    printf ("tgammaf(%d) = %23.16e\ntgamma(%d)  = %23.16e\n", a, r, a, rr);
}

int main (void)
{
    kernel<<<1,1>>>(50);
    cudaDeviceSynchronize();
    return EXIT_SUCCESS;
}

The result of the above program should look like so:

tgammaf(50) =                     inf
tgamma(50)  =  6.0828186403426752e+62