From: Douglas <dt...@us...> - 2011-11-10 12:07:46
|
This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "Maxima, A Computer Algebra System". The branch, master has been updated via 2c908019e696320c206bacaacc689da19fdf92a5 (commit) from 2eb7099928f9e9d7e48040fe0b883ec7f934bb3d (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 2c908019e696320c206bacaacc689da19fdf92a5 Author: Douglas Crosher <dt...@us...> Date: Thu Nov 10 22:53:04 2011 +1100 o Better support for compiling with the 'flonum type as 'long-float or 'double-double on CMUCL. Change float constants using the #\d exponent tag to use the #\e tag so they are read in the default float format which will be the 'flonum type. Change references to 'double-float to 'flonum. diff --git a/share/hypergeometric/nfloat.lisp b/share/hypergeometric/nfloat.lisp index 487c028..1b0ffe2 100644 --- a/share/hypergeometric/nfloat.lisp +++ b/share/hypergeometric/nfloat.lisp @@ -213,7 +213,7 @@ ;; (list (bigfloat::to e) 0)) ((maxima::complex-number-p e #'(lambda (s) (or (maxima::$ratnump s) (maxima::$numberp s)))) - (setq e (bigfloat::to (if (> bits #.(float-digits 1.0d0)) (maxima::$bfloat e) (maxima::$float e)))) + (setq e (bigfloat::to (if (> bits #.(float-digits 1.0e0)) (maxima::$bfloat e) (maxima::$float e)))) (list e (abs e))) ((and (atom e) (maxima::mget e '$numer)) diff --git a/share/lapack/dgemm.lisp b/share/lapack/dgemm.lisp index a6b67da..c6d4b4c 100644 --- a/share/lapack/dgemm.lisp +++ b/share/lapack/dgemm.lisp @@ -2,7 +2,7 @@ ;;; multiplication. (in-package :maxima) -(defun %%dgemm (a b &key (c nil cp) transpose_a transpose_b (alpha 1d0) (beta 0d0 betap)) +(defun %%dgemm (a b &key (c nil cp) transpose_a transpose_b (alpha 1e0) (beta 0e0 betap)) (flet ((maybe-transpose-dims (transp row col) (if transp (values col row) @@ -26,7 +26,7 @@ ;; take this to mean that we just want to add C, ;; implying that beta = 1. (when (and cp (not betap)) - (setf beta 0d0)) + (setf beta 0e0)) ;; Now for some error checking. This is a bit redundant ;; since dgemm does some error checking, but I think we ;; prefer not to see messages from dgemm. It's better @@ -57,8 +57,8 @@ ;; Force beta to be zero to tell LAPACK ;; not to add C. But we still need to ;; create a matrix. - (setf beta 0d0) - (make-array (* c-nrows c-ncols) :element-type 'double-float)))) + (setf beta 0e0) + (make-array (* c-nrows c-ncols) :element-type 'flonum)))) (trans-a (if transpose_a "t" "n")) (trans-b (if transpose_b "t" "n"))) (blas::dgemm trans-a trans-b diff --git a/share/lapack/dgeqrf.lisp b/share/lapack/dgeqrf.lisp index f773f6b..d3be4f5 100644 --- a/share/lapack/dgeqrf.lisp +++ b/share/lapack/dgeqrf.lisp @@ -47,7 +47,7 @@ ((r-nrow a-nrow) (r-ncol a-ncol) (r-mat (make-array (* r-ncol r-nrow) :element-type 'flonum - :initial-element (coerce 0 'flonum)))) + :initial-element 0e0))) (dotimes (j a-ncol) (dotimes (i (1+ j)) (if (< i r-nrow) @@ -59,9 +59,9 @@ (defun dgeqrf-unpack-q (a-nrow a-ncol a-mat tau) (let ((q-mat (make-array (* a-nrow a-nrow) :element-type 'flonum - :initial-element (coerce 0 'flonum)))) + :initial-element 0e0))) (dotimes (i a-nrow) - (setf (aref q-mat (+ (* i a-nrow) i)) (coerce 1 'flonum))) + (setf (aref q-mat (+ (* i a-nrow) i)) 1e0)) (dotimes (i (min a-nrow a-ncol)) (let ((h-mat-i (dgeqrf-h i tau a-nrow a-mat))) (dgeqrf-multiply-into a-nrow q-mat h-mat-i))) @@ -86,7 +86,7 @@ (t 0))) (defun dgeqrf-multiply-into (n mat-1 mat-2) - (let ((row (make-array n :element-type 'flonum :initial-element (coerce 0 'flonum)))) + (let ((row (make-array n :element-type 'flonum :initial-element 0e0))) (dotimes (i n) (dotimes (j n) (setf (aref row j) (dgeqrf-inner-product n mat-1 mat-2 i j))) diff --git a/share/lbfgs/maxima-lbfgs.lisp b/share/lbfgs/maxima-lbfgs.lisp index 4cfd2e8..5fec58c 100644 --- a/share/lbfgs/maxima-lbfgs.lisp +++ b/share/lbfgs/maxima-lbfgs.lisp @@ -84,19 +84,19 @@ estimates : lbfgs ('[F(a, b, c), [F1(a, b, c), F2(a, b, c), F3(a, b, c)]], (m $lbfgs_ncorrections) (nwork (+ (* n (+ (* 2 m) 1)) (* 2 m))) - (xtol double-float-epsilon) + (xtol flonum-epsilon) (iflag 0) - (scache (make-array n :element-type 'double-float)) - (w (make-array nwork :element-type 'double-float)) - (diag (make-array n :element-type 'double-float)) - (g (make-array n :element-type 'double-float)) - (x (make-array n :element-type 'double-float)) + (scache (make-array n :element-type 'flonum)) + (w (make-array nwork :element-type 'flonum)) + (diag (make-array n :element-type 'flonum)) + (g (make-array n :element-type 'flonum)) + (x (make-array n :element-type 'flonum)) (iprint (make-array 2 :element-type 'f2cl-lib:integer4)) (diagco nil) (return-value '((mlist))) - (f 0.0d0) + (f 0.0) FOM-function FOM-grad-lambda FOM-grad-expr FOM-grad-function) diff --git a/src/bessel.lisp b/src/bessel.lisp index 137f305..0c53b9b 100644 --- a/src/bessel.lisp +++ b/src/bessel.lisp @@ -311,7 +311,7 @@ ;; We have a real arg and order > 0 and order not 0 or 1 ;; for this case we can call the function dbesj (multiple-value-bind (n alpha) (floor (float order)) - (let ((jvals (make-array (1+ n) :element-type 'double-float))) + (let ((jvals (make-array (1+ n) :element-type 'flonum))) (slatec:dbesj (abs (float arg)) alpha (1+ n) jvals 0) (cond ((>= arg 0) (aref jvals n)) @@ -347,8 +347,8 @@ (* 0.5 (+ (hankel-1 order arg) (hankel-2 order arg)))) (t (multiple-value-bind (n alpha) (floor (float order)) - (let ((cyr (make-array (1+ n) :element-type 'double-float)) - (cyi (make-array (1+ n) :element-type 'double-float))) + (let ((cyr (make-array (1+ n) :element-type 'flonum)) + (cyi (make-array (1+ n) :element-type 'flonum))) (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz v-ierr) (slatec:zbesj (float (realpart arg)) @@ -526,7 +526,7 @@ ($float t) (order ($float order)) (arg ($float arg)) - (dpi (coerce pi 'double-float))) + (dpi (coerce pi 'flonum))) ($float ($rectform (add @@ -665,7 +665,7 @@ (t result)))))) (t (multiple-value-bind (n alpha) (floor (float order)) - (let ((jvals (make-array (1+ n) :element-type 'double-float))) + (let ((jvals (make-array (1+ n) :element-type 'flonum))) ;; First we do the calculation for an positive argument. (slatec:dbesy (abs (float arg)) alpha (1+ n) jvals) @@ -673,7 +673,7 @@ (cond ((>= arg 0) (aref jvals n)) (t - (let* ((dpi (coerce pi 'double-float)) + (let* ((dpi (coerce pi 'flonum)) (s1 (cis (- (* order dpi)))) (s2 (* #c(0 2) (cos (* order dpi))))) (let ((result (+ (* s1 (aref jvals n)) @@ -700,10 +700,10 @@ (complex 0 2))) (t (multiple-value-bind (n alpha) (floor (float order)) - (let ((cyr (make-array (1+ n) :element-type 'double-float)) - (cyi (make-array (1+ n) :element-type 'double-float)) - (cwrkr (make-array (1+ n) :element-type 'double-float)) - (cwrki (make-array (1+ n) :element-type 'double-float))) + (let ((cyr (make-array (1+ n) :element-type 'flonum)) + (cyi (make-array (1+ n) :element-type 'flonum)) + (cwrkr (make-array (1+ n) :element-type 'flonum)) + (cwrki (make-array (1+ n) :element-type 'flonum))) (multiple-value-bind (v-zr v-zi v-fnu v-kode v-n v-cyr v-cyi v-nz v-cwrkr v-cwrki v-ierr) (slatec::zbesy (float (realpart arg)) @@ -971,7 +971,7 @@ (t ;; Now the case order > 0 and arg >= 0 (multiple-value-bind (n alpha) (floor (float order)) - (let ((jvals (make-array (1+ n) :element-type 'double-float))) + (let ((jvals (make-array (1+ n) :element-type 'flonum))) (slatec:dbesi (float (realpart arg)) alpha 1 (1+ n) jvals 0) (aref jvals n))))))) (t @@ -995,7 +995,7 @@ (cond ((minusp order) ;; I(-a,z) = I(a,z) + (2/pi)*sin(pi*a)*K(a,z) (+ (complex (aref cyr n) (aref cyi n)) - (let ((dpi (coerce pi 'double-float))) + (let ((dpi (coerce pi 'flonum))) (* (/ 2.0 dpi) (sin (* dpi (- order))) (bessel-k (- order) arg))))) @@ -1254,7 +1254,7 @@ ;; This is the extension for negative arg. ;; We use the following formula for evaluation: ;; K[v](-z) = exp(-i*pi*v) * K[n][z]-i * pi *I[n](z) - (let* ((dpi (coerce pi 'double-float)) + (let* ((dpi (coerce pi 'flonum)) (s1 (cis (* dpi (- (abs order))))) (s2 (* (complex 0 -1) dpi)) (result (+ (* s1 (bessel-k (abs order) (- arg))) @@ -1278,7 +1278,7 @@ ;; From A&S 9.6.6, K(-v,z) = K(v,z), so take the ;; absolute value of the order. (multiple-value-bind (n alpha) (floor (abs (float order))) - (let ((jvals (make-array (1+ n) :element-type 'double-float))) + (let ((jvals (make-array (1+ n) :element-type 'flonum))) (slatec:dbesk (float arg) alpha 1 (1+ n) jvals 0) (aref jvals n))))))) (t diff --git a/src/clmacs.lisp b/src/clmacs.lisp index 39d14c2..858a8f1 100644 --- a/src/clmacs.lisp +++ b/src/clmacs.lisp @@ -334,7 +334,7 @@ (setf ext::*intexp-maximum-exponent* 100000) ;;;; Setup the mapping from the Maxima 'flonum float type to a CL float type. ;;;; -;;;; Add :flonum-log to *features* if you want flonum to be a +;;;; Add :flonum-long to *features* if you want flonum to be a ;;;; long-float. Or add :flonum-double-double if you want flonum to ;;;; be a double-double (currently only for CMUCL). Otherwise, you ;;;; get double-float as the flonum type. @@ -377,7 +377,7 @@ #+flonum-long (progn ;;;; The Maxima 'flonum can be a CL 'long-float on the Scieneer CL or CLISP, -;;;; but should be the same a 'double-float on other CL implementations. +;;;; but should be the same as 'double-float on other CL implementations. (eval-when (:compile-toplevel :load-toplevel :execute) (setq *read-default-float-format* 'long-float)) diff --git a/src/combin.lisp b/src/combin.lisp index cc12e11..b157f71 100644 --- a/src/combin.lisp +++ b/src/combin.lisp @@ -653,7 +653,7 @@ (rational (setf s (float s))) ((complex rational) - (setf s (coerce s '(complex double-float))))) + (setf s (coerce s '(complex flonum))))) (cond ((bigfloat:< (bigfloat:abs (bigfloat:* s s)) (bigfloat:epsilon s)) ;; abs(s)^2 < epsilon, use the expansion zeta(s) = -1/2-1/2*log(2*%pi)*s (bigfloat:+ -1/2 diff --git a/src/comm.lisp b/src/comm.lisp index 051e13d..ae3a8bf 100644 --- a/src/comm.lisp +++ b/src/comm.lisp @@ -1194,7 +1194,7 @@ ;; bits in n. Then log(n) = log(2^m) + log(n/2^m). ;; n/2^m is approximately 1, so converting that to a ;; float is no problem. log(2^m) = m * log(2). - (+ (* m (log 2d0)) + (+ (* m (log 2e0)) (log (float (/ n (ash 1 m))))))))) (($ratnump n) ;; float(log(n/m)) where n and m are integers. Try computing @@ -1215,7 +1215,7 @@ (let* ((size (max (integer-length re) (integer-length im))) (scale (ash 1 size))) - (+ (* size (log 2d0)) + (+ (* size (log 2e0)) (log (complex (float (/ re scale)) (float (/ im scale)))))))))) (t diff --git a/src/cpoly.lisp b/src/cpoly.lisp index 0ecfd5b..b267ae6 100644 --- a/src/cpoly.lisp +++ b/src/cpoly.lisp @@ -1029,7 +1029,7 @@ (do ((dx x) (df) (f)) - ((fpgreaterp (intofp 5d-3) + ((fpgreaterp (intofp 5e-3) (fpabs (fpquotient dx x))) x) (setq f (aref *shr-sl* 0) diff --git a/src/csimp2.lisp b/src/csimp2.lisp index e9c89f3..0c512f7 100644 --- a/src/csimp2.lisp +++ b/src/csimp2.lisp @@ -561,7 +561,7 @@ (let ((result (let ((c (* (sqrt (* 2 (float pi))) (exp (slatec::d9lgmc a))))) - (let ((v (expt a (- (* .5d0 a) 0.25d0)))) + (let ((v (expt a (- (* .5e0 a) 0.25e0)))) (* v (/ v (exp a)) c))))) diff --git a/src/ellipt.lisp b/src/ellipt.lisp index 8a35100..e29b953 100644 --- a/src/ellipt.lisp +++ b/src/ellipt.lisp @@ -208,7 +208,7 @@ c (/ (- a b) 2)))) ;; WARNING: This seems to have accuracy problems when u is complex. I -;; (rtoy) do not know why. For example (jacobi-agm #c(1d0 1d0) .7d0) +;; (rtoy) do not know why. For example (jacobi-agm #c(1e0 1e0) .7e0) ;; returns ;; ;; #C(1.134045970915582 0.3522523454566013) @@ -237,7 +237,7 @@ (first agm-data) (declare (ignore b c)) (* a u (ash 1 n)))) - (phi1 0d0)) + (phi1 0e0)) (dolist (agm agm-data) (destructuring-bind (n a b c) agm @@ -2731,7 +2731,7 @@ first kind: (cond ((= m 0) (if (maxima::$bfloatp m) (maxima::$bfloat (maxima::div 'maxima::$%pi 2)) - (float (/ pi 2) 1d0))) + (float (/ pi 2) 1e0))) (t (bf-rf 0 (- 1 m) 1)))) @@ -2760,11 +2760,11 @@ first kind: (cond ((= m 0) (if (typep m 'bigfloat) (bigfloat (maxima::$bfloat (maxima::div 'maxima::$%pi 2))) - (float (/ pi 2) 1d0))) + (float (/ pi 2) 1e0))) ((= m 1) (if (typep m 'bigfloat) (bigfloat 1) - 1d0)) + 1e0)) (t (let ((m1 (- 1 m))) (- (bf-rf 0 m1 1) diff --git a/src/gamma.lisp b/src/gamma.lisp index 879d13b..822ffc3 100644 --- a/src/gamma.lisp +++ b/src/gamma.lisp @@ -295,7 +295,7 @@ (expt (/ 2 pival) (/ (- 1 (cos (* pival z))) 4)) - (expt 2d0 (/ z 2)) + (expt 2e0 (/ z 2)) (gamma-lanczos (+ 1 (/ z 2)))))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; @@ -728,7 +728,7 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defvar *gamma-incomplete-maxit* 10000) -(defvar *gamma-incomplete-eps* (* 2 double-float-epsilon)) +(defvar *gamma-incomplete-eps* (* 2 flonum-epsilon)) (defvar *gamma-incomplete-min* 1.0e-32) (defvar *gamma-radius* 1.0 @@ -2806,7 +2806,7 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defvar *fresnel-maxit* 1000) -(defvar *fresnel-eps* (* 2 double-float-epsilon)) +(defvar *fresnel-eps* (* 2 flonum-epsilon)) (defvar *fresnel-min* 1e-32) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; diff --git a/src/nparse.lisp b/src/nparse.lisp index dcda64f..cbdbb21 100644 --- a/src/nparse.lisp +++ b/src/nparse.lisp @@ -332,7 +332,7 @@ ;; For example, 1.234b1000 is converted by computing ;; bfloat(1234)*10b0^(1000-3). Extra precision is used ;; during the bfloat computations. - (let* ((extra-prec (+ *fast-bfloat-extra-bits* (ceiling (log exp 2d0)))) + (let* ((extra-prec (+ *fast-bfloat-extra-bits* (ceiling (log exp 2e0)))) (fpprec (+ fpprec extra-prec)) (mant (+ (* int-part (expt 10 frac-len)) frac-part)) (bf-mant (bcons (intofp mant))) diff --git a/src/numeric.lisp b/src/numeric.lisp index 5a9ec0f..bd136f7 100644 --- a/src/numeric.lisp +++ b/src/numeric.lisp @@ -5,7 +5,7 @@ ;;; support double-float, (complex double-float) and Maxima bfloat and ;;; complex bfloat arithmetic, without having to write separate ;;; versions for each. Of course, specially written versions for -;;; double-float and (cmoplex double-float) will probably be much +;;; double-float and (complex double-float) will probably be much ;;; faster, but this allows users to write just one routine that can ;;; handle all of the data types in a more "natural" Lisp style. @@ -1736,7 +1736,7 @@ (defmethod float ((x bigfloat) (y cl:float)) - (if (typep y 'double-float) + (if (typep y 'maxima::flonum) (maxima::fp2flo (real-value x)) (fp2single (real-value x)))) @@ -1941,6 +1941,9 @@ ;; (coerce bigfloat foo) (cond ((subtypep type 'cl:float) (float obj (cl:coerce 0 type))) + ((subtypep type '(cl:complex long-float)) + (cl:complex (float (realpart obj) 1l0) + (float (imagpart obj) 1l0))) ((subtypep type '(cl:complex double-float)) (cl:complex (float (realpart obj) 1d0) (float (imagpart obj) 1d0))) @@ -1951,15 +1954,18 @@ ;; What should we do here? Return a ;; complex-bigfloat? A complex double-float? ;; complex single-float? I arbitrarily select - ;; complex double-float for now. - (cl:complex (float (realpart obj) 1d0) - (float (imagpart obj) 1d0))) + ;; complex maxima:flonum for now. + (cl:complex (float (realpart obj) 1.0) + (float (imagpart obj) 1.0))) (t (coerce-error)))) ((typep obj 'complex-bigfloat) ;; (coerce complex-bigfloat foo) (cond ((subtypep type 'complex-bigfloat) obj) + ((subtypep type '(cl:complex long-float)) + (cl:complex (float (realpart obj) 1l0) + (float (imagpart obj) 1l0))) ((subtypep type '(cl:complex double-float)) (cl:complex (float (realpart obj) 1d0) (float (imagpart obj) 1d0))) diff --git a/src/numerical/f2cl-lib.lisp b/src/numerical/f2cl-lib.lisp index 4705118..0e1ee5d 100644 --- a/src/numerical/f2cl-lib.lisp +++ b/src/numerical/f2cl-lib.lisp @@ -380,7 +380,13 @@ is not included") (double-float (truncate (the (double-float #.(float most-negative-fixnum 1d0) #.(float most-positive-fixnum 1d0)) - x))))) + x))) + #+clisp + (long-float + (truncate (the (long-float #.(float most-negative-fixnum 1l0) + #.(float most-positive-fixnum 1l0)) + x))) + )) #+(or cmu scl) (defun int (x) @@ -400,7 +406,20 @@ is not included") (the integer4 (truncate (the (double-float #.(float (- (ash 1 31)) 1d0) #.(float (1- (ash 1 31)) 1d0)) - x)))))) + x)))) + #+scl + (long-float + (the integer4 + (truncate (the (long-float #.(float (- (ash 1 31)) 1l0) + #.(float (1- (ash 1 31)) 1l0)) + x)))) + #+(and cmu double-double) + (kernel:double-double-float + (the integer4 + (truncate (the (kernel:double-double-float #.(float (- (ash 1 31)) 1w0) + #.(float (1- (ash 1 31)) 1w0)) + x)))) + )) (defun ifix (x) @@ -512,7 +531,13 @@ is not included") (if (> r 0) (- r 1) (+ r 1)) - r))))) + r))) + #+double-double + (kernel:double-double-float + (locally + (declare (optimize (space 0) (speed 3))) + (values (ftruncate (the kernel:double-double-float x))))) + )) #-cmu @@ -526,7 +551,11 @@ is not included") (double-float (locally (declare (optimize (space 0) (speed 3))) - (values (ftruncate (the double-float x))))))) + (values (ftruncate (the double-float x))))) + (long-float + (locally + (declare (optimize (space 0) (speed 3))) + (values (ftruncate (the long-float x))))))) (defun dint (x) (aint x)) @@ -784,6 +813,12 @@ is not included") (sqrt (the (or (single-float (0f0)) (member 0f0)) x))) (double-float (sqrt (the (or (double-float (0d0)) (member 0d0)) x))) + #+(or scl clisp) + (long-float + (sqrt (the (or (long-float (0l0)) (member 0l0)) x))) + #+(and cmu double-double) + (kernel:double-double-float + (sqrt (the (or (kernel:double-double-float (0w0)) (member 0w0)) x))) (t (sqrt x)))) @@ -793,6 +828,12 @@ is not included") (log (the (or (single-float (0f0)) (member 0f0)) x))) (double-float (log (the (or (double-float (0d0)) (member 0d0)) x))) + #+(or scl clisp) + (long-float + (log (the (or (long-float (0l0)) (member 0l0)) x))) + #+(and cmu double-double) + (kernel:double-double-float + (log (the (or (kernel:double-double-float (0w0)) (member 0w0)) x))) (t (log x)))) @@ -830,9 +871,21 @@ is not included") (log (the (or (single-float (0.0f0)) (member 0f0)) x) 10f0)) (double-float (log (the (or (double-float (0.0d0)) (member 0d0)) x) 10d0)) + #+(or scl clisp) + (long-float + (log (the (or (long-float (0.0l0)) (member 0l0)) x) 10l0)) + #+(and cmu double-double) + (kernel:double-double-float + (log (the (or (kernel:double-double-float (0.0w0)) (member 0w0)) x) 10w0)) (t (/ (log x) (typecase x + #+(and cmu double-double) + ((complex kernel:double-double-float) + 10w0) + #+(or scl clisp) + ((complex long-float) + 10l0) ((complex double-float) 10d0) ((complex single-float) diff --git a/src/rat3e.lisp b/src/rat3e.lisp index f23f0d6..792db5d 100644 --- a/src/rat3e.lisp +++ b/src/rat3e.lisp @@ -863,7 +863,7 @@ (defmvar $ratprint t) -(defmvar $ratepsilon 2d-15) +(defmvar $ratepsilon 2e-15) ;; This control of conversion from float to rational appears to be explained ;; nowhere. - RJF diff --git a/src/specfn.lisp b/src/specfn.lisp index 51e16ad..97fb6ab 100644 --- a/src/specfn.lisp +++ b/src/specfn.lisp @@ -539,15 +539,15 @@ ;; From http://homepages.physik.uni-muenchen.de/~Winitzki/papers/ (defun init-lambert-w (z) - (let ((A 2.344d0) (B 0.8842d0) (C 0.9294d0) (D 0.5106d0) (E -1.213d0) + (let ((a 2.344e0) (b 0.8842e0) (c 0.9294e0) (d 0.5106e0) (e -1.213e0) (y (sqrt (+ (* 2 %e-val z ) 2)) ) ) ; y=sqrt(2*%e*z+2) ; w = (2*log(1+B*y)-log(1+C*log(1+D*y))+E)/(1+1/(2*log(1+B*y)+2*A) (/ (+ (* 2 (log (+ 1 (* b y)))) - (* -1 (log (+ 1 (* C (log (+ 1 (* D y))))))) - E) + (* -1 (log (+ 1 (* c (log (+ 1 (* d y))))))) + e) (+ 1 - (/ 1 (+ (* 2 (log (+ 1 (* B y)))) (* 2 A))))))) + (/ 1 (+ (* 2 (log (+ 1 (* b y)))) (* 2 a))))))) ;; Algorithm based in part on ;; @@ -573,7 +573,7 @@ ;; number of digits correct at step k is roughly 3 times the number ;; which were correct at step k-1. -(defun lambert-w (z &key (maxiter 100) (prec 1d-14)) +(defun lambert-w (z &key (maxiter 100) (prec 1e-14)) (let ((w (init-lambert-w z))) (dotimes (k maxiter) (let* ((we (* w (exp w))) ----------------------------------------------------------------------- Summary of changes: share/hypergeometric/nfloat.lisp | 2 +- share/lapack/dgemm.lisp | 8 ++-- share/lapack/dgeqrf.lisp | 8 ++-- share/lbfgs/maxima-lbfgs.lisp | 14 ++++---- src/bessel.lisp | 28 +++++++++--------- src/clmacs.lisp | 4 +- src/combin.lisp | 2 +- src/comm.lisp | 4 +- src/cpoly.lisp | 2 +- src/csimp2.lisp | 2 +- src/ellipt.lisp | 10 +++--- src/gamma.lisp | 6 ++-- src/nparse.lisp | 2 +- src/numeric.lisp | 16 +++++++--- src/numerical/f2cl-lib.lisp | 61 +++++++++++++++++++++++++++++++++++-- src/rat3e.lisp | 2 +- src/specfn.lisp | 10 +++--- 17 files changed, 120 insertions(+), 61 deletions(-) hooks/post-receive -- Maxima, A Computer Algebra System |