d9pak.f
1.29 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
double precision function d9pak (y, n)
c december 1979 edition. w. fullerton, c3, los alamos scientific lab.
c
c pack a base 2 exponent into floating point number x. this routine is
c almost the inverse of d9upak. it is not exactly the inverse, because
c dabs(x) need not be between 0.5 and 1.0. if both d9pak and 2.d0**n
c were known to be in range we could compute
c d9pak = x * 2.0d0**n
c
double precision y, aln2b, aln210, d1mach
external d1mach, i1mach
data nmin, nmax / 2*0 /
data aln210 / 3.321928094 8873623478 7031942948 9 d0 /
c
if (nmin.ne.0) go to 10
aln2b = 1.0d0
if (i1mach(10).ne.2) aln2b = d1mach(5)*aln210
nmin = aln2b*dble(float(i1mach(15)))
nmax = aln2b*dble(float(i1mach(16)))
c
10 call d9upak (y, d9pak, ny)
c
nsum = n + ny
if (nsum.lt.nmin) go to 40
if (nsum.gt.nmax) call seteru (
1 31hd9pak packed number overflows, 31, 1, 2)
c
if (nsum.eq.0) return
if (nsum.gt.0) go to 30
c
20 d9pak = 0.5d0*d9pak
nsum = nsum + 1
if (nsum.ne.0) go to 20
return
c
30 d9pak = 2.0d0*d9pak
nsum = nsum - 1
if (nsum.ne.0) go to 30
return
c
40 call seteru (32hd9pak packed number underflows, 32, 1, 0)
d9pak = 0.0d0
return
c
end