← Back to branch summary

~showard314/ubuntu/karmic/r-base/remove_start_comments

« back to all changes in this revision

Viewing changes to src/nmath/dnchisq.c

  • Committer: Bazaar Package Importer
  • Author(s): Dirk Eddelbuettel
  • Date: 2009-01-19 12:40:24 UTC
  • mfrom: (5.1.4 sid)
  • Revision ID: james.westby@ubuntu.com-20090119124024-abxsf4e0y7713w9m
Tags: 2.8.1-2
http://bugs.debian.org/512324
debian/control: Add another Build-Depends: exclusion for the 
'kfreebsd-i386 kfreebsd-amd64 hurd-i386' architecture to openjdk-6-jdk.
Thanks to Petr Salinger for the heads-up.               (Closes: 512324)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nmath/dnchisq.c
2
2
 *  Mathlib : A C Library of Special Functions
3
3
 *  Copyright (C) 1998 Ross Ihaka
4
4
 *  Copyright (C) 2000-8 The R Development Core Team
5
 
 *  Copyright (C) 2004   The R Foundation
 
5
 *  Copyright (C) 2004-8 The R Foundation
6
6
 *
7
7
 *  This program is free software; you can redistribute it and/or modify
8
8
 *  it under the terms of the GNU General Public License as published by
31
31
{
32
32
    const static double eps = 5e-15;
33
33
 
34
 
    double i, ncp2, q, mid, dfmid, imax, errorbound;
 
34
    double i, ncp2, q, mid, dfmid, imax;
35
35
    LDOUBLE sum, term;
36
36
 
37
37
#ifdef IEEE_754
48
48
        return ML_POSINF;
49
49
    if(ncp == 0)
50
50
        return dchisq(x, df, give_log);
 
51
    if(x == ML_POSINF) return R_D__0;
51
52
 
52
53
    ncp2 = 0.5 * ncp;
53
54
 
54
55
    /* find max element of sum */
55
56
    imax = ceil((-(2+df) +sqrt((2-df) * (2-df) + 4 * ncp * x))/4);
56
57
    if (imax < 0) imax = 0;
57
 
    dfmid = df + 2 * imax;
58
 
    mid = dpois_raw(imax, ncp2, FALSE) * dchisq(x, dfmid, FALSE);
 
58
    if(R_FINITE(imax)) {
 
59
        dfmid = df + 2 * imax;
 
60
        mid = dpois_raw(imax, ncp2, FALSE) * dchisq(x, dfmid, FALSE);
 
61
    } else /* imax = Inf */
 
62
        mid = 0;
 
63
 
 
64
    if(mid == 0) {
 
65
        /* underflow to 0 -- maybe numerically correct; maybe can be more accurate,
 
66
         * particularly when  give_log = TRUE */
 
67
        /* Use  central-chisq approximation formula when appropriate;
 
68
         * ((FIXME: the optimal cutoff also depends on (x,df);  use always here? )) */
 
69
        if(give_log || ncp > 1000.) {
 
70
            double nl = df + ncp, ic = nl/(nl + ncp);/* = "1/(1+b)" Abramowitz & St.*/
 
71
            return dchisq(x*ic, nl*ic, give_log);
 
72
        } else
 
73
            return R_D__0;
 
74
    }
59
75
 
60
76
    sum = mid;
 
77
 
 
78
    /* errorbound := term * q / (1-q)  now subsumed in while() / if() below: */
 
79
 
61
80
    /* upper tail */
62
 
    term = mid;
63
 
    i = imax;
64
 
    df = dfmid;
 
81
    term = mid; df = dfmid; i = imax;
65
82
    do {
66
83
        i++;
67
84
        q = x * ncp2 / i / df;
68
85
        df += 2;
69
 
        term = q * term;
 
86
        term *= q;
70
87
        sum += term;
71
 
        errorbound = term * q / (1-q);
72
 
    } while (errorbound > eps || q >= 1);
 
88
    } while (q >= 1 || term * q > (1-q)*eps);
73
89
    /* lower tail */
74
 
    term = mid;
75
 
    df = dfmid;
76
 
    i = imax;
 
90
    term = mid; df = dfmid; i = imax;
77
91
    while ( i ){
78
92
        df -= 2;
79
93
        q = i * df / x / ncp2;
80
94
        i--;
81
 
        term = q * term;
 
95
        term *= q;
82
96
        sum += term;
83
 
        errorbound = term * q / (1-q);
84
 
        if (errorbound <= eps && q < 1) break;
 
97
        if (q < 1 && term * q <= (1-q)*eps) break;
85
98
    }
86
99
    return R_D_val(sum);
87
100
}