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