share/numeric/newton.mac

   1
   2 newton(expr,guess):=block([y,prec:10^(fpprec/-2.0b0),var,dif, listconstvars : false],
   3        var:listofvars(expr)[1],
   4        dif:diff(expr,var),
   5        guess:bfloat(guess),
   6        do (y:sublis([var=guess],dif),y:expand(y),if abs(y) < prec then error("derivative at",guess,"is zero"),
   7            guess:guess-sublis([var=(guess)],expr)/y,guess:rectform(expand(guess)),
   8            if abs(sublis([var=(guess)],expr)) < prec then return(guess)));
   9
  10