TMB Documentation  v1.9.11
bessel.h
1 /*
2  * R : A Computer Language for Statistical Data Analysis
3  * Copyright (C) 2001-2014 R Core Team
4  *
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
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, a copy is available at
17  * https://www.R-project.org/Licenses/
18  */
19 
20 /* Constants und Documentation that apply to several of the
21  * ./bessel_[ijky].c files */
22 
23 /* *******************************************************************
24 
25  Explanation of machine-dependent constants
26 
27  beta = Radix for the floating-point system
28  minexp = Smallest representable power of beta
29  maxexp = Smallest power of beta that overflows
30  it = p = Number of bits (base-beta digits) in the mantissa
31  (significand) of a working precision (floating-point) variable
32  NSIG = Decimal significance desired. Should be set to
33  INT(LOG10(2)*it+1). Setting NSIG lower will result
34  in decreased accuracy while setting NSIG higher will
35  increase CPU time without increasing accuracy. The
36  truncation error is limited to a relative error of
37  T=.5*10^(-NSIG).
38  ENTEN = 10 ^ K, where K is the largest int such that
39  ENTEN is machine-representable in working precision
40  ENSIG = 10 ^ NSIG
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
43  XINF = Largest positive machine number; approximately beta ^ maxexp
44  == DBL_MAX (defined in #include <float.h>)
45  SQXMIN = Square root of beta ^ minexp = sqrt(DBL_MIN)
46 
47  EPS = The smallest positive floating-point number such that 1.0+EPS > 1.0
48  = beta ^ (-p) == DBL_EPSILON
49 
50 
51  For I :
52 
53  EXPARG = Largest working precision argument that the library
54  EXP routine can handle and upper limit on the
55  magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)
56 
57  For I and J :
58 
59  xlrg_IJ = xlrg_BESS_IJ (was = XLARGE). Upper limit on the magnitude of X
60  (when IZE=2 for I()). Bear in mind that if floor(abs(x)) =: N, then
61  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  to 10 ^ 5 (= 1e5).
64 
65  For j :
66  XMIN_J = Smallest acceptable argument for RBESY; approximately
67  max(2*beta ^ minexp, 2/XINF), rounded up
68 
69  For Y :
70 
71  xlrg_Y = (was = XLARGE). Upper bound on X;
72  approximately 1/DEL, because the sine and cosine functions
73  have lost about half of their precision at that point.
74 
75  EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS).
76  THRESH = Lower bound for use of the asymptotic form;
77  approximately AINT(-LOG10(EPS/2.0))+1.0
78 
79 
80  For K :
81 
82  xmax_k = (was = XMAX). Upper limit on the magnitude of X when ize = 1;
83  i.e. maximal x for UNscaled answer.
84 
85  Solution to equation:
86  W(X) * (1 -1/8 X + 9/128 X^2) = beta ^ minexp
87  where W(X) = EXP(-X)*SQRT(PI/2X)
88 
89  --------------------------------------------------------------------
90 
91  Approximate values for some important machines are:
92 
93  beta minexp maxexp it NSIG ENTEN ENSIG RTNSIG ENMTEN EXPARG
94  IEEE (IBM/XT,
95  SUN, etc.) (S.P.) 2 -126 128 24 8 1e38 1e8 1e-2 4.70e-38 88
96  IEEE (...) (D.P.) 2 -1022 1024 53 16 1e308 1e16 1e-4 8.90e-308 709
97  CRAY-1 (S.P.) 2 -8193 8191 48 15 1e2465 1e15 1e-4 1.84e-2466 5677
98  Cyber 180/855
99  under NOS (S.P.) 2 -975 1070 48 15 1e322 1e15 1e-4 1.25e-293 741
100  IBM 3033 (D.P.) 16 -65 63 14 5 1e75 1e5 1e-2 2.16e-78 174
101  VAX (S.P.) 2 -128 127 24 8 1e38 1e8 1e-2 1.17e-38 88
102  VAX D-Format (D.P.) 2 -128 127 56 17 1e38 1e17 1e-5 1.17e-38 88
103  VAX G-Format (D.P.) 2 -1024 1023 53 16 1e307 1e16 1e-4 2.22e-308 709
104 
105 
106 And routine specific :
107 
108  xlrg_IJ xlrg_Y xmax_k EPS_SINC XMIN_J XINF THRESH
109  IEEE (IBM/XT,
110  SUN, etc.) (S.P.) 1e4 1e4 85.337 1e-4 2.36e-38 3.40e38 8.
111  IEEE (...) (D.P.) 1e4 1e8 705.342 1e-8 4.46e-308 1.79e308 16.
112  CRAY-1 (S.P.) 1e4 2e7 5674.858 5e-8 3.67e-2466 5.45e2465 15.
113  Cyber 180/855
114  under NOS (S.P.) 1e4 2e7 672.788 5e-8 6.28e-294 1.26e322 15.
115  IBM 3033 (D.P.) 1e4 1e8 177.852 1e-8 2.77e-76 7.23e75 17.
116  VAX (S.P.) 1e4 1e4 86.715 1e-4 1.18e-38 1.70e38 8.
117  VAX e-Format (D.P.) 1e4 1e9 86.715 1e-9 1.18e-38 1.70e38 17.
118  VAX G-Format (D.P.) 1e4 1e8 706.728 1e-8 2.23e-308 8.98e307 16.
119 
120 */
121 #define nsig_BESS 16
122 #define ensig_BESS 1e16
123 #define rtnsig_BESS 1e-4
124 #define enmten_BESS 8.9e-308
125 #define enten_BESS 1e308
126 
127 #define exparg_BESS 709.
128 #define xlrg_BESS_IJ 1e5
129 #define xlrg_BESS_Y 1e8
130 #define thresh_BESS_Y 16.
131 
132 #define xmax_BESS_K 705.342/* maximal x for UNscaled answer */
133 
134 
135 /* sqrt(DBL_MIN) = 1.491668e-154 */
136 #define sqxmin_BESS_K 1.49e-154
137 
138 /* x < eps_sinc <==> sin(x)/x == 1 (particularly "==>");
139  Linux (around 2001-02) gives 2.14946906753213e-08
140  Solaris 2.5.1 gives 2.14911933289084e-08
141 */
142 #define M_eps_sinc 2.149e-8
License: GPL v2