/*
  Mathlib : A C Library of Special Functions
  Copyright (C) 1999-2007  The R Development Core Team

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or (at
  your option) any later version.

  This program is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, a copy is available at
  http://www.r-project.org/Licenses/

  SYNOPSIS

    #include <Rmath.h>
    double dwilcox(double x, double m, double n, int give_log)
    double pwilcox(double x, double m, double n, int lower_tail, int log_p)
    double qwilcox(double x, double m, double n, int lower_tail, int log_p);
    double rwilcox(double m, double n)

  DESCRIPTION

    dwilcox	The density of the Wilcoxon distribution.
    pwilcox	The distribution function of the Wilcoxon distribution.
    qwilcox	The quantile function of the Wilcoxon distribution.
    rwilcox	Random variates from the Wilcoxon distribution.

*/

/* 
   Note: the checks here for R_CheckInterrupt also do stack checking.

   calloc/free are remapped for use in R, so allocation checks are done there.
   freeing is completed by an on.exit action in the R wrappers.
*/

#include "nmath.h"
#include "dpq.h"

#ifndef MATHLIB_STANDALONE
#include <R_ext/Utils.h>
#endif

static double *w; 
static int allocated_m, allocated_n;

static void
w_free()
{
    if(!w) return;

    free((void *) w);
    w = 0;
    allocated_m = allocated_n = 0;
}

void wilcox_free()
{
    w_free();
}


static void
w_init_maybe(int m, int n)
{
    int c= (m*n+1)/2;


    if( w){
        if( m == allocated_m && n == allocated_n ) return;
        else w_free();
    }

    if( !w){
        w = (double *) calloc(c + 1, sizeof(double));
#ifdef MATHLIB_STANDALONE
        if (!w) MATHLIB_ERROR(_("wilcox allocation error %d"), 1);
#endif
        allocated_m = m; allocated_n = n;
    }
}

/*
   This counts the number of choices with statistic = k

   More correctly, it returns the number of choices with statistic = k
   but it also evaluates number of choices with for every statistic k

   This is because the function implements Harding's algorithm.
   (Harding, E.F. (1984): An Efficient, Minimal-storage Procedure
   for Calculating the Mann-Whitney U, Generalized U and Similar
   Distributions, App. Statist., 33, 1-6)
*/
static double
cwilcox(int k, int m, int n)
{
    int c, u, i, j;
    int q, step;

#ifndef MATHLIB_STANDALONE
    R_CheckUserInterrupt();
#endif

    u = m * n;
    if (k < 0 || k > u)
        return 0;

    c = (int)(u / 2);
    if (k > c)
        k = u - k;

    if( w[k] != 0 )
        return w[k];

    for( i=0; i < n+1; ++i)
        w[i]=1.;

    for( i=2; i < m+1; ++i){
        step = n+i;
        q = c;
        while( q-step > -1){
            w[q] = w[q] - w[q-step];
            --q;
        }
        for( j=i; j<c+1; ++j)
            w[j] += w[j-i];
    }

    return w[k];
}

double dwilcox(double x, double m, double n, int give_log)
{
    double d;

#ifdef IEEE_754
    /* NaNs propagated correctly */
    if (ISNAN(x) || ISNAN(m) || ISNAN(n))
        return(x + m + n);
#endif
    m = floor(m + 0.5);
    n = floor(n + 0.5);
    if (m <= 0 || n <= 0)
        ML_ERR_return_NAN;

    if (fabs(x - floor(x + 0.5)) > 1e-7)
        return(R_D__0);

    x = floor(x + 0.5);
    if ((x < 0) || (x > m * n))
        return(R_D__0);

    w_init_maybe(m, n);
    d = give_log ?
        log(cwilcox(x, m, n)) - lchoose(m + n, n) :
        cwilcox(x, m, n)  / choose(m + n, n);

    return(d);
}

/* args have the same meaning as R function pwilcox */
double pwilcox(double q, double m, double n, int lower_tail, int log_p)
{
    int i;
    double c, p;

#ifdef IEEE_754
    if (ISNAN(q) || ISNAN(m) || ISNAN(n))
        return(q + m + n);
#endif
    if (!R_FINITE(m) || !R_FINITE(n))
        ML_ERR_return_NAN;
    m = floor(m + 0.5);
    n = floor(n + 0.5);
    if (m <= 0 || n <= 0)
        ML_ERR_return_NAN;

    q = floor(q + 1e-7);

    if (q < 0.0)
        return(R_DT_0);
    if (q >= m * n)
        return(R_DT_1);

    w_init_maybe(m, n);
    c = choose(m + n, n);
    p = 0;
    /* Use summation of probs over the shorter range */
    if (q <= (m * n / 2)) {
        for (i = q; i > -1; i--)
            p += cwilcox(i, m, n) / c;
    }
    else {
        q = m * n - q;
        for (i = q-1; i > -1; i--)
            p += cwilcox(i, m, n) / c;
        lower_tail = !lower_tail; /* p = 1 - p; */
    }

    return(R_DT_val(p));
} /* pwilcox */

/* x is 'p' in R function qwilcox */

double qwilcox(double x, double m, double n, int lower_tail, int log_p)
{
    double c, p, q;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(m) || ISNAN(n))
        return(x + m + n);
#endif
    if(!R_FINITE(x) || !R_FINITE(m) || !R_FINITE(n))
        ML_ERR_return_NAN;
    R_Q_P01_check(x);

    m = floor(m + 0.5);
    n = floor(n + 0.5);
    if (m <= 0 || n <= 0)
        ML_ERR_return_NAN;

    if (x == R_DT_0)
        return(0);
    if (x == R_DT_1)
        return(m * n);

    if(log_p || !lower_tail)
        x = R_DT_qIv(x); /* lower_tail,non-log "p" */

    w_init_maybe(m, n);
    c = choose(m + n, n);
    p = 0;
    q = 0;
    if (x <= 0.5) {
        x = x - 10 * DBL_EPSILON;
        for (;;) {
            p += cwilcox(q, m, n) / c;
            if (p >= x)
                break;
            q++;
        }
    }
    else {
        x = 1 - x + 10 * DBL_EPSILON;
        for (;;) {
            p += cwilcox(q, m, n) / c;
            if (p > x) {
                q = m * n - q;
                break;
            }
            q++;
        }
    }

    return(q);
}

double rwilcox(double m, double n)
{
    int i, j, k, *x;
    double r;

#ifdef IEEE_754
    /* NaNs propagated correctly */
    if (ISNAN(m) || ISNAN(n))
        return(m + n);
#endif
    m = floor(m + 0.5);
    n = floor(n + 0.5);
    if ((m < 0) || (n < 0))
        ML_ERR_return_NAN;

    if ((m == 0) || (n == 0))
        return(0);

    r = 0.0;
    k = (int) (m + n);
    x = (int *) calloc(k, sizeof(int));
#ifdef MATHLIB_STANDALONE
    if (!x) MATHLIB_ERROR(_("wilcox allocation error %d"), 4);
#endif
    for (i = 0; i < k; i++)
        x[i] = i;
    for (i = 0; i < n; i++) {
        j = floor(k * unif_rand());
        r += x[j];
        x[j] = x[--k];
    }
    free(x);
    return(r - n * (n - 1) / 2);
}