function sandom(r)
      integer ia1,ia0,ia1ma0,ic,ix1,ix0,iy0,iy1
      data ia1, ia0, ia1ma0 /1536,1029,507/
      data ic /1731/
      data ix1, ix0 /0,0/
c
c           a*x = 2**22*ia1*ix1 + 2**11*(ia1*ix1 + (ia1-ia0)*(ix0-ix1)
c                 + ia0*ix0) + ia0*ix0
c
      iy0=ia0*ix0
      iy1=ia1*ix1+ia1ma0*(ix0-ix1)+iy0
      iy0=iy0+ic
      ix0=mod(iy0,2048)
      iy1=iy1+(iy0-ix0)/2048
      ix1=mod(iy1,2048)
c
      entry sandsv(r)
   10 sandom=ix1*2048+ix0
      sandom=sandom/4194304.
      r=sandom
      return
c
      entry sandst(r)
      ix1=amod(r,1.)*4194304.+0.5
      ix0=mod(ix1,2048)
      ix1=(ix1-ix0)/2048
      go to 10
c
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*