docs: def _erf(x):
      """
      Port of cephes ``ndtr.c`` ``erf`` function.
  
      See https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtr.c
      """
      T = [
          9.60497373987051638749E0,
          9.00260197203842689217E1,
          2.23200534594684319226E3,
          7.00332514112805075473E3,
          5.55923013010394962768E4,
      ]
  
      U = [
          3.35617141647503099647E1,
          5.21357949780152679795E2,
          4.59432382970980127987E3,
          2.26290000613890934246E4,
          4.92673942608635921086E4,
      ]
  
      # Shorcut special cases
      if x == 0:
          return 0
      if x >= MAXVAL:
          return 1
      if x <= -MAXVAL:
          return -1
  
      if abs(x) > 1:
          return 1 - erfc(x)
  
      z = x * x
      return x * _polevl(z, T, 4) / _p1evl(z, U, 5)