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