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 143