function asin (x)
c may 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension asincs(20)
      external csevl, inits, r1mach, sqrt
c
c series for asin       on the interval  0.          to  5.00000d-01
c                                        with weighted error   1.60e-17
c                                         log weighted error  16.79
c                               significant figures required  15.67
c                                    decimal places required  17.45
c
      data asincs( 1) /    .1024639175 3227159e0 /
      data asincs( 2) /    .0549464872 21245833e0 /
      data asincs( 3) /    .0040806303 92544969e0 /
      data asincs( 4) /    .0004078900 68546044e0 /
      data asincs( 5) /    .0000469853 67432203e0 /
      data asincs( 6) /    .0000058809 75813970e0 /
      data asincs( 7) /    .0000007773 23124627e0 /
      data asincs( 8) /    .0000001067 74233400e0 /
      data asincs( 9) /    .0000000150 92399536e0 /
      data asincs(10) /    .0000000021 80972408e0 /
      data asincs(11) /    .0000000003 20759842e0 /
      data asincs(12) /    .0000000000 47855369e0 /
      data asincs(13) /    .0000000000 07225128e0 /
      data asincs(14) /    .0000000000 01101833e0 /
      data asincs(15) /    .0000000000 00169476e0 /
      data asincs(16) /    .0000000000 00026261e0 /
      data asincs(17) /    .0000000000 00004095e0 /
      data asincs(18) /    .0000000000 00000642e0 /
      data asincs(19) /    .0000000000 00000101e0 /
      data asincs(20) /    .0000000000 00000016e0 /
c
      data pi2 / 1.5707963267 9489661923e0/
      data nterms, sqeps /0, 0./
c
      if (nterms.ne.0) go to 10
      nterms = inits (asincs, 20, 0.1*r1mach(3))
      sqeps = sqrt (6.0*r1mach(3))
c
 10   y = abs(x)
      if (y.gt.1.0) go to 20
c
      z = 0.
      if (y.gt.sqeps) z = y*y
      if (z.le.0.5) asin = x*(1. + csevl (4.*z-1., asincs, nterms))
      if (z.gt.0.5) asin = pi2 - sqrt(1.-z)*(1. + csevl (3.-4.*z,
     1  asincs, nterms))
      if (x.ne.0.) asin = sign (asin, x)
      return
c
 20   if (y.gt.1.0+sqeps) call seteru (
     1  33hasin    abs(x) greater than 1+eta, 33, 2, 2)
      call seteru (33hasin    abs(x) gt 1, but lt 1+eta, 33, 1, 1)
      asin = sign (pi2, x)
      return
c
      end