3Â¥: Originally posted by Â̱éɽԭ at 2013-11-12 16:23:55
Ò²¿ÉÒÔÓÃimsl

2Â¥: Originally posted by Â̱éɽԭ at 2013-11-12 16:21:46
cvfûÓÐÎó²îº¯Êý£¬×Ô¼ºÕÒÒ»¸ö¼Óµ½Ô´´úÂëÀïÃæ¡£

4Â¥: Originally posted by sun805869434 at 2013-11-12 16:36:39
ÇëÎÊÔ´´úÂëÄÄÀïÄÜÕÒµ½ÄØ£¿...

function r8_erf ( x )

!*****************************************************************************80
!
!! R8_ERF evaluates the error function of an R8 argument.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    10 September 2011
!
!  Author:
!
!    Original FORTRAN77 version by Wayne Fullerton.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Wayne Fullerton,
!    Portable Special Function Routines,
!    in Portability of Numerical Software,
!    edited by Wayne Cowell,
!    Lecture Notes in Computer Science, Volume 57,
!    Springer 1977,
!    ISBN: 978-3-540-08446-4,
!    LC: QA297.W65.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument.
!
!    Output, real ( kind = 8 ) R8_ERF, the error function of X.
!
  implicit none

  real ( kind = 8 ) erfcs(21)
  integer ( kind = 4 ) nterf
  real ( kind = 8 ) r8_csevl
  real ( kind = 8 ) r8_erf
  real ( kind = 8 ) r8_erfc
  integer ( kind = 4 ) r8_inits
  real ( kind = 8 ) r8_mach
  real ( kind = 8 ), save :: sqeps = 0.0D+00
  real ( kind = 8 ) sqrtpi
  real ( kind = 8 ) value
  real ( kind = 8 ) x
  real ( kind = 8 ) xbig
  real ( kind = 8 ) y

  save erfcs
  save nterf
  save sqrtpi
  save xbig

  data erfcs(  1) / -0.49046121234691808039984544033376D-01 /
  data erfcs(  2) / -0.14226120510371364237824741899631D+00 /
  data erfcs(  3) / +0.10035582187599795575754676712933D-01 /
  data erfcs(  4) / -0.57687646997674847650827025509167D-03 /
  data erfcs(  5) / +0.27419931252196061034422160791471D-04 /
  data erfcs(  6) / -0.11043175507344507604135381295905D-05 /
  data erfcs(  7) / +0.38488755420345036949961311498174D-07 /
  data erfcs(  8) / -0.11808582533875466969631751801581D-08 /
  data erfcs(  9) / +0.32334215826050909646402930953354D-10 /
  data erfcs( 10) / -0.79910159470045487581607374708595D-12 /
  data erfcs( 11) / +0.17990725113961455611967245486634D-13 /
  data erfcs( 12) / -0.37186354878186926382316828209493D-15 /
  data erfcs( 13) / +0.71035990037142529711689908394666D-17 /
  data erfcs( 14) / -0.12612455119155225832495424853333D-18 /
  data erfcs( 15) / +0.20916406941769294369170500266666D-20 /
  data erfcs( 16) / -0.32539731029314072982364160000000D-22 /
  data erfcs( 17) / +0.47668672097976748332373333333333D-24 /
  data erfcs( 18) / -0.65980120782851343155199999999999D-26 /
  data erfcs( 19) / +0.86550114699637626197333333333333D-28 /
  data erfcs( 20) / -0.10788925177498064213333333333333D-29 /
  data erfcs( 21) / +0.12811883993017002666666666666666D-31 /

  data nterf / 0 /
  data sqrtpi / 1.77245385090551602729816748334115D+00 /
  data xbig / 0.0D+00 /

  if ( nterf == 0 ) then
    nterf = r8_inits ( erfcs, 21, 0.1D+00 * r8_mach ( 3 ) )
    xbig = sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) )
    sqeps = sqrt ( 2.0D+00 * r8_mach ( 3 ) )
  end if

  y = abs ( x )

  if ( y <= sqeps ) then
    value = 2.0D+00 * x / sqrtpi
  else if ( y <= 1.0D+00 ) then
    value = x * ( 1.0D+00 &
      + r8_csevl ( 2.0D+00 * x * x - 1.0D+00, erfcs, nterf ) )
  else if ( y <= xbig ) then
    value = 1.0D+00 - r8_erfc ( y )
    if ( x < 0.0D+00 ) then
      value = - value
    end if
  else
    value = 1.0D+00
    if ( x < 0.0D+00 ) then
      value = - value
    end if
  end if

  r8_erf = value

  return
end