/[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.55 by rtoy, Thu Jan 31 19:12:40 2008 UTC revision 1.55.4.2 by rtoy, Thu Dec 18 21:50:18 2008 UTC
# Line 59  Line 59 
59  ;;; Please refer to the Unix man pages for details about these routines.  ;;; Please refer to the Unix man pages for details about these routines.
60    
61  ;;; Trigonometric.  ;;; Trigonometric.
62  #-x86 (def-math-rtn "sin" 1)  #-(and x86 (not sse2))
63  #-x86 (def-math-rtn "cos" 1)  (progn
64  #-x86 (def-math-rtn "tan" 1)    ;; For x86 (without sse2), we can use x87 instructions to implement
65      ;; these.  With sse2, we don't currently support that, so these
66      ;; should be disabled.
67      (def-math-rtn "sin" 1)
68      (def-math-rtn "cos" 1)
69      (def-math-rtn "tan" 1)
70      (def-math-rtn "atan" 1)
71      (def-math-rtn "atan2" 2))
72  (def-math-rtn "asin" 1)  (def-math-rtn "asin" 1)
73  (def-math-rtn "acos" 1)  (def-math-rtn "acos" 1)
 #-x86 (def-math-rtn "atan" 1)  
 #-x86 (def-math-rtn "atan2" 2)  
74  (def-math-rtn "sinh" 1)  (def-math-rtn "sinh" 1)
75  (def-math-rtn "cosh" 1)  (def-math-rtn "cosh" 1)
76  (def-math-rtn "tanh" 1)  (def-math-rtn "tanh" 1)
# Line 74  Line 79 
79  (def-math-rtn "atanh" 1)  (def-math-rtn "atanh" 1)
80    
81  ;;; Exponential and Logarithmic.  ;;; Exponential and Logarithmic.
82  #-x86 (def-math-rtn "exp" 1)  #-(and x86 (not sse2))
83  #-x86 (def-math-rtn "log" 1)  (progn
84  #-x86 (def-math-rtn "log10" 1)    (def-math-rtn "exp" 1)
85      (def-math-rtn "log" 1)
86      (def-math-rtn "log10" 1))
87    
88  (def-math-rtn "pow" 2)  (def-math-rtn "pow" 2)
89  #-(or x86 sparc-v7 sparc-v8 sparc-v9) (def-math-rtn "sqrt" 1)  #-(or x86 sparc-v7 sparc-v8 sparc-v9)
90    (def-math-rtn "sqrt" 1)
91  (def-math-rtn "hypot" 2)  (def-math-rtn "hypot" 2)
 #-(or hpux x86) (def-math-rtn "log1p" 1)  
92    
93  #+x86 ;; These are needed for use by byte-compiled files.  ;; Don't want log1p to use the x87 instruction.
94    #-(or hpux (and x86 (not sse2)))
95    (def-math-rtn "log1p" 1)
96    
97    ;; These are needed for use by byte-compiled files.  But don't use
98    ;; these with sse2 since we don't support using the x87 instructions
99    ;; here.
100    #+(and x86 (not sse2))
101  (progn  (progn
102    #+nil    #+nil
103    (defun %sin (x)    (defun %sin (x)
# Line 204  Line 219 
219    
220  )  )
221    
222  #+ppc  #+(or ppc sse2)
223  (progn  (progn
224  (declaim (inline %%sin %%cos %%tan))  (declaim (inline %%sin %%cos %%tan))
225  (macrolet ((frob (alien-name lisp-name)  (macrolet ((frob (alien-name lisp-name)
# Line 266  Line 281 
281                      (if (evenp n)                      (if (evenp n)
282                          (,tan reduced)                          (,tan reduced)
283                          (- (/ (,tan reduced)))))))))))                          (- (/ (,tan reduced)))))))))))
284    #+x86    ;; Don't want %sin-quick and friends with sse2.
285      #+(and x86 (not sse2))
286    (frob %sin-quick %cos-quick %tan-quick)    (frob %sin-quick %cos-quick %tan-quick)
287    #+ppc    #+(or ppc sse2)
288    (frob %%sin %%cos %%tan))    (frob %%sin %%cos %%tan))
289    
290    
# Line 287  Line 303 
303    
304  (defparameter *intexp-maximum-exponent* 10000)  (defparameter *intexp-maximum-exponent* 10000)
305    
306    (define-condition intexp-limit-error (error)
307      ((base :initarg :base :reader intexp-base)
308       (power :initarg :power :reader intexp-power))
309      (:report (lambda (condition stream)
310                 (format stream "The absolute value of ~S exceeds limit ~S."
311                         (intexp-power condition)
312                         *intexp-maximum-exponent*))))
313    
314  ;;; This function precisely calculates base raised to an integral power.  It  ;;; This function precisely calculates base raised to an integral power.  It
315  ;;; 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
316  ;;; 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 300  Line 324 
324      (return-from intexp base))      (return-from intexp base))
325    
326    (when (> (abs power) *intexp-maximum-exponent*)    (when (> (abs power) *intexp-maximum-exponent*)
327      (cerror "Continue with calculation."      ;; Allow user the option to continue with calculation, possibly
328              "The absolute value of ~S exceeds ~S."      ;; increasing the limit to the given power.
329              power '*intexp-maximum-exponent* base power))      (restart-case
330            (error 'intexp-limit-error
331                   :base base
332                   :power power)
333          (continue ()
334            :report "Continue with calculation")
335          (new-limit ()
336            :report "Continue with calculation, update limit"
337            (setq *intexp-maximum-exponent* power))))
338    (cond ((minusp power)    (cond ((minusp power)
339           (/ (intexp base (- power))))           (/ (intexp base (- power))))
340          ((eql base 2)          ((eql base 2)
# Line 455  Line 487 
487                                     (coerce (* pow (%cos y*pi)) rtype)                                     (coerce (* pow (%cos y*pi)) rtype)
488                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
489        (declare (inline real-expt))        (declare (inline real-expt))
490          ;; This is really messy and should be cleaned up.  The easiest
491          ;; way to see if we're doing what we should is the macroexpand
492          ;; the number-dispatch and check each branch.
493          ;;
494          ;; We try to apply the rule of float precision contagion (CLHS
495          ;; 12.1.4.4): the result has the same precision has the most
496          ;; precise argument.
497        (number-dispatch ((base number) (power number))        (number-dispatch ((base number) (power number))
498          (((foreach fixnum (or bignum ratio) (complex rational)) integer)          (((foreach fixnum (or bignum ratio) (complex rational))
499              integer)
500           (intexp base power))           (intexp base power))
501          (((foreach single-float double-float) rational)          (((foreach single-float double-float)
502              rational)
503           (real-expt base power '(dispatch-type base)))           (real-expt base power '(dispatch-type base)))
504          (((foreach fixnum (or bignum ratio) single-float)          (((foreach fixnum (or bignum ratio) single-float)
505            (foreach ratio single-float))            (foreach ratio single-float))
# Line 469  Line 510 
510          ((double-float single-float)          ((double-float single-float)
511           (real-expt base power 'double-float))           (real-expt base power 'double-float))
512          #+double-double          #+double-double
513          (((foreach fixnum (or bignum ratio) single-float double-float double-double-float)          (((foreach fixnum (or bignum ratio) single-float double-float
514                       double-double-float)
515            double-double-float)            double-double-float)
516           (dd-%pow (coerce base 'double-double-float) power))           (dd-%pow (coerce base 'double-double-float) power))
517          #+double-double          #+double-double
518          ((double-double-float          ((double-double-float
519            (foreach fixnum (or bignum ratio) single-float double-float))            (foreach fixnum (or bignum ratio) single-float double-float))
520           (dd-%pow base (coerce power 'double-double-float)))           (dd-%pow base (coerce power 'double-double-float)))
521          (((foreach (complex rational) (complex float)) rational)          (((foreach (complex rational) (complex single-float) (complex double-float)
522                       #+double-double (complex double-double-float))
523              rational)
524           (* (expt (abs base) power)           (* (expt (abs base) power)
525              (cis (* power (phase base)))))              (cis (* power (phase base)))))
526          (((foreach fixnum (or bignum ratio) single-float double-float          #+double-double
527                     #+double-double double-double-float)          ((double-double-float
528            complex)            complex)
529           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
530               (* base power)               (* base power)
531                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
532            (((foreach fixnum (or bignum ratio) single-float double-float)
533              (foreach (complex double-float)))
534             ;; Result should have double-float accuracy.  Use log2 in
535             ;; case the base won't fit in a double-float.
536             (if (and (zerop base) (plusp (realpart power)))
537                 (* base power)
538                 (exp (* power (* (log2 base) (log 2d0))))))
539            ((double-float
540              (foreach (complex rational) (complex single-float)))
541             (if (and (zerop base) (plusp (realpart power)))
542                 (* base power)
543                 (exp (* power (log base)))))
544            #+double-double
545            (((foreach fixnum (or bignum ratio) single-float double-float)
546              (foreach (complex double-double-float)))
547             ;; Result should have double-double-float accuracy.  Use log2
548             ;; in case the base won't fit in a double-float.
549             (if (and (zerop base) (plusp (realpart power)))
550                 (* base power)
551                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
552            (((foreach fixnum (or bignum ratio) single-float)
553              (foreach (complex single-float)))
554             (if (and (zerop base) (plusp (realpart power)))
555                 (* base power)
556                 (exp (* power (log base)))))
557            (((foreach (complex rational) (complex single-float))
558              (foreach single-float (complex single-float)))
559             (if (and (zerop base) (plusp (realpart power)))
560                 (* base power)
561                 (exp (* power (log base)))))
562            (((foreach (complex rational) (complex single-float))
563              (foreach double-float (complex double-float)))
564             (if (and (zerop base) (plusp (realpart power)))
565                 (* base power)
566                 (exp (* power (log (coerce base '(complex double-float)))))))
567            #+double-double
568            (((foreach (complex rational) (complex single-float))
569              (foreach double-double-float (complex double-double-float)))
570             (if (and (zerop base) (plusp (realpart power)))
571                 (* base power)
572                 (exp (* power (log (coerce base '(complex double-double-float)))))))
573            (((foreach (complex double-float))
574              (foreach single-float double-float (complex single-float)
575                       (complex double-float)))
576             (if (and (zerop base) (plusp (realpart power)))
577                 (* base power)
578               (exp (* power (log base)))))               (exp (* power (log base)))))
579          (((foreach (complex float) (complex rational))          #+double-double
580            (foreach complex double-float single-float #+double-double double-double-float))          (((foreach (complex double-float))
581              (foreach double-double-float (complex double-double-float)))
582             (if (and (zerop base) (plusp (realpart power)))
583                 (* base power)
584                 (exp (* power (log (coerce base '(complex double-double-float)))))))
585            #+double-double
586            (((foreach (complex double-double-float))
587              (foreach float (complex float)))
588           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
589               (* base power)               (* base power)
590               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))

Legend:
Removed from v.1.55  
changed lines
  Added in v.1.55.4.2

  ViewVC Help
Powered by ViewVC 1.1.5