/[cmucl]/src/code/irrat.lisp
ViewVC logotype

Diff of /src/code/irrat.lisp

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

revision 1.45.2.1 by rtoy, Fri Jun 9 16:04:57 2006 UTC revision 1.64 by rtoy, Mon Feb 28 17:14:45 2011 UTC
# Line 16  Line 16 
16  ;;;  ;;;
17    
18  (in-package "KERNEL")  (in-package "KERNEL")
19    (intl:textdomain "cmucl")
20    
21    
22  ;;;; Random constants, utility functions, and macros.  ;;;; Random constants, utility functions, and macros.
# Line 49  Line 50 
50       (,function ,var))       (,function ,var))
51      #+double-double      #+double-double
52      ((double-double-float)      ((double-double-float)
53       ;; A hack until we write double-double-float versions of these       (,(symbolicate "DD-" function) ,var))))
      ;; special functions.  
      (kernel:make-double-double-float (,function (kernel:double-double-hi ,var))  
                                       0d0))))  
54    
55  ); eval-when (compile load eval)  ); eval-when (compile load eval)
56    
# Line 62  Line 60 
60  ;;; Please refer to the Unix man pages for details about these routines.  ;;; Please refer to the Unix man pages for details about these routines.
61    
62  ;;; Trigonometric.  ;;; Trigonometric.
63  #-x86 (def-math-rtn "sin" 1)  #-(and x86 (not sse2))
64  #-x86 (def-math-rtn "cos" 1)  (progn
65  #-x86 (def-math-rtn "tan" 1)    ;; For x86 (without sse2), we can use x87 instructions to implement
66      ;; these.  With sse2, we don't currently support that, so these
67      ;; should be disabled.
68      (def-math-rtn "sin" 1)
69      (def-math-rtn "cos" 1)
70      (def-math-rtn "tan" 1)
71      (def-math-rtn "atan" 1)
72      (def-math-rtn "atan2" 2))
73  (def-math-rtn "asin" 1)  (def-math-rtn "asin" 1)
74  (def-math-rtn "acos" 1)  (def-math-rtn "acos" 1)
 #-x86 (def-math-rtn "atan" 1)  
 #-x86 (def-math-rtn "atan2" 2)  
75  (def-math-rtn "sinh" 1)  (def-math-rtn "sinh" 1)
76  (def-math-rtn "cosh" 1)  (def-math-rtn "cosh" 1)
77  (def-math-rtn "tanh" 1)  (def-math-rtn "tanh" 1)
# Line 77  Line 80 
80  (def-math-rtn "atanh" 1)  (def-math-rtn "atanh" 1)
81    
82  ;;; Exponential and Logarithmic.  ;;; Exponential and Logarithmic.
83  #-x86 (def-math-rtn "exp" 1)  #-(and x86 (not sse2))
84  #-x86 (def-math-rtn "log" 1)  (progn
85  #-x86 (def-math-rtn "log10" 1)    (def-math-rtn "exp" 1)
86      (def-math-rtn "log" 1)
87      (def-math-rtn "log10" 1))
88    
89  (def-math-rtn "pow" 2)  (def-math-rtn "pow" 2)
90  #-(or x86 sparc-v7 sparc-v8 sparc-v9) (def-math-rtn "sqrt" 1)  #-(or x86 sparc-v7 sparc-v8 sparc-v9)
91    (def-math-rtn "sqrt" 1)
92  (def-math-rtn "hypot" 2)  (def-math-rtn "hypot" 2)
 #-(or hpux x86) (def-math-rtn "log1p" 1)  
93    
94  #+x86 ;; These are needed for use by byte-compiled files.  ;; Don't want log1p to use the x87 instruction.
95    #-(or hpux (and x86 (not sse2)))
96    (def-math-rtn "log1p" 1)
97    
98    ;; These are needed for use by byte-compiled files.  But don't use
99    ;; these with sse2 since we don't support using the x87 instructions
100    ;; here.
101    #+(and x86 (not sse2))
102  (progn  (progn
103      #+nil
104    (defun %sin (x)    (defun %sin (x)
105      (declare (double-float x)      (declare (double-float x)
106               (values double-float))               (values double-float))
# Line 95  Line 109 
109      (declare (double-float x)      (declare (double-float x)
110               (values double-float))               (values double-float))
111      (%sin-quick x))      (%sin-quick x))
112      #+nil
113    (defun %cos (x)    (defun %cos (x)
114      (declare (double-float x)      (declare (double-float x)
115               (values double-float))               (values double-float))
# Line 103  Line 118 
118      (declare (double-float x)      (declare (double-float x)
119               (values double-float))               (values double-float))
120      (%cos-quick x))      (%cos-quick x))
121      #+nil
122    (defun %tan (x)    (defun %tan (x)
123      (declare (double-float x)      (declare (double-float x)
124               (values double-float))               (values double-float))
# Line 163  Line 179 
179    
180  ;; As above for x86.  It also seems to be needed to handle  ;; As above for x86.  It also seems to be needed to handle
181  ;; constant-folding in the compiler.  ;; constant-folding in the compiler.
182  #+sparc  #+(or sparc (and x86 sse2))
183  (progn  (progn
184    (defun %sqrt (x)    (defun %sqrt (x)
185      (declare (double-float x)      (declare (double-float x)
# Line 171  Line 187 
187      (%sqrt x))      (%sqrt x))
188    )    )
189    
190    ;;; The standard libm routines for sin, cos, and tan on x86 (Linux)
191    ;;; and ppc are not very accurate for large arguments when compared to
192    ;;; sparc (and maxima).  This is basically caused by the fact that
193    ;;; those libraries do not do an accurate argument reduction.  The
194    ;;; following functions use some routines Sun's free fdlibm library to
195    ;;; do accurate reduction.  Then we call the standard C functions (or
196    ;;; vops for x86) on the reduced argument.  This produces much more
197    ;;; accurate values.
198    
199    #+(or ppc x86)
200    (progn
201    (declaim (inline %%ieee754-rem-pi/2))
202    ;; Basic argument reduction routine.  It returns two values: n and y
203    ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
204    ;; which octant the arg lies.  Y is actually computed in two parts,
205    ;; y[0] and y[1] such that the sum is y, for accuracy.
206    
207    (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
208      (x double-float)
209      (y (* double-float)))
210    
211    ;; Same as above, but instead of needing to pass an array in, the
212    ;; output array is broken up into two output values instead.  This is
213    ;; easier for the user, and we don't have to wrap calls with
214    ;; without-gcing.
215    (declaim (inline %ieee754-rem-pi/2))
216    (alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
217      (x double-float)
218      (y0 double-float :out)
219      (y1 double-float :out))
220    
221    )
222    
223    #+(or ppc sse2)
224    (progn
225    (declaim (inline %%sin %%cos %%tan))
226    (macrolet ((frob (alien-name lisp-name)
227                 `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
228                    (x double-float))))
229      (frob "sin" %%sin)
230      (frob "cos" %%cos)
231      (frob "tan" %%tan))
232    )
233    
234    #+(or ppc x86)
235    (macrolet
236        ((frob (sin cos tan)
237           `(progn
238              ;; In all of the routines below, we just compute the sum of
239              ;; y0 and y1 and use that as the (reduced) argument for the
240              ;; trig functions.  This is slightly less accurate than what
241              ;; fdlibm does, which calls special functions using y0 and
242              ;; y1 separately, for greater accuracy.  This isn't
243              ;; implemented, and some spot checks indicate that what we
244              ;; have here is accurate.
245              ;;
246              ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
247              ;; probably too restrictive.
248              (defun %sin (x)
249                (declare (double-float x))
250                (if (< (abs x) (/ pi 4))
251                    (,sin x)
252                    ;; Argument reduction needed
253                    (multiple-value-bind (n y0 y1)
254                        (%ieee754-rem-pi/2 x)
255                      (let ((reduced (+ y0 y1)))
256                        (case (logand n 3)
257                          (0 (,sin reduced))
258                          (1 (,cos reduced))
259                          (2 (- (,sin reduced)))
260                          (3 (- (,cos reduced))))))))
261              (defun %cos (x)
262                (declare (double-float x))
263                (if (< (abs x) (/ pi 4))
264                    (,cos x)
265                    ;; Argument reduction needed
266                    (multiple-value-bind (n y0 y1)
267                        (%ieee754-rem-pi/2 x)
268                      (let ((reduced (+ y0 y1)))
269                        (case (logand n 3)
270                          (0 (,cos reduced))
271                          (1 (- (,sin reduced)))
272                          (2 (- (,cos reduced)))
273                          (3 (,sin reduced)))))))
274              (defun %tan (x)
275                (declare (double-float x))
276                (if (< (abs x) (/ pi 4))
277                    (,tan x)
278                    ;; Argument reduction needed
279                    (multiple-value-bind (n y0 y1)
280                        (%ieee754-rem-pi/2 x)
281                      (let ((reduced (+ y0 y1)))
282                        (if (evenp n)
283                            (,tan reduced)
284                            (- (/ (,tan reduced)))))))))))
285      ;; Don't want %sin-quick and friends with sse2.
286      #+(and x86 (not sse2))
287      (frob %sin-quick %cos-quick %tan-quick)
288      #+(or ppc sse2)
289      (frob %%sin %%cos %%tan))
290    
291    
292    
293  ;;;; Power functions.  ;;;; Power functions.
294    
# Line 186  Line 304 
304    
305  (defparameter *intexp-maximum-exponent* 10000)  (defparameter *intexp-maximum-exponent* 10000)
306    
307    (define-condition intexp-limit-error (error)
308      ((base :initarg :base :reader intexp-base)
309       (power :initarg :power :reader intexp-power))
310      (:report (lambda (condition stream)
311                 (format stream (intl:gettext "The absolute value of ~S exceeds limit ~S.")
312                         (intexp-power condition)
313                         *intexp-maximum-exponent*))))
314    
315  ;;; This function precisely calculates base raised to an integral power.  It  ;;; This function precisely calculates base raised to an integral power.  It
316  ;;; separates the cases by the sign of power, for efficiency reasons, as powers  ;;; separates the cases by the sign of power, for efficiency reasons, as powers
317  ;;; can be calculated more efficiently if power is a positive integer.  Values  ;;; can be calculated more efficiently if power is a positive integer.  Values
# Line 199  Line 325 
325      (return-from intexp base))      (return-from intexp base))
326    
327    (when (> (abs power) *intexp-maximum-exponent*)    (when (> (abs power) *intexp-maximum-exponent*)
328      (cerror "Continue with calculation."      ;; Allow user the option to continue with calculation, possibly
329              "The absolute value of ~S exceeds ~S."      ;; increasing the limit to the given power.
330              power '*intexp-maximum-exponent* base power))      (restart-case
331            (error 'intexp-limit-error
332                   :base base
333                   :power power)
334          (continue ()
335            :report (lambda (stream)
336                      (write-string (intl:gettext "Continue with calculation") stream)))
337          (new-limit ()
338            :report (lambda (stream)
339                      (write-string (intl:gettext "Continue with calculation, update limit") stream))
340            (setq *intexp-maximum-exponent* (abs power)))))
341    (cond ((minusp power)    (cond ((minusp power)
342           (/ (intexp base (- power))))           (/ (intexp base (- power))))
343          ((eql base 2)          ((eql base 2)
# Line 353  Line 489 
489                                    (complex                                    (complex
490                                     (coerce (* pow (%cos y*pi)) rtype)                                     (coerce (* pow (%cos y*pi)) rtype)
491                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
492        (declare (inline real-expt))        ;; This is really messy and should be cleaned up.  The easiest
493          ;; way to see if we're doing what we should is the macroexpand
494          ;; the number-dispatch and check each branch.
495          ;;
496          ;; We try to apply the rule of float precision contagion (CLHS
497          ;; 12.1.4.4): the result has the same precision has the most
498          ;; precise argument.
499        (number-dispatch ((base number) (power number))        (number-dispatch ((base number) (power number))
500          (((foreach fixnum (or bignum ratio) (complex rational)) integer)          (((foreach fixnum (or bignum ratio) (complex rational))
501              integer)
502           (intexp base power))           (intexp base power))
503          (((foreach single-float double-float) rational)          (((foreach single-float double-float)
504              rational)
505           (real-expt base power '(dispatch-type base)))           (real-expt base power '(dispatch-type base)))
506          (((foreach fixnum (or bignum ratio) single-float)          (((foreach fixnum (or bignum ratio) single-float)
507            (foreach ratio single-float))            (foreach ratio single-float))
# Line 367  Line 511 
511           (real-expt base power 'double-float))           (real-expt base power 'double-float))
512          ((double-float single-float)          ((double-float single-float)
513           (real-expt base power 'double-float))           (real-expt base power 'double-float))
514          (((foreach (complex rational) (complex float)) rational)          #+double-double
515            (((foreach fixnum (or bignum ratio) single-float double-float
516                       double-double-float)
517              double-double-float)
518             (dd-%pow (coerce base 'double-double-float) power))
519            #+double-double
520            ((double-double-float
521              (foreach fixnum (or bignum ratio) single-float double-float))
522             (dd-%pow base (coerce power 'double-double-float)))
523            (((foreach (complex rational) (complex single-float) (complex double-float)
524                       #+double-double (complex double-double-float))
525              rational)
526           (* (expt (abs base) power)           (* (expt (abs base) power)
527              (cis (* power (phase base)))))              (cis (* power (phase base)))))
528          (((foreach fixnum (or bignum ratio) single-float double-float)          #+double-double
529            ((double-double-float
530            complex)            complex)
531           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
532               (* base power)               (* base power)
533                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
534            (((foreach fixnum (or bignum ratio) single-float double-float)
535              (foreach (complex double-float)))
536             ;; Result should have double-float accuracy.  Use log2 in
537             ;; case the base won't fit in a double-float.
538             (if (and (zerop base) (plusp (realpart power)))
539                 (* base power)
540                 (exp (* power (* (log2 base) (log 2d0))))))
541            ((double-float
542              (foreach (complex rational) (complex single-float)))
543             (if (and (zerop base) (plusp (realpart power)))
544                 (* base power)
545                 (exp (* power (log base)))))
546            #+double-double
547            (((foreach fixnum (or bignum ratio) single-float double-float)
548              (foreach (complex double-double-float)))
549             ;; Result should have double-double-float accuracy.  Use log2
550             ;; in case the base won't fit in a double-float.
551             (if (and (zerop base) (plusp (realpart power)))
552                 (* base power)
553                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
554            (((foreach fixnum (or bignum ratio) single-float)
555              (foreach (complex single-float)))
556             (if (and (zerop base) (plusp (realpart power)))
557                 (* base power)
558                 (exp (* power (log base)))))
559            (((foreach (complex rational) (complex single-float))
560              (foreach single-float (complex single-float)))
561             (if (and (zerop base) (plusp (realpart power)))
562                 (* base power)
563                 (exp (* power (log base)))))
564            (((foreach (complex rational) (complex single-float))
565              (foreach double-float (complex double-float)))
566             (if (and (zerop base) (plusp (realpart power)))
567                 (* base power)
568                 (exp (* power (log (coerce base '(complex double-float)))))))
569            #+double-double
570            (((foreach (complex rational) (complex single-float))
571              (foreach double-double-float (complex double-double-float)))
572             (if (and (zerop base) (plusp (realpart power)))
573                 (* base power)
574                 (exp (* power (log (coerce base '(complex double-double-float)))))))
575            (((foreach (complex double-float))
576              (foreach single-float double-float (complex single-float)
577                       (complex double-float)))
578             (if (and (zerop base) (plusp (realpart power)))
579                 (* base power)
580               (exp (* power (log base)))))               (exp (* power (log base)))))
581          (((foreach (complex float) (complex rational))          #+double-double
582            (foreach complex double-float single-float))          (((foreach (complex double-float))
583              (foreach double-double-float (complex double-double-float)))
584             (if (and (zerop base) (plusp (realpart power)))
585                 (* base power)
586                 (exp (* power (log (coerce base '(complex double-double-float)))))))
587            #+double-double
588            (((foreach (complex double-double-float))
589              (foreach float (complex float)))
590           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
591               (* base power)               (* base power)
592               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))
593    
594  ;; Compute the base 2 log of an integer  ;; Log base 2 of a real number.  The result is a either a double-float
595  (defun log2 (x)  ;; or double-double-float number (real or complex, as appropriate),
596    ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n + log2(f).  ;; depending on the type of FLOAT-TYPE.
597    ;;  (defun log2 (x &optional (float-type 1d0))
598    ;; So we grab the top few bits of x and scale that appropriately,    (labels ((log-of-2 (f)
599    ;; take the log of it and add it to n.               ;; log(2), with the precision specified by the type of F
600    (let ((n (integer-length x)))               (number-dispatch ((f real))
601      (if (< n vm:double-float-digits)                 ((double-float)
602          (log (coerce x 'double-float) 2d0)                  #.(log 2d0))
603          (let ((exp (min vm:double-float-digits n))                 #+double-double
604                (f (ldb (byte vm:double-float-digits                 ((double-double-float)
605                              (max 0 (- n vm:double-float-digits)))                  #.(log 2w0))))
606                        x)))             (log-2-pi (f)
607            (+ n (log (scale-float (float f 1d0) (- exp))               ;; log(pi), with the precision specified by the type of F
608                      2d0))))))               (number-dispatch ((f real))
609                   ((double-float)
610                    #.(/ pi (log 2d0)))
611                   #+double-double
612                   ((double-double-float)
613                    #.(/ dd-pi (log 2w0)))))
614               (log1p (x)
615                 ;; log(1+x), with the precision specified by the type of
616                 ;; X
617                 (number-dispatch ((x real))
618                   (((foreach single-float double-float))
619                    (%log1p (float x 1d0)))
620                   #+double-double
621                   ((double-double-float)
622                    (dd-%log1p x))))
623               (log2-bignum (bignum)
624                 ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
625                 ;; + log2(f).
626                 ;;
627                 ;; So we grab the top few bits of x and scale that
628                 ;; appropriately, take the log of it and add it to n.
629                 ;;
630                 ;; Return n and log2(f) separately.
631                 (if (minusp bignum)
632                     (multiple-value-bind (n frac)
633                         (log2-bignum (abs bignum))
634                       (values n (complex frac (log-2-pi float-type))))
635                     (let ((n (integer-length bignum))
636                           (float-bits (float-digits float-type)))
637                       (if (< n float-bits)
638                           (values 0 (log (float bignum float-type)
639                                          (float 2 float-type)))
640                           (let ((exp (min float-bits n))
641                                 (f (ldb (byte float-bits
642                                               (max 0 (- n float-bits)))
643                                         bignum)))
644                             (values n (log (scale-float (float f float-type) (- exp))
645                                            (float 2 float-type)))))))))
646        (etypecase x
647          (float
648           (/ (log (float x float-type)) (log-of-2 float-type)))
649          (ratio
650           (let ((top (numerator x))
651                 (bot (denominator x)))
652             ;; If the number of bits in the numerator and
653             ;; denominator are different, just use the fact
654             ;; log(x/y) = log(x) - log(y).  But to preserve
655             ;; accuracy, we actually do
656             ;; (log2(x)-log2(y))/log2(e)).
657             ;;
658             ;; However, if the numerator and denominator have the
659             ;; same number of bits, implying the quotient is near
660             ;; one, we use log1p(x) = log(1+x). Since the number is
661             ;; rational, we don't lose precision subtracting 1 from
662             ;; it, and converting it to double-float is accurate.
663             (if (= (integer-length top)
664                    (integer-length bot))
665                 (/ (log1p (float (- x 1) float-type))
666                    (log-of-2 float-type))
667                 (multiple-value-bind (top-n top-frac)
668                     (log2-bignum top)
669                   (multiple-value-bind (bot-n bot-frac)
670                       (log2-bignum bot)
671                     (+ (- top-n bot-n)
672                        (- top-frac bot-frac)))))))
673          (integer
674           (multiple-value-bind (n frac)
675               (log2-bignum x)
676             (+ n frac))))))
677    
678  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
679    "Return the logarithm of NUMBER in the base BASE, which defaults to e."    "Return the logarithm of NUMBER in the base BASE, which defaults to e."
# Line 403  Line 681 
681        (cond ((zerop base)        (cond ((zerop base)
682               ;; ANSI spec               ;; ANSI spec
683               base)               base)
             ((and (integerp number) (integerp base)  
                   (plusp number) (plusp base))  
              ;; Let's try to do something nice when both the number  
              ;; and the base are positive integers.  Use the rule that  
              ;; log_b(x) = log_2(x)/log_2(b)  
              (coerce (/ (log2 number) (log2 base)) 'single-float))  
684              ((and (realp number) (realp base))              ((and (realp number) (realp base))
685               ;; CLHS 12.1.4.1 says               ;; CLHS 12.1.4.1 says
686               ;;               ;;
# Line 423  Line 695 
695               ;; This makes (log 17 10.0) = (log 17.0 10) and so on.               ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
696               (number-dispatch ((number real) (base real))               (number-dispatch ((number real) (base real))
697                 ((double-float                 ((double-float
698                   (foreach double-float single-float fixnum bignum ratio))                   (foreach double-float single-float))
699                  (/ (log number) (log (coerce base 'double-float))))                  (/ (log2 number) (log2 base)))
700                 (((foreach single-float fixnum bignum ratio)                 (((foreach fixnum bignum ratio)
701                     (foreach fixnum bignum ratio single-float))
702                    (let* ((result (/ (log2 number) (log2 base))))
703                      ;; Figure out the right result type
704                      (if (realp result)
705                          (coerce result 'single-float)
706                          (coerce result '(complex single-float)))))
707                   (((foreach fixnum bignum ratio)
708                   double-float)                   double-float)
709                    (/ (log2 number) (log2 base)))
710                   ((single-float
711                     (foreach fixnum bignum ratio))
712                    (let* ((result (/ (log2 number) (log2 base))))
713                      ;; Figure out the right result type
714                      (if (realp result)
715                          (coerce result 'single-float)
716                          (coerce result '(complex single-float)))))
717                   ((double-float
718                     (foreach fixnum bignum ratio))
719                    (/ (log2 number) (log2 base)))
720                   ((single-float double-float)
721                  (/ (log (coerce number 'double-float)) (log base)))                  (/ (log (coerce number 'double-float)) (log base)))
722                 (((foreach single-float fixnum bignum ratio)                 #+double-double
723                   (foreach single-float fixnum bignum ratio))                 ((double-double-float
724                     (foreach fixnum bignum ratio))
725                    (/ (log2 number 1w0) (log2 base 1w0)))
726                   #+double-double
727                   ((double-double-float
728                     (foreach double-double-float double-float single-float))
729                    (/ (log number) (log (coerce base 'double-double-float))))
730                   #+double-double
731                   (((foreach fixnum bignum ratio)
732                     double-double-float)
733                    (/ (log2 number 1w0) (log2 base 1w0)))
734                   #+double-double
735                   (((foreach double-float single-float)
736                     double-double-float)
737                    (/ (log (coerce number 'double-double-float)) (log base)))
738                   (((foreach single-float)
739                     (foreach single-float))
740                  ;; Converting everything to double-float helps the                  ;; Converting everything to double-float helps the
741                  ;; cases like (log 17 10) = (/ (log 17) (log 10)).                  ;; cases like (log 17 10) = (/ (log 17) (log 10)).
742                  ;; This is usually handled above, but if we compute (/                  ;; This is usually handled above, but if we compute (/
# Line 489  Line 796 
796                       '(dispatch-type number))))                       '(dispatch-type number))))
797          #+double-double          #+double-double
798          ((double-double-float)          ((double-double-float)
          ;; Hack!  
799           (let ((hi (kernel:double-double-hi number)))           (let ((hi (kernel:double-double-hi number)))
800             (if (< (float-sign hi)             (if (< (float-sign hi) 0d0)
801                    (coerce 0 '(dispatch-type number)))                 (complex (dd-%log (- number)) dd-pi)
802                 (complex (coerce (log (- hi)) 'kernel:double-double-float)                 (dd-%log number))))
                         (coerce pi '(dispatch-type number)))  
                (coerce (%log hi) '(dispatch-type number)))))  
803          ((complex)          ((complex)
804           (complex-log number)))))           (complex-log number)))))
805    
# Line 513  Line 817 
817                   '(dispatch-type number))))                   '(dispatch-type number))))
818      #+double-double      #+double-double
819      ((double-double-float)      ((double-double-float)
820       (multiple-value-bind (hi lo)       (if (minusp number)
821           (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))           (dd-complex-sqrt number)
822         (kernel:make-double-double-float hi lo)))           (multiple-value-bind (hi lo)
823                 (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
824               (kernel:%make-double-double-float hi lo))))
825      ((complex)      ((complex)
826       (complex-sqrt number))))       (complex-sqrt number))))
827    
# Line 525  Line 831 
831  (defun abs (number)  (defun abs (number)
832    "Returns the absolute value of the number."    "Returns the absolute value of the number."
833    (number-dispatch ((number number))    (number-dispatch ((number number))
834      (((foreach single-float double-float fixnum rational))      (((foreach single-float double-float fixnum rational
835                   #+double-double double-double-float))
836       (abs number))       (abs number))
     #+(and nil double-double)  
     ((double-double-float)  
      ;; This is a hack until abs deftransform is working  
      (multiple-value-bind (hi lo)  
          (c::abs-dd (kernel:double-double-hi number) (kernel:double-double-lo number))  
        (kernel:make-double-double-float hi lo)))  
     #+double-double  
     ((double-double-float)  
      ;; This is a hack until abs deftransform is working  
      (let ((hi (kernel:double-double-hi number))  
            (lo (kernel:double-double-lo number)))  
        (declare (double-float hi lo))  
        (when (minusp hi)  
          (setf hi (- hi))  
          (setf lo (- lo)))  
        (kernel:make-double-double-float hi lo)))  
837      ((complex)      ((complex)
838       (let ((rx (realpart number))       (let ((rx (realpart number))
839             (ix (imagpart number)))             (ix (imagpart number)))
# Line 557  Line 848 
848            (%hypot rx ix))            (%hypot rx ix))
849           #+double-double           #+double-double
850           (double-double-float           (double-double-float
851            (error "abs complex double-double-float not implemented!")))))))            (multiple-value-bind (abs^2 scale)
852                  (dd-cssqs number)
853                (scale-float (sqrt abs^2) scale))))))))
854    
855  (defun phase (number)  (defun phase (number)
856    "Returns the angle part of the polar representation of a complex number.    "Returns the angle part of the polar representation of a complex number.
# Line 580  Line 873 
873      #+double-double      #+double-double
874      (double-double-float      (double-double-float
875       (if (minusp (float-sign number))       (if (minusp (float-sign number))
876           (coerce pi 'double-double-float)           dd-pi
877           (kernel:make-double-double-float 0d0 0d0)))           0w0))
878      (complex      (complex
879       (atan (imagpart number) (realpart number)))))       (atan (imagpart number) (realpart number)))))
880    
# Line 616  Line 909 
909  (defun cis (theta)  (defun cis (theta)
910    "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."    "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
911    (if (complexp theta)    (if (complexp theta)
912        (error "Argument to CIS is complex: ~S" theta)        (error (intl:gettext "Argument to CIS is complex: ~S") theta)
913        (complex (cos theta) (sin theta))))        (complex (cos theta) (sin theta))))
914    
915  (defun asin (number)  (defun asin (number)
# Line 632  Line 925 
925                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
926           (coerce (%asin (coerce number 'double-float))           (coerce (%asin (coerce number 'double-float))
927                   '(dispatch-type number))                   '(dispatch-type number))
928                   (complex-asin number)))           (complex-asin number)))
929        #+double-double
930        ((double-double-float)
931         (if (or (float-nan-p number)
932                 (and (<= number 1w0)
933                      (>= number -1w0)))
934             (dd-%asin number)
935             (dd-complex-asin number)))
936      ((complex)      ((complex)
937       (complex-asin number))))       (complex-asin number))))
938    
# Line 649  Line 949 
949                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
950           (coerce (%acos (coerce number 'double-float))           (coerce (%acos (coerce number 'double-float))
951                   '(dispatch-type number))                   '(dispatch-type number))
952                   (complex-acos number)))           (complex-acos number)))
953        #+double-double
954        ((double-double-float)
955         (if (or (float-nan-p number)
956                 (and (<= number 1w0)
957                      (>= number -1w0)))
958             (dd-%acos number)
959             (complex-acos number)))
960      ((complex)      ((complex)
961       (complex-acos number))))       (complex-acos number))))
962    
# Line 678  Line 985 
985            (((foreach single-float fixnum bignum ratio)            (((foreach single-float fixnum bignum ratio)
986              (foreach single-float fixnum bignum ratio))              (foreach single-float fixnum bignum ratio))
987             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
988                     'single-float))))                     'single-float))
989              #+double-double
990              ((double-double-float
991                (foreach double-double-float double-float single-float fixnum bignum ratio))
992               (dd-%atan2 y (coerce x 'double-double-float)))
993              #+double-double
994              (((foreach double-float single-float fixnum bignum ratio)
995                double-double-float)
996               (dd-%atan2 (coerce y 'double-double-float) x))))
997        (number-dispatch ((y number))        (number-dispatch ((y number))
998          (handle-reals %atan y)          (handle-reals %atan y)
999          ((complex)          ((complex)
# Line 731  Line 1046 
1046           (complex-acosh number)           (complex-acosh number)
1047           (coerce (%acosh (coerce number 'double-float))           (coerce (%acosh (coerce number 'double-float))
1048                   '(dispatch-type number))))                   '(dispatch-type number))))
1049        #+double-double
1050        ((double-double-float)
1051         (if (< number 1w0)
1052             (complex-acosh number)
1053             (dd-%acosh number)))
1054      ((complex)      ((complex)
1055       (complex-acosh number))))       (complex-acosh number))))
1056    
# Line 748  Line 1068 
1068           (complex-atanh number)           (complex-atanh number)
1069           (coerce (%atanh (coerce number 'double-float))           (coerce (%atanh (coerce number 'double-float))
1070                   '(dispatch-type number))))                   '(dispatch-type number))))
1071        #+double-double
1072        ((double-double-float)
1073         (if (or (> number 1w0)
1074                 (< number -1w0))
1075             (complex-atanh number)
1076             (dd-%atanh (coerce number 'double-double-float))))
1077      ((complex)      ((complex)
1078       (complex-atanh number))))       (complex-atanh number))))
1079    
# Line 817  Line 1143 
1143    
1144  (declaim (inline square))  (declaim (inline square))
1145  (defun square (x)  (defun square (x)
1146    (declare (double-float x))    (declare (float x))
1147    (* x x))    (* x x))
1148    
1149  ;; If you have these functions in libm, perhaps they should be used  ;; If you have these functions in libm, perhaps they should be used
# Line 828  Line 1154 
1154  (defun scalb (x n)  (defun scalb (x n)
1155    "Compute 2^N * X without compute 2^N first (use properties of the    "Compute 2^N * X without compute 2^N first (use properties of the
1156  underlying floating-point format"  underlying floating-point format"
1157    (declare (type double-float x)    (declare (type float x)
1158             (type double-float-exponent n))             (type double-float-exponent n))
1159    (scale-float x n))    (scale-float x n))
1160    
# Line 836  underlying floating-point format" Line 1162  underlying floating-point format"
1162  (defun logb-finite (x)  (defun logb-finite (x)
1163    "Same as logb but X is not infinity and non-zero and not a NaN, so    "Same as logb but X is not infinity and non-zero and not a NaN, so
1164  that we can always return an integer"  that we can always return an integer"
1165    (declare (type double-float x))    (declare (type float x))
1166    (multiple-value-bind (signif expon sign)    (multiple-value-bind (signif expon sign)
1167        (decode-float x)        (decode-float x)
1168      (declare (ignore signif sign))      (declare (ignore signif sign))
# Line 853  For the special cases, the following val Line 1179  For the special cases, the following val
1179     +/- infinity   +infinity     +/- infinity   +infinity
1180     0              -infinity     0              -infinity
1181  "  "
1182    (declare (type double-float x))    (declare (type float x))
1183    (cond ((float-nan-p x)    (cond ((float-nan-p x)
1184           x)           x)
1185          ((float-infinity-p x)          ((float-infinity-p x)
# Line 861  For the special cases, the following val Line 1187  For the special cases, the following val
1187          ((zerop x)          ((zerop x)
1188           ;; The answer is negative infinity, but we are supposed to           ;; The answer is negative infinity, but we are supposed to
1189           ;; signal divide-by-zero, so do the actual division           ;; signal divide-by-zero, so do the actual division
1190           (/ -1.0d0 x)           (/ -1 x)
1191           )           )
1192          (t          (t
1193           (logb-finite x))))           (logb-finite x))))
# Line 928  and Y are coerced to single-float." Line 1254  and Y are coerced to single-float."
1254    
1255  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1256    (declare (number z))    (declare (number z))
1257      #+double-double
1258      (when (typep z '(or double-double-float (complex double-double-float)))
1259        (return-from complex-sqrt (dd-complex-sqrt z)))
1260    (multiple-value-bind (rho k)    (multiple-value-bind (rho k)
1261        (cssqs z)        (cssqs z)
1262      (declare (type (or (member 0d0) (double-float 0d0)) rho)      (declare (type (or (member 0d0) (double-float 0d0)) rho)
# Line 942  Z may be any number, but the result is a Line 1271  Z may be any number, but the result is a
1271            ;; space 0 to get maybe-inline functions inlined.            ;; space 0 to get maybe-inline functions inlined.
1272            (declare (optimize (speed 3) (space 0)))            (declare (optimize (speed 3) (space 0)))
1273    
1274          (if (not (float-nan-p x))          (if (not (locally (declare (optimize (inhibit-warnings 3)))
1275                       (float-nan-p x)))
1276              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1277    
1278          (cond ((oddp k)          (cond ((oddp k)
# Line 1004  This is for use with J /= 0 only when |z Line 1334  This is for use with J /= 0 only when |z
1334    
1335  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1336    (declare (number z))    (declare (number z))
1337      #+double-double
1338      (when (typep z '(or double-double-float (complex double-double-float)))
1339        (return-from complex-log (dd-complex-log-scaled z 0)))
1340    (complex-log-scaled z 0))    (complex-log-scaled z 0))
1341    
1342  ;; Let us note the following "strange" behavior.  atanh 1.0d0 is  ;; Let us note the following "strange" behavior.  atanh 1.0d0 is
# Line 1014  Z may be any number, but the result is a Line 1347  Z may be any number, but the result is a
1347  (defun complex-atanh (z)  (defun complex-atanh (z)
1348    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
1349    (declare (number z))    (declare (number z))
1350      #+double-double
1351      (when (typep z '(or double-double-float (complex double-double-float)))
1352        (return-from complex-atanh (dd-complex-atanh z)))
1353    
1354    (if (and (realp z) (< z -1))    (if (and (realp z) (< z -1))
1355        ;; atanh is continuous in quadrant III in this case.        ;; atanh is continuous in quadrant III in this case.
1356        (complex-atanh (complex z -0f0))        (complex-atanh (complex z -0f0))
# Line 1074  Z may be any number, but the result is a Line 1411  Z may be any number, but the result is a
1411  (defun complex-tanh (z)  (defun complex-tanh (z)
1412    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"
1413    (declare (number z))    (declare (number z))
1414      #+double-double
1415      (when (typep z '(or double-double-float (complex double-double-float)))
1416        (return-from complex-tanh (dd-complex-tanh z)))
1417    
1418    (let ((x (float (realpart z) 1.0d0))    (let ((x (float (realpart z) 1.0d0))
1419          (y (float (imagpart z) 1.0d0)))          (y (float (imagpart z) 1.0d0)))
1420      (locally      (locally
# Line 1157  Z may be any number, but the result is a Line 1498  Z may be any number, but the result is a
1498    
1499  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1500    (declare (number z))    (declare (number z))
1501      #+double-double
1502      (when (typep z '(or double-double-float (complex double-double-float)))
1503        (return-from complex-acos (dd-complex-acos z)))
1504    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1505        ;; acos is continuous in quadrant IV in this case.        ;; acos is continuous in quadrant IV in this case.
1506        (complex-acos (complex z -0f0))        (complex-acos (complex z -0f0))
# Line 1187  Z may be any number, but the result is a Line 1531  Z may be any number, but the result is a
1531    
1532  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1533    (declare (number z))    (declare (number z))
1534      #+double-double
1535      (when (typep z '(or double-double-float (complex double-double-float)))
1536        (return-from complex-asin (dd-complex-asin z)))
1537    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1538        ;; asin is continuous in quadrant IV in this case.        ;; asin is continuous in quadrant IV in this case.
1539        (complex-asin (complex z -0f0))        (complex-asin (complex z -0f0))
# Line 1204  Z may be any number, but the result is a Line 1551  Z may be any number, but the result is a
1551  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1552    (declare (number z))    (declare (number z))
1553    ;; asinh z = -i * asin (i*z)    ;; asinh z = -i * asin (i*z)
1554      #+double-double
1555      (when (typep z '(or double-double-float (complex double-double-float)))
1556        (return-from complex-asinh (dd-complex-asinh z)))
1557    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1558           (result (complex-asin iz)))           (result (complex-asin iz)))
1559      (complex (imagpart result)      (complex (imagpart result)
# Line 1215  Z may be any number, but the result is a Line 1565  Z may be any number, but the result is a
1565  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1566    (declare (number z))    (declare (number z))
1567    ;; atan z = -i * atanh (i*z)    ;; atan z = -i * atanh (i*z)
1568      #+double-double
1569      (when (typep z '(or double-double-float (complex double-double-float)))
1570        (return-from complex-atan (dd-complex-atan z)))
1571    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1572           (result (complex-atanh iz)))           (result (complex-atanh iz)))
1573      (complex (imagpart result)      (complex (imagpart result)
# Line 1226  Z may be any number, but the result is a Line 1579  Z may be any number, but the result is a
1579  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1580    (declare (number z))    (declare (number z))
1581    ;; tan z = -i * tanh(i*z)    ;; tan z = -i * tanh(i*z)
1582      #+double-double
1583      (when (typep z '(or double-double-float (complex double-double-float)))
1584        (return-from complex-tan (dd-complex-tan z)))
1585    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1586           (result (complex-tanh iz)))           (result (complex-tanh iz)))
1587      (complex (imagpart result)      (complex (imagpart result)

Legend:
Removed from v.1.45.2.1  
changed lines
  Added in v.1.64

  ViewVC Help
Powered by ViewVC 1.1.5