double precision function dasin (x)
c may 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, asincs(39), pi2, sqeps, y, z, dcsevl,
     1  d1mach, dsqrt
      external d1mach, dcsevl, dsqrt, initds
c
c series for asin       on the interval  0.          to  5.00000e-01
c                                        with weighted error   1.62e-32
c                                         log weighted error  31.79
c                               significant figures required  30.67
c                                    decimal places required  32.59
c
      data asincs(  1) / +.1024639175 3227159336 5731483057 85 d+0     /
      data asincs(  2) / +.5494648722 1245833306 0111959029 24 d-1     /
      data asincs(  3) / +.4080630392 5449692851 3070561492 46 d-2     /
      data asincs(  4) / +.4078900685 4604435455 5988239056 12 d-3     /
      data asincs(  5) / +.4698536743 2203691616 0485301362 18 d-4     /
      data asincs(  6) / +.5880975813 9708058986 4543855520 74 d-5     /
      data asincs(  7) / +.7773231246 2777632750 5575281637 95 d-6     /
      data asincs(  8) / +.1067742334 0082039235 0475049565 87 d-6     /
      data asincs(  9) / +.1509239953 6022808262 3864344010 64 d-7     /
      data asincs( 10) / +.2180972408 0055385496 6096147139 30 d-8     /
      data asincs( 11) / +.3207598426 2789614433 2619596673 76 d-9     /
      data asincs( 12) / +.4785536964 6781034461 4931339189 53 d-10    /
      data asincs( 13) / +.7225128736 2910432263 8487545371 12 d-11    /
      data asincs( 14) / +.1101833474 2255783705 3727013349 87 d-11    /
      data asincs( 15) / +.1694763253 9203354877 4237456510 78 d-12    /
      data asincs( 16) / +.2626155866 7348224162 2832415024 16 d-13    /
      data asincs( 17) / +.4095829981 3281178408 8280692911 10 d-14    /
      data asincs( 18) / +.6424479310 8803655891 7279448870 91 d-15    /
      data asincs( 19) / +.1012814219 8228221693 9733612220 41 d-15    /
      data asincs( 20) / +.1603922189 7380787560 0505974647 46 d-16    /
      data asincs( 21) / +.2550350135 5807141715 2987896763 73 d-17    /
      data asincs( 22) / +.4070140379 7862382855 4871656721 06 d-18    /
      data asincs( 23) / +.6517267171 2881144437 8892675754 66 d-19    /
      data asincs( 24) / +.1046745303 7096796954 2448917162 66 d-19    /
      data asincs( 25) / +.1685872556 3380328094 9890951850 66 d-20    /
      data asincs( 26) / +.2722193630 5040227625 1643412479 99 d-21    /
      data asincs( 27) / +.4405929390 0347550617 1268300799 99 d-22    /
      data asincs( 28) / +.7146668524 3375937853 0631680000 00 d-23    /
      data asincs( 29) / +.1161579334 3859516051 7989717333 33 d-23    /
      data asincs( 30) / +.1891523455 2354685801 1841877333 33 d-24    /
      data asincs( 31) / +.3085577204 4244342399 8279680000 00 d-25    /
      data asincs( 32) / +.5041636602 2162453412 9704959999 99 d-26    /
      data asincs( 33) / +.8250272550 2400865081 7536000000 00 d-27    /
      data asincs( 34) / +.1352003263 1020947208 0554666666 66 d-27    /
      data asincs( 35) / +.2218432687 6541720216 6442666666 66 d-28    /
      data asincs( 36) / +.3644249405 4085079212 5781333333 33 d-29    /
      data asincs( 37) / +.5992021855 8643813307 7333333333 33 d-30    /
      data asincs( 38) / +.9858481205 9573785810 2613333333 33 d-31    /
      data asincs( 39) / +.1622250116 6399014393 1733333333 33 d-31    /
c
      data pi2 / 1.5707963267 9489661923 1321691639 75 d0 /
      data nterms, sqeps / 0, 0.d0 /
c
      if (nterms.ne.0) go to 10
      nterms = initds (asincs, 39, 0.1*sngl(d1mach(3)) )
      sqeps = dsqrt (6.d0*d1mach(3))
c
 10   y = dabs (x)
      if (y.gt.1.0d0) go to 20
c
      z = 0.d0
      if (y.gt.sqeps) z = y*y
      if (z.le.0.5d0) dasin = x*(1.0d0 + dcsevl (4.d0*z-1.d0, asincs,
     1  nterms))
      if (z.gt.0.5d0) dasin = pi2 - dsqrt(1.d0-z) * (1.d0 +
     1  dcsevl (3.0d0-4.0d0*z, asincs, nterms))
c
      if (x.ne.0.d0) dasin = dsign (dasin, x)
      return
c
 20   if (y.gt.1.0d0+sqeps) call seteru (
     1  33hdasin   abs(x) greater than 1+eta, 33, 2, 2)
      call seteru (33hdasin   abs(x) gt 1, but lt 1+eta, 33, 1, 1)
      dasin = dsign (pi2, x)
      return
c
      end