func.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 (*/

#include <jwsc/config.h>

#include <inttypes.h>
#include <math.h>

#include "jwsc/math/func.h"

double log_gamma_d(double x)
{
    return (x + 0.5) * log(x + 5.5) - x - 5.5 + 
             log(2.5066282746310005 * 
                (
                    1.000000000190015 + 
                    76.18009172947146     /(1 + x) +
                    -86.50532032941677    /(2 + x) +
                    24.01409824083091     /(3 + x) +
                    -1.231739572450155    /(4 + x) +
                    0.1208650973866179e-2 /(5 + x) + 
                    -0.5395239384953e-5   /(6 + x)
                ) / x);
}


double beta_d(double x, double y)
{
    // beta must be computed in log space for large x, y, 
    return exp(log_beta_d(x, y));
}

double mult_beta_d(const double* const x, uint32_t num_elts)
{
    // beta must be computed in log space for large values
    return exp(log_mult_beta_d(x, num_elts));
}

double log_beta_d(double x, double y)
{
    return log_gamma_d(x) + log_gamma_d(y) - log_gamma_d(x + y);
}

double log_mult_beta_d(const double* const x, uint32_t num_elts)
{
    double numerator = 1;
    double denominator = 0;

    uint32_t i;
    for(i = 0; i < num_elts; i++)
    {
        numerator += log_gamma_d(x[i]);
        denominator += x[i];
    }

    denominator = log_gamma_d(denominator);

    return numerator - denominator;
}