/[emacs]/emacs/lisp/calc/calcalg3.el
ViewVC logotype

Diff of /emacs/lisp/calc/calcalg3.el

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4.4.2 by miles, Tue Oct 14 23:35:49 2003 UTC revision 1.4.4.3 by miles, Wed Dec 8 23:36:21 2004 UTC
# Line 3  Line 3 
3  ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.  ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
4    
5  ;; Author: David Gillespie <daveg@synaptics.com>  ;; Author: David Gillespie <daveg@synaptics.com>
6  ;; Maintainers: D. Goel <deego@gnufans.org>  ;; Maintainer: Jay Belanger <belanger@truman.edu>
 ;;              Colin Walters <walters@debian.org>  
7    
8  ;; This file is part of GNU Emacs.  ;; This file is part of GNU Emacs.
9    
# Line 28  Line 27 
27  ;;; Code:  ;;; Code:
28    
29  ;; This file is autoloaded from calc-ext.el.  ;; This file is autoloaded from calc-ext.el.
 (require 'calc-ext)  
30    
31    (require 'calc-ext)
32  (require 'calc-macs)  (require 'calc-macs)
33    
 (defun calc-Need-calc-alg-3 () nil)  
   
   
34  (defun calc-find-root (var)  (defun calc-find-root (var)
35    (interactive "sVariable(s) to solve for: ")    (interactive "sVariable(s) to solve for: ")
36    (calc-slow-wrapper    (calc-slow-wrapper
# Line 99  Line 95 
95         (calc-enter-result 1 "poli" (list 'calcFunc-polint data         (calc-enter-result 1 "poli" (list 'calcFunc-polint data
96                                           (calc-top 1)))))))                                           (calc-top 1)))))))
97    
98    ;; The variables calc-curve-nvars, calc-curve-varnames, calc-curve-model and calc-curve-coefnames are local to calc-curve-fit, but are
99    ;; used by calc-get-fit-variables which is called by calc-curve-fit.
100    (defvar calc-curve-nvars)
101    (defvar calc-curve-varnames)
102    (defvar calc-curve-model)
103    (defvar calc-curve-coefnames)
104    
105  (defun calc-curve-fit (arg &optional model coefnames varnames)  (defun calc-curve-fit (arg &optional calc-curve-model
106                               calc-curve-coefnames calc-curve-varnames)
107    (interactive "P")    (interactive "P")
108    (calc-slow-wrapper    (calc-slow-wrapper
109     (setq calc-aborted-prefix nil)     (setq calc-aborted-prefix nil)
# Line 108  Line 111 
111                   (if (calc-is-hyperbolic) 'calcFunc-efit                   (if (calc-is-hyperbolic) 'calcFunc-efit
112                     'calcFunc-fit)))                     'calcFunc-fit)))
113           key (which 0)           key (which 0)
114           n nvars temp data           n calc-curve-nvars temp data
115           (homog nil)           (homog nil)
116           (msgs '( "(Press ? for help)"           (msgs '( "(Press ? for help)"
117                    "1 = linear or multilinear"                    "1 = linear or multilinear"
# Line 120  Line 123 
123                    "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"                    "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
124                    "h prefix = homogeneous model (no constant term)"                    "h prefix = homogeneous model (no constant term)"
125                    "' = alg entry, $ = stack, u = Model1, U = Model2")))                    "' = alg entry, $ = stack, u = Model1, U = Model2")))
126       (while (not model)       (while (not calc-curve-model)
127         (message "Fit to model: %s:%s"         (message "Fit to model: %s:%s"
128                  (nth which msgs)                  (nth which msgs)
129                  (if homog " h" ""))                  (if homog " h" ""))
# Line 150  Line 153 
153                        (t (error "Bad prefix argument")))                        (t (error "Bad prefix argument")))
154                  (or (math-matrixp data) (not (cdr (cdr data)))                  (or (math-matrixp data) (not (cdr (cdr data)))
155                      (error "Data matrix is not a matrix!"))                      (error "Data matrix is not a matrix!"))
156                  (setq nvars (- (length data) 2)                  (setq calc-curve-nvars (- (length data) 2)
157                        coefnames nil                        calc-curve-coefnames nil
158                        varnames nil)                        calc-curve-varnames nil)
159                  nil))                  nil))
160               ((= key ?1)  ; linear or multilinear               ((= key ?1)  ; linear or multilinear
161                (calc-get-fit-variables nvars (1+ nvars) (and homog 0))                (calc-get-fit-variables calc-curve-nvars
162                (setq model (math-mul coefnames                                        (1+ calc-curve-nvars) (and homog 0))
163                                      (cons 'vec (cons 1 (cdr varnames))))))                (setq calc-curve-model (math-mul calc-curve-coefnames
164                                        (cons 'vec (cons 1 (cdr calc-curve-varnames))))))
165               ((and (>= key ?2) (<= key ?9))   ; polynomial               ((and (>= key ?2) (<= key ?9))   ; polynomial
166                (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))                (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
167                (setq model (math-build-polynomial-expr (cdr coefnames)                (setq calc-curve-model
168                                                        (nth 1 varnames))))                      (math-build-polynomial-expr (cdr calc-curve-coefnames)
169                                                    (nth 1 calc-curve-varnames))))
170               ((= key ?i)  ; exact polynomial               ((= key ?i)  ; exact polynomial
171                (calc-get-fit-variables 1 (1- (length (nth 1 data)))                (calc-get-fit-variables 1 (1- (length (nth 1 data)))
172                                        (and homog 0))                                        (and homog 0))
173                (setq model (math-build-polynomial-expr (cdr coefnames)                (setq calc-curve-model
174                                                        (nth 1 varnames))))                      (math-build-polynomial-expr (cdr calc-curve-coefnames)
175                                                    (nth 1 calc-curve-varnames))))
176               ((= key ?p)  ; power law               ((= key ?p)  ; power law
177                (calc-get-fit-variables nvars (1+ nvars) (and homog 1))                (calc-get-fit-variables calc-curve-nvars
178                (setq model (math-mul (nth 1 coefnames)                                        (1+ calc-curve-nvars) (and homog 1))
179                  (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
180                                      (calcFunc-reduce                                      (calcFunc-reduce
181                                       '(var mul var-mul)                                       '(var mul var-mul)
182                                       (calcFunc-map                                       (calcFunc-map
183                                        '(var pow var-pow)                                        '(var pow var-pow)
184                                        varnames                                        calc-curve-varnames
185                                        (cons 'vec (cdr (cdr coefnames))))))))                                        (cons 'vec (cdr (cdr calc-curve-coefnames))))))))
186               ((= key ?^)  ; exponential law               ((= key ?^)  ; exponential law
187                (calc-get-fit-variables nvars (1+ nvars) (and homog 1))                (calc-get-fit-variables calc-curve-nvars
188                (setq model (math-mul (nth 1 coefnames)                                        (1+ calc-curve-nvars) (and homog 1))
189                  (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
190                                      (calcFunc-reduce                                      (calcFunc-reduce
191                                       '(var mul var-mul)                                       '(var mul var-mul)
192                                       (calcFunc-map                                       (calcFunc-map
193                                        '(var pow var-pow)                                        '(var pow var-pow)
194                                        (cons 'vec (cdr (cdr coefnames)))                                        (cons 'vec (cdr (cdr calc-curve-coefnames)))
195                                        varnames)))))                                        calc-curve-varnames)))))
196               ((memq key '(?e ?E))               ((memq key '(?e ?E))
197                (calc-get-fit-variables nvars (1+ nvars) (and homog 1))                (calc-get-fit-variables calc-curve-nvars
198                (setq model (math-mul (nth 1 coefnames)                                        (1+ calc-curve-nvars) (and homog 1))
199                  (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
200                                      (calcFunc-reduce                                      (calcFunc-reduce
201                                       '(var mul var-mul)                                       '(var mul var-mul)
202                                       (calcFunc-map                                       (calcFunc-map
# Line 198  Line 207 
207                                            (^ 10 (var a var-a))))                                            (^ 10 (var a var-a))))
208                                        (calcFunc-map                                        (calcFunc-map
209                                         '(var mul var-mul)                                         '(var mul var-mul)
210                                         (cons 'vec (cdr (cdr coefnames)))                                         (cons 'vec (cdr (cdr calc-curve-coefnames)))
211                                         varnames))))))                                         calc-curve-varnames))))))
212               ((memq key '(?x ?X))               ((memq key '(?x ?X))
213                (calc-get-fit-variables nvars (1+ nvars) (and homog 0))                (calc-get-fit-variables calc-curve-nvars
214                (setq model (math-mul coefnames                                        (1+ calc-curve-nvars) (and homog 0))
215                                      (cons 'vec (cons 1 (cdr varnames)))))                (setq calc-curve-model (math-mul calc-curve-coefnames
216                (setq model (if (eq key ?x)                                      (cons 'vec (cons 1 (cdr calc-curve-varnames)))))
217                                (list 'calcFunc-exp model)                (setq calc-curve-model (if (eq key ?x)
218                              (list '^ 10 model))))                                (list 'calcFunc-exp calc-curve-model)
219                                (list '^ 10 calc-curve-model))))
220               ((memq key '(?l ?L))               ((memq key '(?l ?L))
221                (calc-get-fit-variables nvars (1+ nvars) (and homog 0))                (calc-get-fit-variables calc-curve-nvars
222                (setq model (math-mul coefnames                                        (1+ calc-curve-nvars) (and homog 0))
223                  (setq calc-curve-model (math-mul calc-curve-coefnames
224                                      (cons 'vec                                      (cons 'vec
225                                            (cons 1 (cdr (calcFunc-map                                            (cons 1 (cdr (calcFunc-map
226                                                          (if (eq key ?l)                                                          (if (eq key ?l)
227                                                              '(var ln var-ln)                                                              '(var ln var-ln)
228                                                            '(var log10                                                            '(var log10
229                                                                  var-log10))                                                                  var-log10))
230                                                          varnames)))))))                                                          calc-curve-varnames)))))))
231               ((= key ?q)               ((= key ?q)
232                (calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0))                (calc-get-fit-variables calc-curve-nvars
233                (let ((c coefnames)                                        (1+ (* 2 calc-curve-nvars)) (and homog 0))
234                      (v varnames))                (let ((c calc-curve-coefnames)
235                  (setq model (nth 1 c))                      (v calc-curve-varnames))
236                    (setq calc-curve-model (nth 1 c))
237                  (while (setq v (cdr v) c (cdr (cdr c)))                  (while (setq v (cdr v) c (cdr (cdr c)))
238                    (setq model (math-add                    (setq calc-curve-model (math-add
239                                 model                                 calc-curve-model
240                                 (list '*                                 (list '*
241                                       (car c)                                       (car c)
242                                       (list '^                                       (list '^
243                                             (list '- (car v) (nth 1 c))                                             (list '- (car v) (nth 1 c))
244                                             2)))))))                                             2)))))))
245               ((= key ?g)               ((= key ?g)
246                (setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")                (setq calc-curve-model
247                      varnames '(vec (var XFit var-XFit))                      (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
248                      coefnames '(vec (var AFit var-AFit)                      calc-curve-varnames '(vec (var XFit var-XFit))
249                        calc-curve-coefnames '(vec (var AFit var-AFit)
250                                      (var BFit var-BFit)                                      (var BFit var-BFit)
251                                      (var CFit var-CFit)))                                      (var CFit var-CFit)))
252                (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))                (calc-get-fit-variables 1 (1- (length calc-curve-coefnames))
253                                          (and homog 1)))
254               ((memq key '(?\$ ?\' ?u ?U))               ((memq key '(?\$ ?\' ?u ?U))
255                (let* ((defvars nil)                (let* ((defvars nil)
256                       (record-entry nil))                       (record-entry nil))
# Line 244  Line 258 
258                      (let* ((calc-dollar-values calc-arg-values)                      (let* ((calc-dollar-values calc-arg-values)
259                             (calc-dollar-used 0)                             (calc-dollar-used 0)
260                             (calc-hashes-used 0))                             (calc-hashes-used 0))
261                        (setq model (calc-do-alg-entry "" "Model formula: "))                        (setq calc-curve-model (calc-do-alg-entry "" "Model formula: "))
262                        (if (/= (length model) 1)                        (if (/= (length calc-curve-model) 1)
263                            (error "Bad format"))                            (error "Bad format"))
264                        (setq model (car model)                        (setq calc-curve-model (car calc-curve-model)
265                              record-entry t)                              record-entry t)
266                        (if (> calc-dollar-used 0)                        (if (> calc-dollar-used 0)
267                            (setq coefnames                            (setq calc-curve-coefnames
268                                  (cons 'vec                                  (cons 'vec
269                                        (nthcdr (- (length calc-arg-values)                                        (nthcdr (- (length calc-arg-values)
270                                                   calc-dollar-used)                                                   calc-dollar-used)
271                                                (reverse calc-arg-values))))                                                (reverse calc-arg-values))))
272                          (if (> calc-hashes-used 0)                          (if (> calc-hashes-used 0)
273                              (setq coefnames                              (setq calc-curve-coefnames
274                                    (cons 'vec (calc-invent-args                                    (cons 'vec (calc-invent-args
275                                                calc-hashes-used))))))                                                calc-hashes-used))))))
276                    (progn                    (progn
277                      (setq model (cond ((eq key ?u)                      (setq calc-curve-model (cond ((eq key ?u)
278                                         (calc-var-value 'var-Model1))                                         (calc-var-value 'var-Model1))
279                                        ((eq key ?U)                                        ((eq key ?U)
280                                         (calc-var-value 'var-Model2))                                         (calc-var-value 'var-Model2))
281                                        (t (calc-top 1))))                                        (t (calc-top 1))))
282                      (or model (error "User model not yet defined"))                      (or calc-curve-model (error "User model not yet defined"))
283                      (if (math-vectorp model)                      (if (math-vectorp calc-curve-model)
284                          (if (and (memq (length model) '(3 4))                          (if (and (memq (length calc-curve-model) '(3 4))
285                                   (not (math-objvecp (nth 1 model)))                                   (not (math-objvecp (nth 1 calc-curve-model)))
286                                   (math-vectorp (nth 2 model))                                   (math-vectorp (nth 2 calc-curve-model))
287                                   (or (null (nth 3 model))                                   (or (null (nth 3 calc-curve-model))
288                                       (math-vectorp (nth 3 model))))                                       (math-vectorp (nth 3 calc-curve-model))))
289                              (setq varnames (nth 2 model)                              (setq calc-curve-varnames (nth 2 calc-curve-model)
290                                    coefnames (or (nth 3 model)                                    calc-curve-coefnames
291                                                  (cons 'vec                                    (or (nth 3 calc-curve-model)
292                                                        (math-all-vars-but                                        (cons 'vec
293                                                         model varnames)))                                              (math-all-vars-but
294                                    model (nth 1 model))                                               calc-curve-model calc-curve-varnames)))
295                                      calc-curve-model (nth 1 calc-curve-model))
296                            (error "Incorrect model specifier")))))                            (error "Incorrect model specifier")))))
297                  (or varnames                  (or calc-curve-varnames
298                      (let ((with-y (eq (car-safe model) 'calcFunc-eq)))                      (let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq)))
299                        (if coefnames                        (if calc-curve-coefnames
300                            (calc-get-fit-variables (if with-y (1+ nvars) nvars)                            (calc-get-fit-variables
301                                                    (1- (length coefnames))                             (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
302                                                    (math-all-vars-but                             (1- (length calc-curve-coefnames))
303                                                     model coefnames)                             (math-all-vars-but
304                                                    nil with-y)                              calc-curve-model calc-curve-coefnames)
305                          (let* ((coefs (math-all-vars-but model nil))                             nil with-y)
306                            (let* ((coefs (math-all-vars-but calc-curve-model nil))
307                                 (vars nil)                                 (vars nil)
308                                 (n (- (length coefs) nvars (if with-y 2 1)))                                 (n (- (length coefs) calc-curve-nvars (if with-y 2 1)))
309                                 p)                                 p)
310                            (if (< n 0)                            (if (< n 0)
311                                (error "Not enough variables in model"))                                (error "Not enough variables in model"))
312                            (setq p (nthcdr n coefs))                            (setq p (nthcdr n coefs))
313                            (setq vars (cdr p))                            (setq vars (cdr p))
314                            (setcdr p nil)                            (setcdr p nil)
315                            (calc-get-fit-variables (if with-y (1+ nvars) nvars)                            (calc-get-fit-variables
316                                                    (length coefs)                             (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
317                                                    vars coefs with-y)))))                             (length coefs)
318                               vars coefs with-y)))))
319                  (if record-entry                  (if record-entry
320                      (calc-record (list 'vec model varnames coefnames)                      (calc-record (list 'vec calc-curve-model
321                                           calc-curve-varnames calc-curve-coefnames)
322                                   "modl"))))                                   "modl"))))
323               (t (beep))))               (t (beep))))
324       (let ((calc-fit-to-trail t))       (let ((calc-fit-to-trail t))
325         (calc-enter-result n (substring (symbol-name func) 9)         (calc-enter-result n (substring (symbol-name func) 9)
326                            (list func model                            (list func calc-curve-model
327                                  (if (= (length varnames) 2)                                  (if (= (length calc-curve-varnames) 2)
328                                      (nth 1 varnames)                                      (nth 1 calc-curve-varnames)
329                                    varnames)                                    calc-curve-varnames)
330                                  (if (= (length coefnames) 2)                                  (if (= (length calc-curve-coefnames) 2)
331                                      (nth 1 coefnames)                                      (nth 1 calc-curve-coefnames)
332                                    coefnames)                                    calc-curve-coefnames)
333                                  data))                                  data))
334         (if (consp calc-fit-to-trail)         (if (consp calc-fit-to-trail)
335             (calc-record (calc-normalize calc-fit-to-trail) "parm"))))))             (calc-record (calc-normalize calc-fit-to-trail) "parm"))))))
# Line 340  Line 358 
358        (calc-invent-variables num but t base))))        (calc-invent-variables num but t base))))
359    
360  (defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)  (defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
361    (or (= nv (if with-y (1+ nvars) nvars))    (or (= nv (if with-y (1+ calc-curve-nvars) calc-curve-nvars))
362        (error "Wrong number of data vectors for this type of model"))        (error "Wrong number of data vectors for this type of model"))
363    (if (integerp defv)    (if (integerp defv)
364        (setq homog defv        (setq homog defv
# Line 388  Line 406 
406          (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))          (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
407      (if homog      (if homog
408          (setq coefs (cons 'vec (cons homog (cdr coefs)))))          (setq coefs (cons 'vec (cons homog (cdr coefs)))))
409      (if varnames      (if calc-curve-varnames
410          (setq model (math-multi-subst model (cdr varnames) (cdr vars))))          (setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-varnames) (cdr vars))))
411      (if coefnames      (if calc-curve-coefnames
412          (setq model (math-multi-subst model (cdr coefnames) (cdr coefs))))          (setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-coefnames) (cdr coefs))))
413      (setq varnames vars      (setq calc-curve-varnames vars
414            coefnames coefs)))            calc-curve-coefnames coefs)))
415    
416    
417    
# Line 401  Line 419 
419  ;;; The following algorithms are from Numerical Recipes chapter 9.  ;;; The following algorithms are from Numerical Recipes chapter 9.
420    
421  ;;; "rtnewt" with safety kludges  ;;; "rtnewt" with safety kludges
422    
423    (defvar var-DUMMY)
424    
425  (defun math-newton-root (expr deriv guess orig-guess limit)  (defun math-newton-root (expr deriv guess orig-guess limit)
426    (math-working "newton" guess)    (math-working "newton" guess)
427    (let* ((var-DUMMY guess)    (let* ((var-DUMMY guess)
# Line 494  Line 515 
515                                     low vlow high vhigh))))))                                     low vlow high vhigh))))))
516    
517  ;;; Search for a root in an interval with no overt zero crossing.  ;;; Search for a root in an interval with no overt zero crossing.
518    
519    ;; The variable math-root-widen is local to math-find-root, but
520    ;; is used by math-search-root, which is called (directly and
521    ;; indirectly) by math-find-root.
522    (defvar math-root-widen)
523    
524  (defun math-search-root (expr deriv low vlow high vhigh)  (defun math-search-root (expr deriv low vlow high vhigh)
525    (let (found)    (let (found)
526      (if root-widen      (if math-root-widen
527          (let ((iters 0)          (let ((iters 0)
528                (iterlim (if (eq root-widen 'point)                (iterlim (if (eq math-root-widen 'point)
529                             (+ calc-internal-prec 10)                             (+ calc-internal-prec 10)
530                           20))                           20))
531                (factor (if (eq root-widen 'point)                (factor (if (eq math-root-widen 'point)
532                            '(float 9 0)                            '(float 9 0)
533                          '(float 16 -1)))                          '(float 16 -1)))
534                (prev nil) vprev waslow                (prev nil) vprev waslow
# Line 600  Line 627 
627      (list 'vec mid vmid)))      (list 'vec mid vmid)))
628    
629  ;;; "mnewt"  ;;; "mnewt"
630    
631    (defvar math-root-vars [(var DUMMY var-DUMMY)])
632    
633  (defun math-newton-multi (expr jacob n guess orig-guess limit)  (defun math-newton-multi (expr jacob n guess orig-guess limit)
634    (let ((m -1)    (let ((m -1)
635          (p guess)          (p guess)
# Line 624  Line 654 
654            (math-reject-arg nil "*Newton's method failed to converge"))            (math-reject-arg nil "*Newton's method failed to converge"))
655        (list 'vec next expr-val))))        (list 'vec next expr-val))))
656    
 (defvar math-root-vars [(var DUMMY var-DUMMY)])  
657    
658  (defun math-find-root (expr var guess root-widen)  (defun math-find-root (expr var guess math-root-widen)
659    (if (eq (car-safe expr) 'vec)    (if (eq (car-safe expr) 'vec)
660        (let ((n (1- (length expr)))        (let ((n (1- (length expr)))
661              (calc-symbolic-mode nil)              (calc-symbolic-mode nil)
# Line 710  Line 739 
739                      var-DUMMY guess                      var-DUMMY guess
740                      vlow (math-evaluate-expr expr)                      vlow (math-evaluate-expr expr)
741                      vhigh vlow                      vhigh vlow
742                      root-widen 'point)                      math-root-widen 'point)
743              (if (eq (car guess) 'intv)              (if (eq (car guess) 'intv)
744                  (progn                  (progn
745                    (or (math-constp guess) (math-reject-arg guess 'constp))                    (or (math-constp guess) (math-reject-arg guess 'constp))
# Line 752  Line 781 
781    
782  ;;; The following algorithms come from Numerical Recipes, chapter 10.  ;;; The following algorithms come from Numerical Recipes, chapter 10.
783    
784    (defvar math-min-vars [(var DUMMY var-DUMMY)])
785    
786  (defun math-min-eval (expr a)  (defun math-min-eval (expr a)
787    (if (Math-vectorp a)    (if (Math-vectorp a)
788        (let ((m -1))        (let ((m -1))
# Line 894  Line 925 
925          (tol (list 'float 1 (- -1 prec)))          (tol (list 'float 1 (- -1 prec)))
926          (zeps (list 'float 1 (- -5 prec)))          (zeps (list 'float 1 (- -5 prec)))
927          (e '(float 0 0))          (e '(float 0 0))
928          u vu xm tol1 tol2 etemp p q r xv xw)          d u vu xm tol1 tol2 etemp p q r xv xw)
929      (while (progn      (while (progn
930               (setq xm (math-mul-float '(float 5 -1)               (setq xm (math-mul-float '(float 5 -1)
931                                        (math-add-float a b))                                        (math-add-float a b))
# Line 1056  Line 1087 
1087      (list (math-add line-p xi) xi (nth 2 res))))      (list (math-add line-p xi) xi (nth 2 res))))
1088    
1089    
 (defvar math-min-vars [(var DUMMY var-DUMMY)])  
   
1090  (defun math-find-minimum (expr var guess min-widen)  (defun math-find-minimum (expr var guess min-widen)
1091    (let* ((calc-symbolic-mode nil)    (let* ((calc-symbolic-mode nil)
1092           (n 0)           (n 0)
# Line 1072  Line 1101 
1101          (math-dimension-error))          (math-dimension-error))
1102      (while (setq var (cdr var) guess (cdr guess))      (while (setq var (cdr var) guess (cdr guess))
1103        (or (eq (car-safe (car var)) 'var)        (or (eq (car-safe (car var)) 'var)
1104            (math-reject-arg (car vg) "*Expected a variable"))            (math-reject-arg (car var) "*Expected a variable"))
1105        (or (math-expr-contains expr (car var))        (or (math-expr-contains expr (car var))
1106            (math-reject-arg (car var)            (math-reject-arg (car var)
1107                             "*Formula does not contain specified variable"))                             "*Formula does not contain specified variable"))
# Line 1314  Line 1343 
1343    
1344    
1345  ;;; Open Romberg method; "qromo" in section 4.4.  ;;; Open Romberg method; "qromo" in section 4.4.
1346    
1347    ;; The variable math-ninteg-temp is local to math-ninteg-romberg,
1348    ;; but is used by math-ninteg-midpoint, which is used by
1349    ;; math-ninteg-romberg.
1350    (defvar math-ninteg-temp)
1351    
1352  (defun math-ninteg-romberg (func expr lo hi mode)  (defun math-ninteg-romberg (func expr lo hi mode)
1353    (let ((curh '(float 1 0))    (let ((curh '(float 1 0))
1354          (h nil)          (h nil)
# Line 1321  Line 1356 
1356          (j 0)          (j 0)
1357          (ss nil)          (ss nil)
1358          (prec calc-internal-prec)          (prec calc-internal-prec)
1359          (integ-temp nil))          (math-ninteg-temp nil))
1360      (math-with-extra-prec 2      (math-with-extra-prec 2
1361        ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.        ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.
1362        (or (while (and (null ss) (<= (setq j (1+ j)) 8))        (or (while (and (null ss) (<= (setq j (1+ j)) 8))
# Line 1332  Line 1367 
1367                    (if (math-lessp (math-abs (nth 1 res))                    (if (math-lessp (math-abs (nth 1 res))
1368                                    (calcFunc-scf (math-abs (car res))                                    (calcFunc-scf (math-abs (car res))
1369                                                  (- prec)))                                                  (- prec)))
1370                        (setq math-ninteg-convergence j                        (setq ss (car res)))))
                             ss (car res)))))  
1371              (if (>= j 5)              (if (>= j 5)
1372                  (setq s (cdr s)                  (setq s (cdr s)
1373                        h (cdr h)))                        h (cdr h)))
# Line 1354  Line 1388 
1388      res))      res))
1389    
1390    
1391  (defun math-ninteg-midpoint (expr lo hi mode)    ; uses "integ-temp"  (defun math-ninteg-midpoint (expr lo hi mode)    ; uses "math-ninteg-temp"
1392    (if (eq mode 'inf)    (if (eq mode 'inf)
1393        (let ((math-infinite-mode t) temp)        (let ((math-infinite-mode t) temp)
1394          (setq temp (math-div 1 lo)          (setq temp (math-div 1 lo)
1395                lo (math-div 1 hi)                lo (math-div 1 hi)
1396                hi temp)))                hi temp)))
1397    (if integ-temp    (if math-ninteg-temp
1398        (let* ((it3 (* 3 (car integ-temp)))        (let* ((it3 (* 3 (car math-ninteg-temp)))
1399               (math-working-step-2 (* 2 (car integ-temp)))               (math-working-step-2 (* 2 (car math-ninteg-temp)))
1400               (math-working-step 0)               (math-working-step 0)
1401               (range (math-sub hi lo))               (range (math-sub hi lo))
1402               (del (math-div range (math-float it3)))               (del (math-div range (math-float it3)))
# Line 1371  Line 1405 
1405               (x (math-add lo (math-mul '(float 5 -1) del)))               (x (math-add lo (math-mul '(float 5 -1) del)))
1406               (sum '(float 0 0))               (sum '(float 0 0))
1407               (j 0) temp)               (j 0) temp)
1408          (while (<= (setq j (1+ j)) (car integ-temp))          (while (<= (setq j (1+ j)) (car math-ninteg-temp))
1409            (setq math-working-step (1+ math-working-step)            (setq math-working-step (1+ math-working-step)
1410                  temp (math-ninteg-evaluate expr x mode)                  temp (math-ninteg-evaluate expr x mode)
1411                  math-working-step (1+ math-working-step)                  math-working-step (1+ math-working-step)
# Line 1379  Line 1413 
1413                                                    expr (math-add x del2)                                                    expr (math-add x del2)
1414                                                    mode)))                                                    mode)))
1415                  x (math-add x del3)))                  x (math-add x del3)))
1416          (setq integ-temp (list it3          (setq math-ninteg-temp (list it3
1417                                 (math-add (math-div (nth 1 integ-temp)                                       (math-add (math-div (nth 1 math-ninteg-temp)
1418                                                     '(float 3 0))                                                           '(float 3 0))
1419                                           (math-mul sum del)))))                                                 (math-mul sum del)))))
1420      (setq integ-temp (list 1 (math-mul      (setq math-ninteg-temp (list 1 (math-mul
1421                                (math-sub hi lo)                                      (math-sub hi lo)
1422                                (math-ninteg-evaluate                                      (math-ninteg-evaluate
1423                                 expr                                       expr
1424                                 (math-mul (math-add lo hi) '(float 5 -1))                                       (math-mul (math-add lo hi) '(float 5 -1))
1425                                 mode)))))                                       mode)))))
1426    (nth 1 integ-temp))    (nth 1 math-ninteg-temp))
1427    
1428    
1429    
# Line 1427  Line 1461 
1461      (math-with-extra-prec 2      (math-with-extra-prec 2
1462        (math-general-fit expr vars coefs data 'full))))        (math-general-fit expr vars coefs data 'full))))
1463    
1464    ;; The variables math-fit-first-var, math-fit-first-coef and
1465    ;; math-fit-new-coefs are local to math-general-fit, but are used by
1466    ;; calcFunc-fitvar, calcFunc-fitparam and calcFunc-fitdummy
1467    ;; (respectively), which are used by math-general-fit.
1468    (defvar math-fit-first-var)
1469    (defvar math-fit-first-coef)
1470    (defvar math-fit-new-coefs)
1471    
1472  (defun math-general-fit (expr vars coefs data mode)  (defun math-general-fit (expr vars coefs data mode)
1473    (let ((calc-simplify-mode nil)    (let ((calc-simplify-mode nil)
1474          (math-dummy-counter math-dummy-counter)          (math-dummy-counter math-dummy-counter)
1475          (math-in-fit 1)          (math-in-fit 1)
1476          (extended (eq mode 'full))          (extended (eq mode 'full))
1477          (first-coef math-dummy-counter)          (math-fit-first-coef math-dummy-counter)
1478          first-var          math-fit-first-var
1479          (plain-expr expr)          (plain-expr expr)
1480          orig-expr          orig-expr
1481          have-sdevs need-chisq chisq          have-sdevs need-chisq chisq
# Line 1441  Line 1483 
1483          (y-filter nil)          (y-filter nil)
1484          y-dummy          y-dummy
1485          (coef-filters nil)          (coef-filters nil)
1486          new-coefs          math-fit-new-coefs
1487          (xy-values nil)          (xy-values nil)
1488          (weights nil)          (weights nil)
1489          (var-YVAL nil) (var-YVALX nil)          (var-YVAL nil) (var-YVALX nil)
# Line 1496  Line 1538 
1538        (setq dummy (math-dummy-variable)        (setq dummy (math-dummy-variable)
1539              expr (math-expr-subst expr (car p)              expr (math-expr-subst expr (car p)
1540                                    (list 'calcFunc-fitparam                                    (list 'calcFunc-fitparam
1541                                          (- math-dummy-counter first-coef)))))                                          (- math-dummy-counter math-fit-first-coef)))))
1542      (setq first-var math-dummy-counter      (setq math-fit-first-var math-dummy-counter
1543            p vars)            p vars)
1544      (while (setq p (cdr p))      (while (setq p (cdr p))
1545        (or (eq (car-safe (car p)) 'var)        (or (eq (car-safe (car p)) 'var)
# Line 1505  Line 1547 
1547        (setq dummy (math-dummy-variable)        (setq dummy (math-dummy-variable)
1548              expr (math-expr-subst expr (car p)              expr (math-expr-subst expr (car p)
1549                                    (list 'calcFunc-fitvar                                    (list 'calcFunc-fitvar
1550                                          (- math-dummy-counter first-var)))))                                          (- math-dummy-counter math-fit-first-var)))))
1551      (if (< math-dummy-counter (+ first-var v))      (if (< math-dummy-counter (+ math-fit-first-var v))
1552          (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed          (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
1553      (setq y-dummy dummy      (setq y-dummy dummy
1554            orig-expr expr)            orig-expr expr)
# Line 1565  Line 1607 
1607                    (setq sigmasqr (math-add (math-sqr (nth 2 xval))                    (setq sigmasqr (math-add (math-sqr (nth 2 xval))
1608                                             (or sigmasqr 0))                                             (or sigmasqr 0))
1609                          xval (nth 1 xval))))                          xval (nth 1 xval))))
1610              (set (nth 2 (aref math-dummy-vars (+ first-var j))) xval)              (set (nth 2 (aref math-dummy-vars (+ math-fit-first-var j))) xval)
1611              (setq j (1+ j)))              (setq j (1+ j)))
1612    
1613            ;; Compute Y value for this data point.            ;; Compute Y value for this data point.
# Line 1656  Line 1698 
1698                      xy-values (cdr xy-values)))))                      xy-values (cdr xy-values)))))
1699    
1700        ;; Convert coefficients back into original terms.        ;; Convert coefficients back into original terms.
1701        (setq new-coefs (copy-sequence beta))        (setq math-fit-new-coefs (copy-sequence beta))
1702        (let* ((bp new-coefs)        (let* ((bp math-fit-new-coefs)
1703               (cp covar)               (cp covar)
1704               (sigdat 1)               (sigdat 1)
1705               (math-in-fit 3)               (math-in-fit 3)
# Line 1673  Line 1715 
1715                            (math-sqrt (math-mul (nth (setq j (1+ j))                            (math-sqrt (math-mul (nth (setq j (1+ j))
1716                                                      (car (setq cp (cdr cp))))                                                      (car (setq cp (cdr cp))))
1717                                                 sigdat))))))                                                 sigdat))))))
1718          (setq new-coefs (math-evaluate-expr coef-filters))          (setq math-fit-new-coefs (math-evaluate-expr coef-filters))
1719          (if calc-fit-to-trail          (if calc-fit-to-trail
1720              (let ((bp new-coefs)              (let ((bp math-fit-new-coefs)
1721                    (cp coefs)                    (cp coefs)
1722                    (vec nil))                    (vec nil))
1723                (while (setq bp (cdr bp) cp (cdr cp))                (while (setq bp (cdr bp) cp (cdr cp))
# Line 1695  Line 1737 
1737                      (setq vec (cons (list 'calcFunc-fitparam n) vec)                      (setq vec (cons (list 'calcFunc-fitparam n) vec)
1738                            n (1- n)))                            n (1- n)))
1739                    vec)                    vec)
1740                  (append (cdr new-coefs) (cdr vars))))                  (append (cdr math-fit-new-coefs) (cdr vars))))
1741    
1742      ;; Package the result.      ;; Package the result.
1743      (math-normalize      (math-normalize
# Line 1719  Line 1761 
1761  (defun calcFunc-fitvar (x)  (defun calcFunc-fitvar (x)
1762    (if (>= math-in-fit 2)    (if (>= math-in-fit 2)
1763        (progn        (progn
1764          (setq x (aref math-dummy-vars (+ first-var x -1)))          (setq x (aref math-dummy-vars (+ math-fit-first-var x -1)))
1765          (or (calc-var-value (nth 2 x)) x))          (or (calc-var-value (nth 2 x)) x))
1766      (math-reject-arg x)))      (math-reject-arg x)))
1767    
1768  (defun calcFunc-fitparam (x)  (defun calcFunc-fitparam (x)
1769    (if (>= math-in-fit 2)    (if (>= math-in-fit 2)
1770        (progn        (progn
1771          (setq x (aref math-dummy-vars (+ first-coef x -1)))          (setq x (aref math-dummy-vars (+ math-fit-first-coef x -1)))
1772          (or (calc-var-value (nth 2 x)) x))          (or (calc-var-value (nth 2 x)) x))
1773      (math-reject-arg x)))      (math-reject-arg x)))
1774    
1775  (defun calcFunc-fitdummy (x)  (defun calcFunc-fitdummy (x)
1776    (if (= math-in-fit 3)    (if (= math-in-fit 3)
1777        (nth x new-coefs)        (nth x math-fit-new-coefs)
1778      (math-reject-arg x)))      (math-reject-arg x)))
1779    
1780  (defun calcFunc-hasfitvars (expr)  (defun calcFunc-hasfitvars (expr)
# Line 1759  Line 1801 
1801      (sort (mapcar 'car vars)      (sort (mapcar 'car vars)
1802            (function (lambda (x y) (string< (nth 1 x) (nth 1 y)))))))            (function (lambda (x y) (string< (nth 1 x) (nth 1 y)))))))
1803    
1804    ;; The variables math-all-vars-vars (the vars for math-all-vars) and
1805    ;; math-all-vars-found are local to math-all-vars-in, but are used by
1806    ;; math-all-vars-rec which is called by math-all-vars-in.
1807    (defvar math-all-vars-vars)
1808    (defvar math-all-vars-found)
1809    
1810  (defun math-all-vars-in (expr)  (defun math-all-vars-in (expr)
1811    (let ((vars nil)    (let ((math-all-vars-vars nil)
1812          found)          math-all-vars-found)
1813      (math-all-vars-rec expr)      (math-all-vars-rec expr)
1814      vars))      math-all-vars-vars))
1815    
1816  (defun math-all-vars-rec (expr)  (defun math-all-vars-rec (expr)
1817    (if (Math-primp expr)    (if (Math-primp expr)
1818        (if (eq (car-safe expr) 'var)        (if (eq (car-safe expr) 'var)
1819            (or (math-const-var expr)            (or (math-const-var expr)
1820                (if (setq found (assoc expr vars))                (if (setq math-all-vars-found (assoc expr math-all-vars-vars))
1821                    (setcdr found (1+ (cdr found)))                    (setcdr math-all-vars-found (1+ (cdr math-all-vars-found)))
1822                  (setq vars (cons (cons expr 1) vars)))))                  (setq math-all-vars-vars (cons (cons expr 1) math-all-vars-vars)))))
1823      (while (setq expr (cdr expr))      (while (setq expr (cdr expr))
1824        (math-all-vars-rec (car expr)))))        (math-all-vars-rec (car expr)))))
1825    
1826    (provide 'calcalg3)
1827    
1828  ;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6  ;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6
1829  ;;; calcalg3.el ends here  ;;; calcalg3.el ends here

Legend:
Removed from v.1.4.4.2  
changed lines
  Added in v.1.4.4.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26