Function: math-find-root

math-find-root is an autoloaded and byte-compiled function defined in calcalg3.el.gz.

Signature

(math-find-root EXPR VAR GUESS ROOT-WIDEN)

Source Code

;; Defined in /usr/src/emacs/lisp/calc/calcalg3.el.gz
(defun math-find-root (expr var guess root-widen)
  (let ((math-root-widen root-widen))
  (if (eq (car-safe expr) 'vec)
      (let ((n (1- (length expr)))
	    (calc-symbolic-mode nil)
	    (var-DUMMY nil)
	    (jacob (list 'vec))
	    p p2 m row)
	(unless (eq (car-safe var) 'vec)
	  (math-reject-arg var 'vectorp))
	(unless (= (length var) (1+ n))
	  (math-dimension-error))
	(setq expr (copy-sequence expr))
	(while (>= n (length math-root-vars))
	  (let ((symb (intern (concat "math-root-v"
				      (int-to-string
				       (length math-root-vars))))))
	    (setq math-root-vars (vconcat math-root-vars
					  (vector (list 'var symb symb))))))
	(setq m -1)
	(while (< (setq m (1+ m)) n)
	  (set (nth 2 (aref math-root-vars m)) nil))
	(setq m -1 p var)
	(while (setq m (1+ m) p (cdr p))
	  (or (eq (car-safe (car p)) 'var)
	      (math-reject-arg var "*Expected a variable"))
	  (setq p2 expr)
	  (while (setq p2 (cdr p2))
	    (setcar p2 (math-expr-subst (car p2) (car p)
					(aref math-root-vars m)))))
	(unless (eq (car-safe guess) 'vec)
	  (math-reject-arg guess 'vectorp))
	(unless (= (length guess) (1+ n))
	  (math-dimension-error))
	(setq guess (copy-sequence guess)
	      p guess)
	(while (setq p (cdr p))
	  (or (Math-numberp (car guess))
	      (math-reject-arg guess 'numberp))
	  (setcar p (math-float (car p))))
	(setq p expr)
	(while (setq p (cdr p))
	  (if (assq (car-safe (car p)) calc-tweak-eqn-table)
	      (setcar p (math-sub (nth 1 (car p)) (nth 2 (car p)))))
	  (setcar p (math-evaluate-expr (car p)))
	  (setq row (list 'vec)
		m -1)
	  (while (< (setq m (1+ m)) n)
	    (nconc row (list (math-evaluate-expr
			      (or (calcFunc-deriv (car p)
						  (aref math-root-vars m)
						  nil t)
				  (math-reject-arg
				   expr
				   "*Formulas must be differentiable"))))))
	  (nconc jacob (list row)))
	(setq m (math-abs-approx guess))
	(math-newton-multi expr jacob n guess guess
			   (if (math-zerop m) '(float 1 3) (math-mul m 10))))
    (unless (eq (car-safe var) 'var)
      (math-reject-arg var "*Expected a variable"))
    (unless (math-expr-contains expr var)
      (math-reject-arg expr "*Formula does not contain specified variable"))
    (if (assq (car expr) calc-tweak-eqn-table)
	(setq expr (math-sub (nth 1 expr) (nth 2 expr))))
    (math-with-extra-prec 2
      (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
      (let* ((calc-symbolic-mode nil)
	     (var-DUMMY nil)
	     (expr (math-evaluate-expr expr))
	     (deriv (calcFunc-deriv expr '(var DUMMY var-DUMMY) nil t))
	     low high vlow vhigh)
	(and deriv (setq deriv (math-evaluate-expr deriv)))
	(setq guess (math-float guess))
	(if (and (math-numberp guess)
		 deriv)
	    (math-newton-root expr deriv guess guess
			      (if (math-zerop guess) '(float 1 6)
				(math-mul (math-abs-approx guess) 100)))
	  (if (Math-realp guess)
	      (setq low guess
		    high guess
		    var-DUMMY guess
		    vlow (math-evaluate-expr expr)
		    vhigh vlow
		    math-root-widen 'point)
	    (if (eq (car guess) 'intv)
		(progn
		  (or (math-constp guess) (math-reject-arg guess 'constp))
		  (setq low (nth 2 guess)
			high (nth 3 guess))
		  (if (memq (nth 1 guess) '(0 1))
		      (setq low (calcFunc-incr low 1 high)))
		  (if (memq (nth 1 guess) '(0 2))
		      (setq high (calcFunc-incr high -1 low)))
		  (setq var-DUMMY low
			vlow (math-evaluate-expr expr)
			var-DUMMY high
			vhigh (math-evaluate-expr expr)))
	      (if (math-complexp guess)
		  (math-reject-arg "*Complex root finder must have derivative")
		(math-reject-arg guess 'realp))))
	  (if (Math-zerop vlow)
	      (list 'vec low vlow)
	    (if (Math-zerop vhigh)
		(list 'vec high vhigh)
	      (if (and deriv (Math-numberp vlow) (Math-numberp vhigh))
		  (math-newton-search-root expr deriv nil nil nil nil
					   low vlow high vhigh)
		(if (or (and (Math-posp vlow) (Math-posp vhigh))
			(and (Math-negp vlow) (Math-negp vhigh))
			(not (Math-numberp vlow))
			(not (Math-numberp vhigh)))
		    (math-search-root expr deriv low vlow high vhigh)
		  (math-bisect-root expr low vlow high vhigh)))))))))))