| Line 1... |
Line 1... |
| 1 |
/*
|
1 |
/*
|
| 2 |
* R : A Computer Language for Statistical Data Analysis
|
2 |
* R : A Computer Language for Statistical Data Analysis
|
| 3 |
* Copyright (C) 2001-2014 R Core Team
|
3 |
* Copyright (C) 2001-2026 R Core Team
|
| 4 |
*
|
4 |
*
|
| 5 |
* This program is free software; you can redistribute it and/or modify
|
5 |
* 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
|
6 |
* it under the terms of the GNU General Public License as published by
|
| 7 |
* the Free Software Foundation; either version 2 of the License, or
|
7 |
* the Free Software Foundation; either version 2 of the License, or
|
| 8 |
* (at your option) any later version.
|
8 |
* (at your option) any later version.
|
| Line 37... |
Line 37... |
| 37 |
T=.5*10^(-NSIG).
|
37 |
T=.5*10^(-NSIG).
|
| 38 |
ENTEN = 10 ^ K, where K is the largest int such that
|
38 |
ENTEN = 10 ^ K, where K is the largest int such that
|
| 39 |
ENTEN is machine-representable in working precision
|
39 |
ENTEN is machine-representable in working precision
|
| 40 |
ENSIG = 10 ^ NSIG
|
40 |
ENSIG = 10 ^ NSIG
|
| 41 |
RTNSIG = 10 ^ (-K) for the smallest int K such that K >= NSIG/4
|
41 |
RTNSIG = 10 ^ (-K) for the smallest int K such that K >= NSIG/4
|
| 42 |
ENMTEN = Smallest ABS(X) such that X/4 does not underflow
|
42 |
ENMTEN = Smallest ABS(X) such that X/4 does not underflow [ ~= 4 * DBL_MIN ]
|
| 43 |
XINF = Largest positive machine number; approximately beta ^ maxexp
|
43 |
XINF = Largest positive machine number; approximately beta ^ maxexp
|
| 44 |
== DBL_MAX (defined in #include <float.h>)
|
44 |
== DBL_MAX (defined in #include <float.h>)
|
| 45 |
SQXMIN = Square root of beta ^ minexp = sqrt(DBL_MIN)
|
45 |
SQXMIN = Square root of beta ^ minexp = sqrt(DBL_MIN)
|
| 46 |
|
46 |
|
| 47 |
EPS = The smallest positive floating-point number such that 1.0+EPS > 1.0
|
47 |
EPS = The smallest positive floating-point number such that 1.0+EPS > 1.0
|
| 48 |
= beta ^ (-p) == DBL_EPSILON
|
48 |
= beta ^ (-p) == DBL_EPSILON
|
| 49 |
|
49 |
|
| 50 |
|
50 |
|
| 51 |
For I :
|
51 |
For I :
|
| 52 |
|
52 |
|
| 53 |
EXPARG = Largest working precision argument that the library
|
53 |
EXPARG = Largest working precision argument that the library
|
| 54 |
EXP routine can handle and upper limit on the
|
54 |
EXP routine can handle and upper limit on the
|
| 55 |
magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)
|
55 |
magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)
|
| 56 |
|
56 |
|
| 57 |
For I and J :
|
57 |
For I and J :
|
| 58 |
|
58 |
|
| 59 |
xlrg_IJ = xlrg_BESS_IJ (was = XLARGE). Upper limit on the magnitude of X
|
59 |
xlrg_BESS_IJ (was = XLARGE) = xlrg_IJ in table below.
|
| - |
|
60 |
Upper limit on the magnitude of X (when IZE=2 for I()).
|
| 60 |
(when IZE=2 for I()). Bear in mind that if floor(abs(x)) =: N, then
|
61 |
Bear in mind that if floor(abs(x)) =: N, then
|
| 61 |
at least N iterations of the backward recursion will be executed.
|
62 |
at least N iterations of the backward recursion will be executed.
|
| 62 |
The value of 10 ^ 4 was used till Feb.2009, when it was increased
|
63 |
The value of 10 ^ 4 was used till Feb.2009, when it was increased
|
| 63 |
to 10 ^ 5 (= 1e5).
|
64 |
to 10 ^ 5 (= 1e5).
|
| 64 |
|
65 |
|
| 65 |
For j :
|
66 |
For J :
|
| - |
|
67 |
|
| 66 |
XMIN_J = Smallest acceptable argument for RBESY; approximately
|
68 |
XMIN_J = Smallest acceptable argument for RBESY; approximately
|
| 67 |
max(2*beta ^ minexp, 2/XINF), rounded up
|
69 |
max(2*beta ^ minexp, 2/XINF), rounded up
|
| 68 |
|
70 |
|
| 69 |
For Y :
|
71 |
For Y :
|
| 70 |
|
72 |
|
| 71 |
xlrg_Y = (was = XLARGE). Upper bound on X;
|
73 |
xlrg_Y = (was = XLARGE). Upper bound on X;
|
| 72 |
approximately 1/DEL, because the sine and cosine functions
|
74 |
approximately 1/DEL, because the sine and cosine functions
|
| 73 |
have lost about half of their precision at that point.
|
75 |
have lost about half of their precision at that point.
|
| 74 |
|
76 |
|
| 75 |
EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS).
|
77 |
EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS).
|
| 76 |
THRESH = Lower bound for use of the asymptotic form;
|
78 |
THRESH = Lower bound for use of the asymptotic form;
|
| 77 |
approximately AINT(-LOG10(EPS/2.0))+1.0
|
79 |
approximately AINT(-LOG10(EPS/2.0))+1.0
|
| 78 |
|
80 |
|
| 79 |
|
81 |
|
| 80 |
For K :
|
82 |
For K :
|
| 81 |
|
83 |
|
| 82 |
xmax_k = (was = XMAX). Upper limit on the magnitude of X when ize = 1;
|
84 |
xmax_k = (was = XMAX). Upper limit on the magnitude of X when ize = 1;
|
| 83 |
i.e. maximal x for UNscaled answer.
|
85 |
i.e. maximal x for UNscaled answer.
|
| 84 |
|
86 |
|
| 85 |
Solution to equation:
|
87 |
Solution to equation:
|
| Line 134... |
Line 136... |
| 134 |
|
136 |
|
| 135 |
/* sqrt(DBL_MIN) = 1.491668e-154 */
|
137 |
/* sqrt(DBL_MIN) = 1.491668e-154 */
|
| 136 |
#define sqxmin_BESS_K 1.49e-154
|
138 |
#define sqxmin_BESS_K 1.49e-154
|
| 137 |
|
139 |
|
| 138 |
/* x < eps_sinc <==> sin(x)/x == 1 (particularly "==>");
|
140 |
/* x < eps_sinc <==> sin(x)/x == 1 (particularly "==>");
|
| 139 |
Linux (around 2001-02) gives 2.14946906753213e-08
|
141 |
Linux (around 2001-02) gives 2.14946906753213e-08; 2026-01: 2.14911933289082e-8 ("=" Solaris 2.5.1 below !)
|
| 140 |
Solaris 2.5.1 gives 2.14911933289084e-08
|
142 |
Solaris 2.5.1 gives 2.14911933289084e-08
|
| 141 |
*/
|
143 |
*/
|
| 142 |
#define M_eps_sinc 2.149e-8
|
144 |
#define M_eps_sinc 2.149119e-8
|
| - |
|
145 |
|
The R Project SVN
R