function si (x)
c december 1980 edition,  w. fullerton, bell labs.
      dimension sics(12)
      external cos, csevl, inits, r1mach, sin, sqrt
c
c series for si   on the interval  0.00000e+00 to  1.60000e+01
c                                        with weighted error   1.22e-17
c                                         log weighted error  16.91
c                               significant figures required  16.37
c                                    decimal places required  17.45
c
      data si  cs(  1) /    -0.1315646598 184841929 e0/
      data si  cs(  2) /    -0.2776578526 973601892 e0/
      data si  cs(  3) /     0.0354414054 866659180 e0/
      data si  cs(  4) /    -0.0025631631 447933978 e0/
      data si  cs(  5) /     0.0001162365 390497009 e0/
      data si  cs(  6) /    -0.0000035904 327241606 e0/
      data si  cs(  7) /     0.0000000802 342123706 e0/
      data si  cs(  8) /    -0.0000000013 562997693 e0/
      data si  cs(  9) /     0.0000000000 179440722 e0/
      data si  cs( 10) /    -0.0000000000 001908387 e0/
      data si  cs( 11) /     0.0000000000 000016670 e0/
      data si  cs( 12) /    -0.0000000000 000000122 e0/
c
      data pi2 / 1.570796326 7948966e0 /
      data nsi, xsml /0, 0.0/
c
      if (nsi.ne.0) go to 10
      nsi = inits (sics, 12, 0.1*r1mach(3))
      xsml = sqrt (eps)
c
 10   absx = abs(x)
      if (absx.gt.4.0) go to 20
      si = x
      if (absx.lt.xsml) return
c
      si = x*(0.75 + csevl ((x*x-8.0)*0.125, sics, nsi))
      return
c
 20   call r9sifg (absx, f, g)
      cosx = cos (absx)
      call erroff
      si = pi2 - f*cosx - g*sin(x)
      if (x.lt.0.0) si = -si
c
      return
      end