src/libsphinxbase/util/slamch.c

   1 /* src/slamch.f -- translated by f2c (version 20050501).
   2    You must link the resulting object file with libf2c:
   3         on Microsoft Windows system, link with libf2c.lib;
   4         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
   5         or, if you install libf2c.a in a standard place, with -lf2c -lm
   6         -- in that order, at the end of the command line, as in
   7                 cc *.o -lf2c -lm
   8         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
   9
  10                 http://www.netlib.org/f2c/libf2c.zip
  11 */
  12
  13 #include "sphinxbase/f2c.h"
  14
  15 #ifdef _MSC_VER
  16 #pragma warning (disable: 4244)
  17 #endif
  18
  19 /* Table of constant values */
  20
  21 static integer c__1 = 1;
  22 static real c_b32 = 0.f;
  23
  24 doublereal
  25 slamch_(char *cmach, ftnlen cmach_len)
  26 {
  27     /* Initialized data */
  28
  29     static logical first = TRUE_;
  30
  31     /* System generated locals */
  32     integer i__1;
  33     real ret_val;
  34
  35     /* Builtin functions */
  36     double pow_ri(real *, integer *);
  37
  38     /* Local variables */
  39     static real t;
  40     static integer it;
  41     static real rnd, eps, base;
  42     static integer beta;
  43     static real emin, prec, emax;
  44     static integer imin, imax;
  45     static logical lrnd;
  46     static real rmin, rmax, rmach;
  47     extern logical lsame_(char *, char *, ftnlen, ftnlen);
  48     static real small, sfmin;
  49     extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real
  50                                         *, integer *, real *, integer *,
  51                                         real *);
  52
  53
  54 /*  -- LAPACK auxiliary routine (version 3.0) -- */
  55 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
  56 /*     Courant Institute, Argonne National Lab, and Rice University */
  57 /*     October 31, 1992 */
  58
  59 /*     .. Scalar Arguments .. */
  60 /*     .. */
  61
  62 /*  Purpose */
  63 /*  ======= */
  64
  65 /*  SLAMCH determines single precision machine parameters. */
  66
  67 /*  Arguments */
  68 /*  ========= */
  69
  70 /*  CMACH   (input) CHARACTER*1 */
  71 /*          Specifies the value to be returned by SLAMCH: */
  72 /*          = 'E' or 'e',   SLAMCH := eps */
  73 /*          = 'S' or 's ,   SLAMCH := sfmin */
  74 /*          = 'B' or 'b',   SLAMCH := base */
  75 /*          = 'P' or 'p',   SLAMCH := eps*base */
  76 /*          = 'N' or 'n',   SLAMCH := t */
  77 /*          = 'R' or 'r',   SLAMCH := rnd */
  78 /*          = 'M' or 'm',   SLAMCH := emin */
  79 /*          = 'U' or 'u',   SLAMCH := rmin */
  80 /*          = 'L' or 'l',   SLAMCH := emax */
  81 /*          = 'O' or 'o',   SLAMCH := rmax */
  82
  83 /*          where */
  84
  85 /*          eps   = relative machine precision */
  86 /*          sfmin = safe minimum, such that 1/sfmin does not overflow */
  87 /*          base  = base of the machine */
  88 /*          prec  = eps*base */
  89 /*          t     = number of (base) digits in the mantissa */
  90 /*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise */
  91 /*          emin  = minimum exponent before (gradual) underflow */
  92 /*          rmin  = underflow threshold - base**(emin-1) */
  93 /*          emax  = largest exponent before overflow */
  94 /*          rmax  = overflow threshold  - (base**emax)*(1-eps) */
  95
  96 /* ===================================================================== */
  97
  98 /*     .. Parameters .. */
  99 /*     .. */
 100 /*     .. Local Scalars .. */
 101 /*     .. */
 102 /*     .. External Functions .. */
 103 /*     .. */
 104 /*     .. External Subroutines .. */
 105 /*     .. */
 106 /*     .. Save statement .. */
 107 /*     .. */
 108 /*     .. Data statements .. */
 109 /*     .. */
 110 /*     .. Executable Statements .. */
 111
 112     if (first) {
 113         first = FALSE_;
 114         slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
 115         base = (real) beta;
 116         t = (real) it;
 117         if (lrnd) {
 118             rnd = 1.f;
 119             i__1 = 1 - it;
 120             eps = pow_ri(&base, &i__1) / 2;
 121         }
 122         else {
 123             rnd = 0.f;
 124             i__1 = 1 - it;
 125             eps = pow_ri(&base, &i__1);
 126         }
 127         prec = eps * base;
 128         emin = (real) imin;
 129         emax = (real) imax;
 130         sfmin = rmin;
 131         small = 1.f / rmax;
 132         if (small >= sfmin) {
 133
 134 /*           Use SMALL plus a bit, to avoid the possibility of rounding */
 135 /*           causing overflow when computing  1/sfmin. */
 136
 137             sfmin = small * (eps + 1.f);
 138         }
 139     }
 140
 141     if (lsame_(cmach, "E", (ftnlen) 1, (ftnlen) 1)) {
 142         rmach = eps;
 143     }
 144     else if (lsame_(cmach, "S", (ftnlen) 1, (ftnlen) 1)) {
 145         rmach = sfmin;
 146     }
 147     else if (lsame_(cmach, "B", (ftnlen) 1, (ftnlen) 1)) {
 148         rmach = base;
 149     }
 150     else if (lsame_(cmach, "P", (ftnlen) 1, (ftnlen) 1)) {
 151         rmach = prec;
 152     }
 153     else if (lsame_(cmach, "N", (ftnlen) 1, (ftnlen) 1)) {
 154         rmach = t;
 155     }
 156     else if (lsame_(cmach, "R", (ftnlen) 1, (ftnlen) 1)) {
 157         rmach = rnd;
 158     }
 159     else if (lsame_(cmach, "M", (ftnlen) 1, (ftnlen) 1)) {
 160         rmach = emin;
 161     }
 162     else if (lsame_(cmach, "U", (ftnlen) 1, (ftnlen) 1)) {
 163         rmach = rmin;
 164     }
 165     else if (lsame_(cmach, "L", (ftnlen) 1, (ftnlen) 1)) {
 166         rmach = emax;
 167     }
 168     else if (lsame_(cmach, "O", (ftnlen) 1, (ftnlen) 1)) {
 169         rmach = rmax;
 170     }
 171
 172     ret_val = rmach;
 173     return ret_val;
 174
 175 /*     End of SLAMCH */
 176
 177 }                               /* slamch_ */
 178
 179
 180 /* *********************************************************************** */
 181
 182 /* Subroutine */ int
 183 slamc1_(integer * beta, integer * t, logical * rnd, logical * ieee1)
 184 {
 185     /* Initialized data */
 186
 187     static logical first = TRUE_;
 188
 189     /* System generated locals */
 190     real r__1, r__2;
 191
 192     /* Local variables */
 193     static real a, b, c__, f, t1, t2;
 194     static integer lt;
 195     static real one, qtr;
 196     static logical lrnd;
 197     static integer lbeta;
 198     static real savec;
 199     static logical lieee1;
 200     extern doublereal slamc3_(real *, real *);
 201
 202
 203 /*  -- LAPACK auxiliary routine (version 3.0) -- */
 204 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
 205 /*     Courant Institute, Argonne National Lab, and Rice University */
 206 /*     October 31, 1992 */
 207
 208 /*     .. Scalar Arguments .. */
 209 /*     .. */
 210
 211 /*  Purpose */
 212 /*  ======= */
 213
 214 /*  SLAMC1 determines the machine parameters given by BETA, T, RND, and */
 215 /*  IEEE1. */
 216
 217 /*  Arguments */
 218 /*  ========= */
 219
 220 /*  BETA    (output) INTEGER */
 221 /*          The base of the machine. */
 222
 223 /*  T       (output) INTEGER */
 224 /*          The number of ( BETA ) digits in the mantissa. */
 225
 226 /*  RND     (output) LOGICAL */
 227 /*          Specifies whether proper rounding  ( RND = .TRUE. )  or */
 228 /*          chopping  ( RND = .FALSE. )  occurs in addition. This may not */
 229 /*          be a reliable guide to the way in which the machine performs */
 230 /*          its arithmetic. */
 231
 232 /*  IEEE1   (output) LOGICAL */
 233 /*          Specifies whether rounding appears to be done in the IEEE */
 234 /*          'round to nearest' style. */
 235
 236 /*  Further Details */
 237 /*  =============== */
 238
 239 /*  The routine is based on the routine  ENVRON  by Malcolm and */
 240 /*  incorporates suggestions by Gentleman and Marovich. See */
 241
 242 /*     Malcolm M. A. (1972) Algorithms to reveal properties of */
 243 /*        floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
 244
 245 /*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
 246 /*        that reveal properties of floating point arithmetic units. */
 247 /*        Comms. of the ACM, 17, 276-277. */
 248
 249 /* ===================================================================== */
 250
 251 /*     .. Local Scalars .. */
 252 /*     .. */
 253 /*     .. External Functions .. */
 254 /*     .. */
 255 /*     .. Save statement .. */
 256 /*     .. */
 257 /*     .. Data statements .. */
 258 /*     .. */
 259 /*     .. Executable Statements .. */
 260
 261     if (first) {
 262         first = FALSE_;
 263         one = 1.f;
 264
 265 /*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA, */
 266 /*        IEEE1, T and RND. */
 267
 268 /*        Throughout this routine  we use the function  SLAMC3  to ensure */
 269 /*        that relevant values are  stored and not held in registers,  or */
 270 /*        are not affected by optimizers. */
 271
 272 /*        Compute  a = 2.0**m  with the  smallest positive integer m such */
 273 /*        that */
 274
 275 /*           fl( a + 1.0 ) = a. */
 276
 277         a = 1.f;
 278         c__ = 1.f;
 279
 280 /* +       WHILE( C.EQ.ONE )LOOP */
 281       L10:
 282         if (c__ == one) {
 283             a *= 2;
 284             c__ = slamc3_(&a, &one);
 285             r__1 = -a;
 286             c__ = slamc3_(&c__, &r__1);
 287             goto L10;
 288         }
 289 /* +       END WHILE */
 290
 291 /*        Now compute  b = 2.0**m  with the smallest positive integer m */
 292 /*        such that */
 293
 294 /*           fl( a + b ) .gt. a. */
 295
 296         b = 1.f;
 297         c__ = slamc3_(&a, &b);
 298
 299 /* +       WHILE( C.EQ.A )LOOP */
 300       L20:
 301         if (c__ == a) {
 302             b *= 2;
 303             c__ = slamc3_(&a, &b);
 304             goto L20;
 305         }
 306 /* +       END WHILE */
 307
 308 /*        Now compute the base.  a and c  are neighbouring floating point */
 309 /*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so */
 310 /*        their difference is beta. Adding 0.25 to c is to ensure that it */
 311 /*        is truncated to beta and not ( beta - 1 ). */
 312
 313         qtr = one / 4;
 314         savec = c__;
 315         r__1 = -a;
 316         c__ = slamc3_(&c__, &r__1);
 317         lbeta = c__ + qtr;
 318
 319 /*        Now determine whether rounding or chopping occurs,  by adding a */
 320 /*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a. */
 321
 322         b = (real) lbeta;
 323         r__1 = b / 2;
 324         r__2 = -b / 100;
 325         f = slamc3_(&r__1, &r__2);
 326         c__ = slamc3_(&f, &a);
 327         if (c__ == a) {
 328             lrnd = TRUE_;
 329         }
 330         else {
 331             lrnd = FALSE_;
 332         }
 333         r__1 = b / 2;
 334         r__2 = b / 100;
 335         f = slamc3_(&r__1, &r__2);
 336         c__ = slamc3_(&f, &a);
 337         if (lrnd && c__ == a) {
 338             lrnd = FALSE_;
 339         }
 340
 341 /*        Try and decide whether rounding is done in the  IEEE  'round to */
 342 /*        nearest' style. B/2 is half a unit in the last place of the two */
 343 /*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit */
 344 /*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change */
 345 /*        A, but adding B/2 to SAVEC should change SAVEC. */
 346
 347         r__1 = b / 2;
 348         t1 = slamc3_(&r__1, &a);
 349         r__1 = b / 2;
 350         t2 = slamc3_(&r__1, &savec);
 351         lieee1 = t1 == a && t2 > savec && lrnd;
 352
 353 /*        Now find  the  mantissa, t.  It should  be the  integer part of */
 354 /*        log to the base beta of a,  however it is safer to determine  t */
 355 /*        by powering.  So we find t as the smallest positive integer for */
 356 /*        which */
 357
 358 /*           fl( beta**t + 1.0 ) = 1.0. */
 359
 360         lt = 0;
 361         a = 1.f;
 362         c__ = 1.f;
 363
 364 /* +       WHILE( C.EQ.ONE )LOOP */
 365       L30:
 366         if (c__ == one) {
 367             ++lt;
 368             a *= lbeta;
 369             c__ = slamc3_(&a, &one);
 370             r__1 = -a;
 371             c__ = slamc3_(&c__, &r__1);
 372             goto L30;
 373         }
 374 /* +       END WHILE */
 375
 376     }
 377
 378     *beta = lbeta;
 379     *t = lt;
 380     *rnd = lrnd;
 381     *ieee1 = lieee1;
 382     return 0;
 383
 384 /*     End of SLAMC1 */
 385
 386 }                               /* slamc1_ */
 387
 388
 389 /* *********************************************************************** */
 390
 391 /* Subroutine */ int
 392 slamc2_(integer * beta, integer * t, logical * rnd, real *
 393         eps, integer * emin, real * rmin, integer * emax, real * rmax)
 394 {
 395     /* Initialized data */
 396
 397     static logical first = TRUE_;
 398     static logical iwarn = FALSE_;
 399
 400     /* Format strings */
 401     static char fmt_9999[] =
 402         "(//\002 WARNING. The value EMIN may be incorre"
 403         "ct:-\002,\002  EMIN = \002,i8,/\002 If, after inspection, the va"
 404         "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
 405         " the IF block as marked within the code of routine\002,\002 SLAM"
 406         "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
 407
 408     /* System generated locals */
 409     integer i__1;
 410     real r__1, r__2, r__3, r__4, r__5;
 411
 412     /* Builtin functions */
 413     double pow_ri(real *, integer *);
 414     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen),
 415         e_wsfe(void);
 416
 417     /* Local variables */
 418     static real a, b, c__;
 419     static integer i__, lt;
 420     static real one, two;
 421     static logical ieee;
 422     static real half;
 423     static logical lrnd;
 424     static real leps, zero;
 425     static integer lbeta;
 426     static real rbase;
 427     static integer lemin, lemax, gnmin;
 428     static real small;
 429     static integer gpmin;
 430     static real third, lrmin, lrmax, sixth;
 431     static logical lieee1;
 432     extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
 433                                         logical *);
 434     extern doublereal slamc3_(real *, real *);
 435     extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
 436         slamc5_(integer *, integer *, integer *, logical *, integer *,
 437                 real *);
 438     static integer ngnmin, ngpmin;
 439
 440     /* Fortran I/O blocks */
 441     static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
 442
 443
 444
 445 /*  -- LAPACK auxiliary routine (version 3.0) -- */
 446 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
 447 /*     Courant Institute, Argonne National Lab, and Rice University */
 448 /*     October 31, 1992 */
 449
 450 /*     .. Scalar Arguments .. */
 451 /*     .. */
 452
 453 /*  Purpose */
 454 /*  ======= */
 455
 456 /*  SLAMC2 determines the machine parameters specified in its argument */
 457 /*  list. */
 458
 459 /*  Arguments */
 460 /*  ========= */
 461
 462 /*  BETA    (output) INTEGER */
 463 /*          The base of the machine. */
 464
 465 /*  T       (output) INTEGER */
 466 /*          The number of ( BETA ) digits in the mantissa. */
 467
 468 /*  RND     (output) LOGICAL */
 469 /*          Specifies whether proper rounding  ( RND = .TRUE. )  or */
 470 /*          chopping  ( RND = .FALSE. )  occurs in addition. This may not */
 471 /*          be a reliable guide to the way in which the machine performs */
 472 /*          its arithmetic. */
 473
 474 /*  EPS     (output) REAL */
 475 /*          The smallest positive number such that */
 476
 477 /*             fl( 1.0 - EPS ) .LT. 1.0, */
 478
 479 /*          where fl denotes the computed value. */
 480
 481 /*  EMIN    (output) INTEGER */
 482 /*          The minimum exponent before (gradual) underflow occurs. */
 483
 484 /*  RMIN    (output) REAL */
 485 /*          The smallest normalized number for the machine, given by */
 486 /*          BASE**( EMIN - 1 ), where  BASE  is the floating point value */
 487 /*          of BETA. */
 488
 489 /*  EMAX    (output) INTEGER */
 490 /*          The maximum exponent before overflow occurs. */
 491
 492 /*  RMAX    (output) REAL */
 493 /*          The largest positive number for the machine, given by */
 494 /*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point */
 495 /*          value of BETA. */
 496
 497 /*  Further Details */
 498 /*  =============== */
 499
 500 /*  The computation of  EPS  is based on a routine PARANOIA by */
 501 /*  W. Kahan of the University of California at Berkeley. */
 502
 503 /* ===================================================================== */
 504
 505 /*     .. Local Scalars .. */
 506 /*     .. */
 507 /*     .. External Functions .. */
 508 /*     .. */
 509 /*     .. External Subroutines .. */
 510 /*     .. */
 511 /*     .. Intrinsic Functions .. */
 512 /*     .. */
 513 /*     .. Save statement .. */
 514 /*     .. */
 515 /*     .. Data statements .. */
 516 /*     .. */
 517 /*     .. Executable Statements .. */
 518
 519     if (first) {
 520         first = FALSE_;
 521         zero = 0.f;
 522         one = 1.f;
 523         two = 2.f;
 524
 525 /*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of */
 526 /*        BETA, T, RND, EPS, EMIN and RMIN. */
 527
 528 /*        Throughout this routine  we use the function  SLAMC3  to ensure */
 529 /*        that relevant values are stored  and not held in registers,  or */
 530 /*        are not affected by optimizers. */
 531
 532 /*        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1. */
 533
 534         slamc1_(&lbeta, &lt, &lrnd, &lieee1);
 535
 536 /*        Start to find EPS. */
 537
 538         b = (real) lbeta;
 539         i__1 = -lt;
 540         a = pow_ri(&b, &i__1);
 541         leps = a;
 542
 543 /*        Try some tricks to see whether or not this is the correct  EPS. */
 544
 545         b = two / 3;
 546         half = one / 2;
 547         r__1 = -half;
 548         sixth = slamc3_(&b, &r__1);
 549         third = slamc3_(&sixth, &sixth);
 550         r__1 = -half;
 551         b = slamc3_(&third, &r__1);
 552         b = slamc3_(&b, &sixth);
 553         b = dabs(b);
 554         if (b < leps) {
 555             b = leps;
 556         }
 557
 558         leps = 1.f;
 559
 560 /* +       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
 561       L10:
 562         if (leps > b && b > zero) {
 563             leps = b;
 564             r__1 = half * leps;
 565 /* Computing 5th power */
 566             r__3 = two, r__4 = r__3, r__3 *= r__3;
 567 /* Computing 2nd power */
 568             r__5 = leps;
 569             r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
 570             c__ = slamc3_(&r__1, &r__2);
 571             r__1 = -c__;
 572             c__ = slamc3_(&half, &r__1);
 573             b = slamc3_(&half, &c__);
 574             r__1 = -b;
 575             c__ = slamc3_(&half, &r__1);
 576             b = slamc3_(&half, &c__);
 577             goto L10;
 578         }
 579 /* +       END WHILE */
 580
 581         if (a < leps) {
 582             leps = a;
 583         }
 584
 585 /*        Computation of EPS complete. */
 586
 587 /*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)). */
 588 /*        Keep dividing  A by BETA until (gradual) underflow occurs. This */
 589 /*        is detected when we cannot recover the previous A. */
 590
 591         rbase = one / lbeta;
 592         small = one;
 593         for (i__ = 1; i__ <= 3; ++i__) {
 594             r__1 = small * rbase;
 595             small = slamc3_(&r__1, &zero);
 596 /* L20: */
 597         }
 598         a = slamc3_(&one, &small);
 599         slamc4_(&ngpmin, &one, &lbeta);
 600         r__1 = -one;
 601         slamc4_(&ngnmin, &r__1, &lbeta);
 602         slamc4_(&gpmin, &a, &lbeta);
 603         r__1 = -a;
 604         slamc4_(&gnmin, &r__1, &lbeta);
 605         ieee = FALSE_;
 606
 607         if (ngpmin == ngnmin && gpmin == gnmin) {
 608             if (ngpmin == gpmin) {
 609                 lemin = ngpmin;
 610 /*            ( Non twos-complement machines, no gradual underflow; */
 611 /*              e.g.,  VAX ) */
 612             }
 613             else if (gpmin - ngpmin == 3) {
 614                 lemin = ngpmin - 1 + lt;
 615                 ieee = TRUE_;
 616 /*            ( Non twos-complement machines, with gradual underflow; */
 617 /*              e.g., IEEE standard followers ) */
 618             }
 619             else {
 620                 lemin = min(ngpmin, gpmin);
 621 /*            ( A guess; no known machine ) */
 622                 iwarn = TRUE_;
 623             }
 624
 625         }
 626         else if (ngpmin == gpmin && ngnmin == gnmin) {
 627             if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
 628                 lemin = max(ngpmin, ngnmin);
 629 /*            ( Twos-complement machines, no gradual underflow; */
 630 /*              e.g., CYBER 205 ) */
 631             }
 632             else {
 633                 lemin = min(ngpmin, ngnmin);
 634 /*            ( A guess; no known machine ) */
 635                 iwarn = TRUE_;
 636             }
 637
 638         }
 639         else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1
 640                  && gpmin == gnmin) {
 641             if (gpmin - min(ngpmin, ngnmin) == 3) {
 642                 lemin = max(ngpmin, ngnmin) - 1 + lt;
 643 /*            ( Twos-complement machines with gradual underflow; */
 644 /*              no known machine ) */
 645             }
 646             else {
 647                 lemin = min(ngpmin, ngnmin);
 648 /*            ( A guess; no known machine ) */
 649                 iwarn = TRUE_;
 650             }
 651
 652         }
 653         else {
 654 /* Computing MIN */
 655             i__1 = min(ngpmin, ngnmin), i__1 = min(i__1, gpmin);
 656             lemin = min(i__1, gnmin);
 657 /*         ( A guess; no known machine ) */
 658             iwarn = TRUE_;
 659         }
 660 /* ** */
 661 /* Comment out this if block if EMIN is ok */
 662         if (iwarn) {
 663             first = TRUE_;
 664             s_wsfe(&io___58);
 665             do_fio(&c__1, (char *) &lemin, (ftnlen) sizeof(integer));
 666             e_wsfe();
 667         }
 668 /* ** */
 669
 670 /*        Assume IEEE arithmetic if we found denormalised  numbers above, */
 671 /*        or if arithmetic seems to round in the  IEEE style,  determined */
 672 /*        in routine SLAMC1. A true IEEE machine should have both  things */
 673 /*        true; however, faulty machines may have one or the other. */
 674
 675         ieee = ieee || lieee1;
 676
 677 /*        Compute  RMIN by successive division by  BETA. We could compute */
 678 /*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during */
 679 /*        this computation. */
 680
 681         lrmin = 1.f;
 682         i__1 = 1 - lemin;
 683         for (i__ = 1; i__ <= i__1; ++i__) {
 684             r__1 = lrmin * rbase;
 685             lrmin = slamc3_(&r__1, &zero);
 686 /* L30: */
 687         }
 688
 689 /*        Finally, call SLAMC5 to compute EMAX and RMAX. */
 690
 691         slamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
 692     }
 693
 694     *beta = lbeta;
 695     *t = lt;
 696     *rnd = lrnd;
 697     *eps = leps;
 698     *emin = lemin;
 699     *rmin = lrmin;
 700     *emax = lemax;
 701     *rmax = lrmax;
 702
 703     return 0;
 704
 705
 706 /*     End of SLAMC2 */
 707
 708 }                               /* slamc2_ */
 709
 710
 711 /* *********************************************************************** */
 712
 713 doublereal
 714 slamc3_(real * a, real * b)
 715 {
 716     /* System generated locals */
 717     real ret_val;
 718
 719
 720 /*  -- LAPACK auxiliary routine (version 3.0) -- */
 721 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
 722 /*     Courant Institute, Argonne National Lab, and Rice University */
 723 /*     October 31, 1992 */
 724
 725 /*     .. Scalar Arguments .. */
 726 /*     .. */
 727
 728 /*  Purpose */
 729 /*  ======= */
 730
 731 /*  SLAMC3  is intended to force  A  and  B  to be stored prior to doing */
 732 /*  the addition of  A  and  B ,  for use in situations where optimizers */
 733 /*  might hold one of these in a register. */
 734
 735 /*  Arguments */
 736 /*  ========= */
 737
 738 /*  A, B    (input) REAL */
 739 /*          The values A and B. */
 740
 741 /* ===================================================================== */
 742
 743 /*     .. Executable Statements .. */
 744
 745     ret_val = *a + *b;
 746
 747     return ret_val;
 748
 749 /*     End of SLAMC3 */
 750
 751 }                               /* slamc3_ */
 752
 753
 754 /* *********************************************************************** */
 755
 756 /* Subroutine */ int
 757 slamc4_(integer * emin, real * start, integer * base)
 758 {
 759     /* System generated locals */
 760     integer i__1;
 761     real r__1;
 762
 763     /* Local variables */
 764     static real a;
 765     static integer i__;
 766     static real b1, b2, c1, c2, d1, d2, one, zero, rbase;
 767     extern doublereal slamc3_(real *, real *);
 768
 769
 770 /*  -- LAPACK auxiliary routine (version 3.0) -- */
 771 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
 772 /*     Courant Institute, Argonne National Lab, and Rice University */
 773 /*     October 31, 1992 */
 774
 775 /*     .. Scalar Arguments .. */
 776 /*     .. */
 777
 778 /*  Purpose */
 779 /*  ======= */
 780
 781 /*  SLAMC4 is a service routine for SLAMC2. */
 782
 783 /*  Arguments */
 784 /*  ========= */
 785
 786 /*  EMIN    (output) EMIN */
 787 /*          The minimum exponent before (gradual) underflow, computed by */
 788 /*          setting A = START and dividing by BASE until the previous A */
 789 /*          can not be recovered. */
 790
 791 /*  START   (input) REAL */
 792 /*          The starting point for determining EMIN. */
 793
 794 /*  BASE    (input) INTEGER */
 795 /*          The base of the machine. */
 796
 797 /* ===================================================================== */
 798
 799 /*     .. Local Scalars .. */
 800 /*     .. */
 801 /*     .. External Functions .. */
 802 /*     .. */
 803 /*     .. Executable Statements .. */
 804
 805     a = *start;
 806     one = 1.f;
 807     rbase = one / *base;
 808     zero = 0.f;
 809     *emin = 1;
 810     r__1 = a * rbase;
 811     b1 = slamc3_(&r__1, &zero);
 812     c1 = a;
 813     c2 = a;
 814     d1 = a;
 815     d2 = a;
 816 /* +    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
 817 /*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP */
 818   L10:
 819     if (c1 == a && c2 == a && d1 == a && d2 == a) {
 820         --(*emin);
 821         a = b1;
 822         r__1 = a / *base;
 823         b1 = slamc3_(&r__1, &zero);
 824         r__1 = b1 * *base;
 825         c1 = slamc3_(&r__1, &zero);
 826         d1 = zero;
 827         i__1 = *base;
 828         for (i__ = 1; i__ <= i__1; ++i__) {
 829             d1 += b1;
 830 /* L20: */
 831         }
 832         r__1 = a * rbase;
 833         b2 = slamc3_(&r__1, &zero);
 834         r__1 = b2 / rbase;
 835         c2 = slamc3_(&r__1, &zero);
 836         d2 = zero;
 837         i__1 = *base;
 838         for (i__ = 1; i__ <= i__1; ++i__) {
 839             d2 += b2;
 840 /* L30: */
 841         }
 842         goto L10;
 843     }
 844 /* +    END WHILE */
 845
 846     return 0;
 847
 848 /*     End of SLAMC4 */
 849
 850 }                               /* slamc4_ */
 851
 852
 853 /* *********************************************************************** */
 854
 855 /* Subroutine */ int
 856 slamc5_(integer * beta, integer * p, integer * emin,
 857         logical * ieee, integer * emax, real * rmax)
 858 {
 859     /* System generated locals */
 860     integer i__1;
 861     real r__1;
 862
 863     /* Local variables */
 864     static integer i__;
 865     static real y, z__;
 866     static integer try__, lexp;
 867     static real oldy;
 868     static integer uexp, nbits;
 869     extern doublereal slamc3_(real *, real *);
 870     static real recbas;
 871     static integer exbits, expsum;
 872
 873
 874 /*  -- LAPACK auxiliary routine (version 3.0) -- */
 875 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
 876 /*     Courant Institute, Argonne National Lab, and Rice University */
 877 /*     October 31, 1992 */
 878
 879 /*     .. Scalar Arguments .. */
 880 /*     .. */
 881
 882 /*  Purpose */
 883 /*  ======= */
 884
 885 /*  SLAMC5 attempts to compute RMAX, the largest machine floating-point */
 886 /*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum */
 887 /*  approximately to a power of 2.  It will fail on machines where this */
 888 /*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
 889 /*  EMAX = 28718).  It will also fail if the value supplied for EMIN is */
 890 /*  too large (i.e. too close to zero), probably with overflow. */
 891
 892 /*  Arguments */
 893 /*  ========= */
 894
 895 /*  BETA    (input) INTEGER */
 896 /*          The base of floating-point arithmetic. */
 897
 898 /*  P       (input) INTEGER */
 899 /*          The number of base BETA digits in the mantissa of a */
 900 /*          floating-point value. */
 901
 902 /*  EMIN    (input) INTEGER */
 903 /*          The minimum exponent before (gradual) underflow. */
 904
 905 /*  IEEE    (input) LOGICAL */
 906 /*          A logical flag specifying whether or not the arithmetic */
 907 /*          system is thought to comply with the IEEE standard. */
 908
 909 /*  EMAX    (output) INTEGER */
 910 /*          The largest exponent before overflow */
 911
 912 /*  RMAX    (output) REAL */
 913 /*          The largest machine floating-point number. */
 914
 915 /* ===================================================================== */
 916
 917 /*     .. Parameters .. */
 918 /*     .. */
 919 /*     .. Local Scalars .. */
 920 /*     .. */
 921 /*     .. External Functions .. */
 922 /*     .. */
 923 /*     .. Intrinsic Functions .. */
 924 /*     .. */
 925 /*     .. Executable Statements .. */
 926
 927 /*     First compute LEXP and UEXP, two powers of 2 that bound */
 928 /*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
 929 /*     approximately to the bound that is closest to abs(EMIN). */
 930 /*     (EMAX is the exponent of the required number RMAX). */
 931
 932     lexp = 1;
 933     exbits = 1;
 934   L10:
 935     try__ = lexp << 1;
 936     if (try__ <= -(*emin)) {
 937         lexp = try__;
 938         ++exbits;
 939         goto L10;
 940     }
 941     if (lexp == -(*emin)) {
 942         uexp = lexp;
 943     }
 944     else {
 945         uexp = try__;
 946         ++exbits;
 947     }
 948
 949 /*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
 950 /*     than or equal to EMIN. EXBITS is the number of bits needed to */
 951 /*     store the exponent. */
 952
 953     if (uexp + *emin > -lexp - *emin) {
 954         expsum = lexp << 1;
 955     }
 956     else {
 957         expsum = uexp << 1;
 958     }
 959
 960 /*     EXPSUM is the exponent range, approximately equal to */
 961 /*     EMAX - EMIN + 1 . */
 962
 963     *emax = expsum + *emin - 1;
 964     nbits = exbits + 1 + *p;
 965
 966 /*     NBITS is the total number of bits needed to store a */
 967 /*     floating-point number. */
 968
 969     if (nbits % 2 == 1 && *beta == 2) {
 970
 971 /*        Either there are an odd number of bits used to store a */
 972 /*        floating-point number, which is unlikely, or some bits are */
 973 /*        not used in the representation of numbers, which is possible, */
 974 /*        (e.g. Cray machines) or the mantissa has an implicit bit, */
 975 /*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
 976 /*        most likely. We have to assume the last alternative. */
 977 /*        If this is true, then we need to reduce EMAX by one because */
 978 /*        there must be some way of representing zero in an implicit-bit */
 979 /*        system. On machines like Cray, we are reducing EMAX by one */
 980 /*        unnecessarily. */
 981
 982         --(*emax);
 983     }
 984
 985     if (*ieee) {
 986
 987 /*        Assume we are on an IEEE machine which reserves one exponent */
 988 /*        for infinity and NaN. */
 989
 990         --(*emax);
 991     }
 992
 993 /*     Now create RMAX, the largest machine number, which should */
 994 /*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
 995
 996 /*     First compute 1.0 - BETA**(-P), being careful that the */
 997 /*     result is less than 1.0 . */
 998
 999     recbas = 1.f / *beta;
1000     z__ = *beta - 1.f;
1001     y = 0.f;
1002     i__1 = *p;
1003     for (i__ = 1; i__ <= i__1; ++i__) {
1004         z__ *= recbas;
1005         if (y < 1.f) {
1006             oldy = y;
1007         }
1008         y = slamc3_(&y, &z__);
1009 /* L20: */
1010     }
1011     if (y >= 1.f) {
1012         y = oldy;
1013     }
1014
1015 /*     Now multiply by BETA**EMAX to get RMAX. */
1016
1017     i__1 = *emax;
1018     for (i__ = 1; i__ <= i__1; ++i__) {
1019         r__1 = y * *beta;
1020         y = slamc3_(&r__1, &c_b32);
1021 /* L30: */
1022     }
1023
1024     *rmax = y;
1025     return 0;
1026
1027 /*     End of SLAMC5 */
1028
1029 }                               /* slamc5_ */