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
21static integer c__1 = 1;
22static real c_b32 = 0.f;
23
24doublereal
25slamch_(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
183slamc1_(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
392slamc2_(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
713doublereal
714slamc3_(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
757slamc4_(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
856slamc5_(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_ */
Definition: f2c.h:46