← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/Rmath/src/pnbeta.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-02-06 17:54:29 UTC
  • mfrom: (1.1.3)
  • Revision ID: package-import@ubuntu.com-20130206175429-13br5kqpkfjqdmre
Tags: 0.0.0+20130206.git32ff5759-1
* New upstream snapshot.
* debian/copyright: reflect upstream changes
* debian/rules: update get-orig-source to reflect upstream changes
   + Don't ship nginx
   + Adapt for new configure-random target in deps/Makefile
* Enable build of Tk wrapper.
   + debian/control: add build dependency on tk-dev
   + debian/rules: add tk rule to build-arch
* debian/julia.install: install VERSION and COMMIT files
* no-webrepl.patch: new patch
* Refresh other patches
* Add source override for config.status file under deps/random/

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/Rmath/src/pnbeta.c
1
1
/*
2
 
 *  Copyright (C) 2000-2008 The R Development Core Team
 
2
 *  Copyright (C) 2000-2012 The R Core Team
3
3
 *
4
4
 *  Algorithm AS 226 Appl. Statist. (1987) Vol. 36, No. 2
 
5
 *  by Russell V. Lenth
5
6
 *  Incorporates modification AS R84 from AS Vol. 39, pp311-2, 1990
6
7
 *                        and AS R95 from AS Vol. 44, pp551-2, 1995
 
8
 *  by H. Frick and Min Long Lam.
7
9
 *  original (C) Royal Statistical Society 1987, 1990, 1995
8
10
 *
9
11
 *  Returns the cumulative probability of x for the non-central
28
30
    const int    itrmax = 10000;  /* 100 is not enough for pf(ncp=200)
29
31
                                     see PR#11277 */
30
32
 
31
 
    double a0, ax, lbeta, c, errbd, temp, x0, tmp_c;
 
33
    double a0, lbeta, c, errbd, x0, temp, tmp_c;
32
34
    int j, ierr;
33
35
 
34
 
    long double ans, gx, q, sumq;
 
36
    long double ans, ax, gx, q, sumq;
35
37
 
36
38
    if (ncp < 0. || a <= 0. || b <= 0.) ML_ERR_return_NAN;
37
39
 
59
61
    ans = ax = q * temp;
60
62
 
61
63
        /* recurse over subsequent terms until convergence is achieved */
62
 
    j = x0;
 
64
    j = (int) x0;
63
65
    do {
64
66
        j++;
65
 
        temp -= gx;
 
67
        temp -= (double) gx;
66
68
        gx *= x * (a + b + j - 1.) / (a + j);
67
69
        q *= c / j;
68
70
        sumq -= q;
69
71
        ax = temp * q;
70
72
        ans += ax;
71
 
        errbd = (temp - gx) * sumq;
 
73
        errbd = (double)((temp - gx) * sumq);
72
74
    }
73
75
    while (errbd > errmax && j < itrmax + x0);
74
76
 
85
87
        /* o_x  == 1 - x  but maybe more accurate */
86
88
        int lower_tail, int log_p)
87
89
{
88
 
    long double ans= pnbeta_raw(x, o_x, a,b, ncp);
 
90
    long double ans = pnbeta_raw(x, o_x, a,b, ncp);
89
91
 
90
92
    /* return R_DT_val(ans), but we want to warn about cancellation here */
91
 
    if(lower_tail) return log_p ? log(ans) : ans;
 
93
    if (lower_tail) return (double) (log_p ? logl(ans) : ans);
92
94
    else {
93
95
        if(ans > 1 - 1e-10) ML_ERROR(ME_PRECISION, "pnbeta");
94
96
        ans = fmin2(ans, 1.0);  /* Precaution */
95
 
        return log_p ? log1p(-ans) : (1 - ans);
 
97
        return (double) (log_p ? log1p((double)-ans) : (1. - ans));
96
98
    }
97
99
}
98
100