/* *******************************************************************

 Explanation of machine-dependent constants

   beta   = Radix for the floating-point system
   minexp = Smallest representable power of beta
   maxexp = Smallest power of beta that overflows
   it = p = Number of bits (base-beta digits) in the mantissa
            (significand) of a working precision (floating-point) variable
   NSIG   = Decimal significance desired.  Should be set to
            INT(LOG10(2)*it+1).  Setting NSIG lower will result
            in decreased accuracy while setting NSIG higher will
            increase CPU time without increasing accuracy.  The
            truncation error is limited to a relative error of
            T=.5*10^(-NSIG).
   ENTEN  = 10 ^ K, where K is the largest long such that
            ENTEN is machine-representable in working precision
   ENSIG  = 10 ^ NSIG
   RTNSIG = 10 ^ (-K) for the smallest long K such that
            K >= NSIG/4
   ENMTEN = Smallest ABS(X) such that X/4 does not underflow
   XINF   = Largest positive machine number; approximately beta ^ maxexp
            == DBL_MAX (defined in  #include <float.h>)
   SQXMIN = Square root of beta ^ minexp = sqrt(DBL_MIN)

   EPS    = The smallest positive floating-point number such that 1.0+EPS > 1.0
          = beta ^ (-p)  == DBL_EPSILON


  For I :

   EXPARG = Largest working precision argument that the library
            EXP routine can handle and upper limit on the
            magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)

  For I and J :

   xlrg_IJ = xlrg_BESS_IJ (was = XLARGE). Upper limit on the magnitude of X
            (when IZE=2 for I()).  Bear in mind that if floor(abs(x)) =: N, then
            at least N iterations of the backward recursion will be executed.
            The value of 10 ^ 4 was used till Feb.2009, when it was increased
            to 10 ^ 5 (= 1e5).

  For j :
   XMIN_J  = Smallest acceptable argument for RBESY; approximately
            max(2*beta ^ minexp, 2/XINF), rounded up

  For Y :

   xlrg_Y =  (was = XLARGE). Upper bound on X;
            approximately 1/DEL, because the sine and cosine functions
            have lost about half of their precision at that point.

   EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS).
   THRESH = Lower bound for use of the asymptotic form;
            approximately AINT(-LOG10(EPS/2.0))+1.0


  For K :

   xmax_k =  (was = XMAX). Upper limit on the magnitude of X when ize = 1;
            i.e. maximal x for UNscaled answer.

            Solution to equation:
               W(X) * (1 -1/8 X + 9/128 X^2) = beta ^ minexp
            where  W(X) = EXP(-X)*SQRT(PI/2X)

 --------------------------------------------------------------------

     Approximate values for some important machines are:

                  beta minexp maxexp it NSIG ENTEN ENSIG RTNSIG ENMTEN   EXPARG
 IEEE (IBM/XT,
   SUN, etc.) (S.P.)  2   -126  128  24   8  1e38   1e8   1e-2  4.70e-38     88
 IEEE   (...) (D.P.)  2  -1022 1024  53  16  1e308  1e16  1e-4  8.90e-308   709
 CRAY-1       (S.P.)  2  -8193 8191  48  15  1e2465 1e15  1e-4  1.84e-2466 5677
 Cyber 180/855
   under NOS  (S.P.)  2   -975 1070  48  15  1e322  1e15  1e-4  1.25e-293   741
 IBM 3033     (D.P.) 16    -65   63  14   5  1e75   1e5   1e-2  2.16e-78    174
 VAX          (S.P.)  2   -128  127  24   8  1e38   1e8   1e-2  1.17e-38     88
 VAX D-Format (D.P.)  2   -128  127  56  17  1e38   1e17  1e-5  1.17e-38     88
 VAX G-Format (D.P.)  2  -1024 1023  53  16  1e307  1e16  1e-4  2.22e-308   709


And routine specific :

                    xlrg_IJ xlrg_Y xmax_k EPS_SINC XMIN_J    XINF   THRESH
 IEEE (IBM/XT,
   SUN, etc.) (S.P.)    1e4  1e4   85.337  1e-4  2.36e-38   3.40e38     8.
 IEEE   (...) (D.P.)    1e4  1e8  705.342  1e-8  4.46e-308  1.79e308   16.
 CRAY-1       (S.P.)    1e4  2e7 5674.858  5e-8  3.67e-2466 5.45e2465  15.
 Cyber 180/855
   under NOS  (S.P.)    1e4  2e7  672.788  5e-8  6.28e-294  1.26e322   15.
 IBM 3033     (D.P.)    1e4  1e8  177.852  1e-8  2.77e-76   7.23e75    17.
 VAX          (S.P.)    1e4  1e4   86.715  1e-4  1.18e-38   1.70e38     8.
 VAX e-Format (D.P.)    1e4  1e9   86.715  1e-9  1.18e-38   1.70e38    17.
 VAX G-Format (D.P.)    1e4  1e8  706.728  1e-8  2.23e-308  8.98e307   16.


   beta - radix for the floating-point representation
   maxexp - the smallest positive power of beta that overflows
   XBIG - the largest argument for which GAMMA(X) is representable
        in the machine, i.e., the solution to the equation
        GAMMA(XBIG) = beta**maxexp
   XINF - the largest machine representable floating-point number;
        approximately beta**maxexp
   EPS  - the smallest positive floating-point number such that  1.0+EPS > 1.0
   XMININ - the smallest positive floating-point number such that
        1/XMININ is machine representable

   Approximate values for some important machines are:

   beta       maxexp         XBIG

   CRAY-1               (S.P.)        2         8191        966.961
   Cyber 180/855
   under NOS    (S.P.)        2         1070        177.803
   IEEE (IBM/XT,
   SUN, etc.)   (S.P.)        2          128        35.040
   IEEE (IBM/XT,
   SUN, etc.)   (D.P.)        2         1024        171.624
   IBM 3033     (D.P.)       16           63        57.574
   VAX D-Format (D.P.)        2          127        34.844
   VAX G-Format (D.P.)        2         1023        171.489

   XINF  EPS        XMININ

   CRAY-1               (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
   Cyber 180/855
   under NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
   IEEE (IBM/XT,
   SUN, etc.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
   IEEE (IBM/XT,
   SUN, etc.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
   IBM 3033     (D.P.)   7.23D+75     2.22D-16    1.39D-76
   VAX D-Format (D.P.)   1.70D+38     1.39D-17    5.88D-39
   VAX G-Format (D.P.)   8.98D+307    1.11D-16    1.12D-308


