home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / dlamc1.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  188 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
  5. *
  6. *  -- LAPACK auxiliary routine (version 2.0) --
  7. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  8. *     Courant Institute, Argonne National Lab, and Rice University
  9. *     October 31, 1992
  10. *
  11. *     .. Scalar Arguments ..
  12.       LOGICAL            IEEE1, RND
  13.       INTEGER            BETA, T
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DLAMC1 determines the machine parameters given by BETA, T, RND, and
  20. *  IEEE1.
  21. *
  22. *  Arguments
  23. *  =========
  24. *
  25. *  BETA    (output) INTEGER
  26. *          The base of the machine.
  27. *
  28. *  T       (output) INTEGER
  29. *          The number of ( BETA ) digits in the mantissa.
  30. *
  31. *  RND     (output) LOGICAL
  32. *          Specifies whether proper rounding  ( RND = .TRUE. )  or
  33. *          chopping  ( RND = .FALSE. )  occurs in addition. This may not
  34. *          be a reliable guide to the way in which the machine performs
  35. *          its arithmetic.
  36. *
  37. *  IEEE1   (output) LOGICAL
  38. *          Specifies whether rounding appears to be done in the IEEE
  39. *          'round to nearest' style.
  40. *
  41. *  Further Details
  42. *  ===============
  43. *
  44. *  The routine is based on the routine  ENVRON  by Malcolm and
  45. *  incorporates suggestions by Gentleman and Marovich. See
  46. *
  47. *     Malcolm M. A. (1972) Algorithms to reveal properties of
  48. *        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
  49. *
  50. *     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
  51. *        that reveal properties of floating point arithmetic units.
  52. *        Comms. of the ACM, 17, 276-277.
  53. *
  54. * =====================================================================
  55. *
  56. *     .. Local Scalars ..
  57.       LOGICAL            FIRST, LIEEE1, LRND
  58.       INTEGER            LBETA, LT
  59.       DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
  60. *     ..
  61. *     .. External Functions ..
  62.       DOUBLE PRECISION   DLAMC3
  63.       EXTERNAL           DLAMC3
  64. *     ..
  65. *     .. Save statement ..
  66.       SAVE               FIRST, LIEEE1, LBETA, LRND, LT
  67. *     ..
  68. *     .. Data statements ..
  69.       DATA               FIRST / .TRUE. /
  70. *     ..
  71. *     .. Executable Statements ..
  72. *
  73.       IF( FIRST ) THEN
  74.          FIRST = .FALSE.
  75.          ONE = 1
  76. *
  77. *        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
  78. *        IEEE1, T and RND.
  79. *
  80. *        Throughout this routine  we use the function  DLAMC3  to ensure
  81. *        that relevant values are  stored and not held in registers,  or
  82. *        are not affected by optimizers.
  83. *
  84. *        Compute  a = 2.0**m  with the  smallest positive integer m such
  85. *        that
  86. *
  87. *           fl( a + 1.0 ) = a.
  88. *
  89.          A = 1
  90.          C = 1
  91. *
  92. *+       WHILE( C.EQ.ONE )LOOP
  93.    10    CONTINUE
  94.          IF( C.EQ.ONE ) THEN
  95.             A = 2*A
  96.             C = DLAMC3( A, ONE )
  97.             C = DLAMC3( C, -A )
  98.             GO TO 10
  99.          END IF
  100. *+       END WHILE
  101. *
  102. *        Now compute  b = 2.0**m  with the smallest positive integer m
  103. *        such that
  104. *
  105. *           fl( a + b ) .gt. a.
  106. *
  107.          B = 1
  108.          C = DLAMC3( A, B )
  109. *
  110. *+       WHILE( C.EQ.A )LOOP
  111.    20    CONTINUE
  112.          IF( C.EQ.A ) THEN
  113.             B = 2*B
  114.             C = DLAMC3( A, B )
  115.             GO TO 20
  116.          END IF
  117. *+       END WHILE
  118. *
  119. *        Now compute the base.  a and c  are neighbouring floating point
  120. *        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
  121. *        their difference is beta. Adding 0.25 to c is to ensure that it
  122. *        is truncated to beta and not ( beta - 1 ).
  123. *
  124.          QTR = ONE / 4
  125.          SAVEC = C
  126.          C = DLAMC3( C, -A )
  127.          LBETA = C + QTR
  128. *
  129. *        Now determine whether rounding or chopping occurs,  by adding a
  130. *        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
  131. *
  132.          B = LBETA
  133.          F = DLAMC3( B / 2, -B / 100 )
  134.          C = DLAMC3( F, A )
  135.          IF( C.EQ.A ) THEN
  136.             LRND = .TRUE.
  137.          ELSE
  138.             LRND = .FALSE.
  139.          END IF
  140.          F = DLAMC3( B / 2, B / 100 )
  141.          C = DLAMC3( F, A )
  142.          IF( ( LRND ) .AND. ( C.EQ.A ) )
  143.      $      LRND = .FALSE.
  144. *
  145. *        Try and decide whether rounding is done in the  IEEE  'round to
  146. *        nearest' style. B/2 is half a unit in the last place of the two
  147. *        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
  148. *        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
  149. *        A, but adding B/2 to SAVEC should change SAVEC.
  150. *
  151.          T1 = DLAMC3( B / 2, A )
  152.          T2 = DLAMC3( B / 2, SAVEC )
  153.          LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
  154. *
  155. *        Now find  the  mantissa, t.  It should  be the  integer part of
  156. *        log to the base beta of a,  however it is safer to determine  t
  157. *        by powering.  So we find t as the smallest positive integer for
  158. *        which
  159. *
  160. *           fl( beta**t + 1.0 ) = 1.0.
  161. *
  162.          LT = 0
  163.          A = 1
  164.          C = 1
  165. *
  166. *+       WHILE( C.EQ.ONE )LOOP
  167.    30    CONTINUE
  168.          IF( C.EQ.ONE ) THEN
  169.             LT = LT + 1
  170.             A = A*LBETA
  171.             C = DLAMC3( A, ONE )
  172.             C = DLAMC3( C, -A )
  173.             GO TO 30
  174.          END IF
  175. *+       END WHILE
  176. *
  177.       END IF
  178. *
  179.       BETA = LBETA
  180.       T = LT
  181.       RND = LRND
  182.       IEEE1 = LIEEE1
  183.       RETURN
  184. *
  185. *     End of DLAMC1
  186. *
  187.       END
  188.