Nerf3(a:FLOAT,ndigs:PI,nterms:PI):FLOAT == 
-- evaluate the error function 
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 7.1.5 of Abramowitz
  table.0:=a
  for i in 1..nterms repeat table.i:=(-1)*a^2*(2*i-1)*table.(i-1)/(i*(2*i+1))
  ans:=(2/sqrt(%pi)) * reduce(+,parts table)
  digits olds
  ans


Nerff(a:FLOAT,ndigs:PI,nterms:PI):FLOAT ==
-- in terms of onefone 7.1.21 of Abramowitz
  olds:= digits ndigs
  ans:= 2*a*exp(-a*a)*onefone(1.0,1.5,a*a,ndigs,nterms)/sqrt(%pi)
  digits olds
  ans



onefone(a:FLOAT,b:FLOAT,z:FLOAT,ndigs:PI,nterms:PI):FLOAT == 
-- the 1F1 confluent hypergeometric function
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 13.1.2 of Abramowitz
  table.0:=1
  table.1:=a*z/b
  for i in 2..nterms repeat table.i:=(a+i-1)*z*table.(i-1)/((b+i-1)*i)
  ans:=reduce(+,parts table)
  digits olds
  ans

expr:=erf(1.0)
tower %
% 1
erfop:=operator %
evl:LIST EXPR FLOAT -> EXPR FLOAT
-- insert an evaluation rule that always uses 50 digits and 500 terms
evl(x) == Nerff(x(1),50,500)
--attach it to the erf operator
evaluate(erfop,evl)$BOP1(EXPR FLOAT)
erf (1.0)
expr
numeric %

integrate(exp(-x^2),x=0..1)
numeric %