The R Project SVN R

Rev

Rev 79203 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 79203 Rev 89909
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 2003-2026 The R Core Team
3 
 *  Copyright (C) 2003-2007     The R Foundation
 4 
 *  Copyright (C) 2003-2007 The R Foundation
4 
 *
 5 
 *
5 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
6 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
7 
 *  the Free Software Foundation; either version 2, or (at your option)
 8 
 *  the Free Software Foundation; either version 2, or (at your option)
8 
 *  any later version.
 9 
 *  any later version.
Line 49... Line 50...
49 
 *  where rN[j] ~ Bin(n, prob[j]) ,  sum_j rN[j] == n,  sum_j prob[j] == 1,
 50 
 *  where rN[j] ~ Bin(n, prob[j]) ,  sum_j rN[j] == n,  sum_j prob[j] == 1,
50 
 */
 51 
 */
51 
{
 52 
{
52 
    int k;
 53 
    int k;
53 
    double pp;
 54 
    double pp;
54 
    LDOUBLE p_tot = 0.;
 - 
 
55 
    /* This calculation is sensitive to exact values, so we try to
 55 
    /* Using Kahan compensated summation for platform-independent accuracy
56 
       ensure that the calculations are as accurate as possible
 56 
       avoids relying on LDOUBLE (long double) which varies across
57 
       so different platforms are more likely to give the same
 57 
       platforms (e.g., 80-bit on x86 vs 64-bit on ARM).		See PR#18693. */
58 
       result. */
 58 
    double p_tot = 0.,
- 
 
 59 
 	   p_comp = 0.; /* Kahan compensation term */
59 
 
 60 
 
60 
#ifdef MATHLIB_STANDALONE
 61 
#ifdef MATHLIB_STANDALONE
61 
    if (K < 1) { ML_WARNING(ME_DOMAIN, "rmultinom"); return;}
 62 
    if (K < 1) { ML_WARNING(ME_DOMAIN, "rmultinom"); return;}
62 
    if (n < 0)  ML_WARN_ret_NAN(0);
 63 
    if (n < 0)  ML_WARN_ret_NAN(0);
63 
#else
 64 
#else
Line 69... Line 70...
69 
     *       Could make loop one shorter and drop that check !
 70 
     *       Could make loop one shorter and drop that check !
70 
     */
 71 
     */
71 
    for(k = 0; k < K; k++) {
 72 
    for(k = 0; k < K; k++) {
72 
	pp = prob[k];
 73 
	pp = prob[k];
73 
	if (!R_FINITE(pp) || pp < 0. || pp > 1.) ML_WARN_ret_NAN(k);
 74 
	if (!R_FINITE(pp) || pp < 0. || pp > 1.) ML_WARN_ret_NAN(k);
- 
 
 75 
	/* Kahan summation: p_tot += pp with compensation */
- 
 
 76 
	double y = pp - p_comp,
- 
 
 77 
	    t = p_tot + y;
- 
 
 78 
	p_comp = (t - p_tot) - y;
74 
	p_tot += pp;
 79 
	p_tot = t;
75 
	rN[k] = 0;
 80 
	rN[k] = 0;
76 
    }
 81 
    }
77 
    if(fabs((double)(p_tot - 1.)) > 1e-7)
 82 
    if(fabs(p_tot - 1.) > 1e-7)
78 
	MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"),
 83 
	MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"), p_tot);
79 
		      (double) p_tot);
 - 
 
80 
    if (n == 0) return;
 84 
    if(n == 0 ||
81 
    if (K == 1 && p_tot == 0.) return;/* trivial border case: do as rbinom */
 85 
       (K == 1 && p_tot == 0.)) /* trivial border case: do as rbinom */
- 
 
 86 
	return;
82 
 
 87 
 
83 
    /* Generate the first K-1 obs. via binomials */
 88 
    /* Generate the first K-1 obs. via binomials */
84 
 
 89 
 
85 
    for(k = 0; k < K-1; k++) { /* (p_tot, n) are for "remaining binomial" */
 90 
    for(k = 0; k < K-1; k++) { /* (p_tot, n) are for "remaining binomial" */
86 
	if(prob[k] != 0.) {
 91 
	if(prob[k] != 0.) {
87 
	    pp = (double)(prob[k] / p_tot);
 92 
	    pp = prob[k] / p_tot;
88 
	    /* printf("[%d] %.17f\n", k+1, pp); */
 93 
	    /* printf("[%d] %.17f\n", k+1, pp); */
89 
	    rN[k] = ((pp < 1.) ? (int) rbinom((double) n,  pp) :
 94 
	    rN[k] = ((pp < 1.) ? (int) rbinom((double) n,  pp) :
90 
		     /*>= 1; > 1 happens because of rounding */
 95 
		     /*>= 1; > 1 happens because of rounding */
91 
		     n);
 96 
		     n);
92 
	    n -= rN[k];
 97 
	    n -= rN[k];
93 
	}
 98 
	}
94 
	else rN[k] = 0;
 99 
	else rN[k] = 0;
95 
	if(n <= 0) /* we have all*/ return;
 100 
	if(n <= 0) /* we have all*/ return;
96 
	p_tot -= prob[k]; /* i.e. = sum(prob[(k+1):K]) */
 101 
	/* Kahan subtraction: p_tot -= prob[k] with compensation */
- 
 
 102 
	double y = -prob[k] - p_comp,
- 
 
 103 
	    t = p_tot + y;
- 
 
 104 
	p_comp = (t - p_tot) - y;
- 
 
 105 
	p_tot = t;
97 
    }
 106 
    }
98 
    rN[K-1] = n;
 107 
    rN[K-1] = n;
99 
    return;
 108 
    return;
100 
}
 109 
}
101 
#undef ML_WARN_ret_NAN
 110 
#undef ML_WARN_ret_NAN