SphinxBase  5prealpha
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_ */
cilist
Definition: f2c.h:45