lang/pc/test/machar.p

   1 { $Id: machar.p,v 2.2 1994/06/24 12:36:47 ceriel Exp $ }
   2
   3 procedure machar (var ibeta , it , irnd , ngrd , machep , negep , iexp,
   4   minexp , maxexp : integer ; var eps , epsneg , xmin , xmax : real ) ;
   5 var trapped:boolean;
   6
   7 procedure encaps(procedure p; procedure q(i:integer)); extern;
   8 procedure trap(i:integer); extern;
   9
  10 procedure catch(i:integer);
  11 const underflo=5;
  12 begin if i=underflo then trapped:=true else trap(i) end;
  13
  14 procedure work;
  15 var
  16
  17
  18 {     This subroutine is intended to determine the characteristics
  19       of the floating-point arithmetic system that are specified
  20       below.  The first three are determined according to an
  21       algorithm due to M. Malcolm, CACM 15 (1972), pp. 949-951,
  22       incorporating some, but not all, of the improvements
  23       suggested by M. Gentleman and S. Marovich, CACM 17 (1974),
  24       pp. 276-277.  The version given here is for single precision.
  25
  26       Latest revision - October 1, 1976.
  27
  28       Author - W. J. Cody
  29                Argonne National Laboratory
  30
  31       Revised for Pascal - R. A. Freak
  32                            University of Tasmania
  33                            Hobart
  34                            Tasmania
  35
  36       ibeta    is the radix of the floating-point representation
  37       it       is the number of base ibeta digits in the floating-point
  38                   significand
  39       irnd     =  0 if the arithmetic chops,
  40                   1 if the arithmetic rounds
  41       ngrd     =  0 if  irnd=1, or if  irnd=0  and only  it  base ibeta
  42                     digits participate in the post normalization shift
  43                     of the floating-point significand in multiplication
  44                   1 if  irnd=0  and more than  it  base  ibeta  digits
  45                     participate in the post normalization shift of the
  46                     floating-point significand in multiplication
  47       machep   is the exponent on the smallest positive floating-point
  48                   number  eps such that  1.0+eps <> 1.0
  49       negeps   is the exponent on the smallest positive fl. pt. no.
  50                   negeps such that  1.0-negeps <> 1.0, except that
  51                   negeps is bounded below by  it-3
  52       iexp     is the number of bits (decimal places if ibeta = 10)
  53                   reserved for the representation of the exponent of
  54                   a floating-point number
  55       minexp   is the exponent of the smallest positive fl. pt. no.
  56                   xmin
  57       maxexp   is the exponent of the largest finite floating-point
  58                   number  xmax
  59       eps      is the smallest positive floating-point number such
  60                   that  1.0+eps <> 1.0. in particular,
  61                   eps = ibeta**machep
  62       epsneg   is the smallest positive floating-point number such
  63                   that  1.0-eps <> 1.0  (except that the exponent
  64                   negeps is bounded below by it-3).  in particular
  65                   epsneg = ibeta**negep
  66       xmin     is the smallest positive floating-point number.  in
  67                   particular,  xmin = ibeta ** minexp
  68       xmax     is the largest finite floating-point number.  in
  69                   particular   xmax = (1.0-epsneg) * ibeta ** maxexp
  70                   note - on some machines  xmax  will be only the
  71                   second, or perhaps third, largest number, being
  72                   too small by 1 or 2 units in the last digit of
  73                   the significand.
  74
  75                                                                     }
  76
  77    i , iz , j , k , mx : integer ;
  78    a , b , beta , betain , betam1 , one , y , z , zero : real ;
  79
  80 begin
  81    irnd := 1 ;
  82    one := ( irnd );
  83    a := one + one ;
  84    b := a ;
  85    zero := 0.0 ;
  86 {
  87       determine ibeta,beta ala Malcolm
  88                                                                     }
  89    while ( ( ( a + one ) - a ) - one = zero ) do begin
  90       a := a + a ;
  91    end ;
  92    while ( ( a + b ) - a = zero ) do begin
  93       b := b + b ;
  94    end ;
  95    ibeta := trunc ( ( a + b ) - a );
  96    beta := ( ibeta );
  97    betam1 := beta - one ;
  98 {
  99       determine irnd,ngrd,it
 100                                                                     }
 101    if ( ( a + betam1 ) - a = zero ) then irnd := 0 ;
 102    it := 0 ;
 103    a := one ;
 104    repeat begin
 105       it := it + 1 ;
 106       a := a * beta ;
 107    end until ( ( ( a + one ) - a ) - one <> zero ) ;
 108 {
 109       determine negep, epsneg
 110                                                                     }
 111    negep := it + 3 ;
 112    a := one ;
 113
 114    for i := 1 to negep do begin
 115       a := a / beta ;
 116    end ;
 117
 118    while ( ( one - a ) - one = zero ) do begin
 119       a := a * beta ;
 120       negep := negep - 1 ;
 121    end ;
 122    negep := - negep ;
 123    epsneg := a ;
 124 {
 125       determine machep, eps
 126                                                                     }
 127    machep := negep ;
 128    while ( ( one + a ) - one = zero ) do begin
 129       a := a * beta ;
 130       machep := machep + 1 ;
 131    end ;
 132    eps := a ;
 133 {
 134       determine ngrd
 135                                                                     }
 136    ngrd := 0 ;
 137    if(( irnd = 0) and((( one + eps) * one - one) <> zero)) then
 138    ngrd := 1 ;
 139 {
 140       determine iexp, minexp, xmin
 141
 142       loop to determine largest i such that
 143           (1/beta) ** (2**(i))
 144       does not underflow
 145       exit from loop is signall by an underflow
 146                                                                     }
 147    i := 0 ;
 148    betain := one / beta ;
 149    z := betain ;
 150    trapped:=false;
 151    repeat begin
 152       y := z ;
 153       z := y * y ;
 154 {
 155       check for underflow
 156                                                                     }
 157       i := i + 1 ;
 158    end until trapped;
 159    i := i - 1;
 160    k := 1 ;
 161 {
 162       determine k such that (1/beta)**k does not underflow
 163
 164       first set k = 2 ** i
 165                                                                     }
 166
 167    for j := 1 to i do begin
 168       k := k + k ;
 169    end ;
 170
 171    iexp := i + 1 ;
 172    mx := k + k ;
 173    if ( ibeta = 10 ) then begin
 174 {
 175       for decimal machines only                                     }
 176       iexp := 2 ;
 177       iz := ibeta ;
 178       while ( k >= iz ) do begin
 179          iz := iz * ibeta ;
 180          iexp := iexp + 1 ;
 181       end ;
 182       mx := iz + iz - 1 ;
 183    end;
 184    trapped:=false;
 185    repeat begin
 186 {
 187       loop to construct xmin
 188       exit from loop is signalled by an underflow
 189                                                                     }
 190       xmin := y ;
 191       y := y * betain ;
 192       k := k + 1 ;
 193    end until trapped;
 194    k := k - 1;
 195    minexp := - k ;
 196 {  determine maxexp, xmax
 197                                                                     }
 198    if ( ( mx <= k + k - 3 ) and ( ibeta <> 10 ) ) then begin
 199       mx := mx + mx ;
 200       iexp := iexp + 1 ;
 201    end;
 202    maxexp := mx + minexp ;
 203 {  adjust for machines with implicit leading
 204    bit in binary significand and machines with
 205    radix point at extreme right of significand
 206                                                                     }
 207    i := maxexp + minexp ;
 208    if ( ( ibeta = 2 ) and ( i = 0 ) ) then maxexp := maxexp - 1 ;
 209    if ( i > 20 ) then maxexp := maxexp - 3 ;
 210    xmax := one - epsneg ;
 211    if ( xmax * one <> xmax ) then xmax := one - beta * epsneg ;
 212    xmax := ( xmax * betain * betain * betain ) / xmin ;
 213    i := maxexp + minexp + 3 ;
 214    if  ( i > 0 ) then begin
 215
 216       for j := 1 to i do begin
 217          xmax := xmax * beta ;
 218       end ;
 219    end;
 220
 221 end;
 222
 223 begin
 224   trapped:=false;
 225   encaps(work,catch);
 226 end;