# Diff of /src/code/irrat.lisp

revision 1.55 by rtoy, Thu Jan 31 19:12:40 2008 UTC revision 1.55.6.1 by rtoy, Sun Nov 2 13:30:01 2008 UTC
# Line 287  Line 287
287
288  (defparameter *intexp-maximum-exponent* 10000)  (defparameter *intexp-maximum-exponent* 10000)
289
290    (define-condition intexp-limit-error (error)
291      ((base :initarg :base :reader intexp-base)
292       (power :initarg :power :reader intexp-power))
293      (:report (lambda (condition stream)
294                 (format stream "The absolute value of ~S exceeds limit ~S."
295                         (intexp-power condition)
296                         *intexp-maximum-exponent*))))
297
298  ;;; This function precisely calculates base raised to an integral power.  It  ;;; This function precisely calculates base raised to an integral power.  It
299  ;;; 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
300  ;;; 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 308
308      (return-from intexp base))      (return-from intexp base))
309
310    (when (> (abs power) *intexp-maximum-exponent*)    (when (> (abs power) *intexp-maximum-exponent*)
311      (cerror "Continue with calculation."      ;; Allow user the option to continue with calculation, possibly
312              "The absolute value of ~S exceeds ~S."      ;; increasing the limit to the given power.
313              power '*intexp-maximum-exponent* base power))      (restart-case
314            (error 'intexp-limit-error
315                   :base base
316                   :power power)
317          (continue ()
318            :report "Continue with calculation")
319          (new-limit ()
320            :report "Continue with calculation, update limit"
321            (setq *intexp-maximum-exponent* power))))
322    (cond ((minusp power)    (cond ((minusp power)
323           (/ (intexp base (- power))))           (/ (intexp base (- power))))
324          ((eql base 2)          ((eql base 2)
# Line 455  Line 471
471                                     (coerce (* pow (%cos y*pi)) rtype)                                     (coerce (* pow (%cos y*pi)) rtype)
472                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))                                     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
473        (declare (inline real-expt))        (declare (inline real-expt))
474          ;; This is really messy and should be cleaned up.  The easiest
475          ;; way to see if we're doing what we should is the macroexpand
476          ;; the number-dispatch and check each branch.
477          ;;
478          ;; We try to apply the rule of float precision contagion (CLHS
479          ;; 12.1.4.4): the result has the same precision has the most
480          ;; precise argument.
481        (number-dispatch ((base number) (power number))        (number-dispatch ((base number) (power number))
482          (((foreach fixnum (or bignum ratio) (complex rational)) integer)          (((foreach fixnum (or bignum ratio) (complex rational))
483              integer)
484           (intexp base power))           (intexp base power))
485          (((foreach single-float double-float) rational)          (((foreach single-float double-float)
486              rational)
487           (real-expt base power '(dispatch-type base)))           (real-expt base power '(dispatch-type base)))
488          (((foreach fixnum (or bignum ratio) single-float)          (((foreach fixnum (or bignum ratio) single-float)
489            (foreach ratio single-float))            (foreach ratio single-float))
# Line 469  Line 494
494          ((double-float single-float)          ((double-float single-float)
495           (real-expt base power 'double-float))           (real-expt base power 'double-float))
496          #+double-double          #+double-double
497          (((foreach fixnum (or bignum ratio) single-float double-float double-double-float)          (((foreach fixnum (or bignum ratio) single-float double-float
498                       double-double-float)
499            double-double-float)            double-double-float)
500           (dd-%pow (coerce base 'double-double-float) power))           (dd-%pow (coerce base 'double-double-float) power))
501          #+double-double          #+double-double
502          ((double-double-float          ((double-double-float
503            (foreach fixnum (or bignum ratio) single-float double-float))            (foreach fixnum (or bignum ratio) single-float double-float))
504           (dd-%pow base (coerce power 'double-double-float)))           (dd-%pow base (coerce power 'double-double-float)))
505          (((foreach (complex rational) (complex float)) rational)          (((foreach (complex rational) (complex single-float) (complex double-float)
506                       #+double-double (complex double-double-float))
507              rational)
508           (* (expt (abs base) power)           (* (expt (abs base) power)
509              (cis (* power (phase base)))))              (cis (* power (phase base)))))
510          (((foreach fixnum (or bignum ratio) single-float double-float          #+double-double
511                     #+double-double double-double-float)          ((double-double-float
512            complex)            complex)
513           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
514               (* base power)               (* base power)
515                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
516            (((foreach fixnum (or bignum ratio) single-float double-float)
517              (foreach (complex double-float)))
518             ;; Result should have double-float accuracy.  Use log2 in
519             ;; case the base won't fit in a double-float.
520             (if (and (zerop base) (plusp (realpart power)))
521                 (* base power)
522                 (exp (* power (* (log2 base) (log 2d0))))))
523            ((double-float
524              (foreach (complex rational) (complex single-float)))
525             (if (and (zerop base) (plusp (realpart power)))
526                 (* base power)
527               (exp (* power (log base)))))               (exp (* power (log base)))))
528          (((foreach (complex float) (complex rational))          #+double-double
529            (foreach complex double-float single-float #+double-double double-double-float))          (((foreach fixnum (or bignum ratio) single-float double-float)
530              (foreach (complex double-double-float)))
531             ;; Result should have double-double-float accuracy.  Use log2
532             ;; in case the base won't fit in a double-float.
533             (if (and (zerop base) (plusp (realpart power)))
534                 (* base power)
535                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
536            (((foreach fixnum (or bignum ratio) single-float)
537              (foreach (complex single-float)))
538             (if (and (zerop base) (plusp (realpart power)))
539                 (* base power)
540                 (exp (* power (log base)))))
541            (((foreach (complex rational) (complex single-float))
542              (foreach single-float (complex single-float)))
543             (if (and (zerop base) (plusp (realpart power)))
544                 (* base power)
545                 (exp (* power (log base)))))
546            (((foreach (complex rational) (complex single-float))
547              (foreach double-float (complex double-float)))
548             (if (and (zerop base) (plusp (realpart power)))
549                 (* base power)
550                 (exp (* power (log (coerce base '(complex double-float)))))))
551            #+double-double
552            (((foreach (complex rational) (complex single-float))
553              (foreach double-double-float (complex double-double-float)))
554             (if (and (zerop base) (plusp (realpart power)))
555                 (* base power)
556                 (exp (* power (log (coerce base '(complex double-double-float)))))))
557            (((foreach (complex double-float))
558              (foreach single-float double-float (complex single-float)
559                       (complex double-float)))
560             (if (and (zerop base) (plusp (realpart power)))
561                 (* base power)
562                 (exp (* power (log base)))))
563            #+double-double
564            (((foreach (complex double-float))
565              (foreach double-double-float (complex double-double-float)))
566             (if (and (zerop base) (plusp (realpart power)))
567                 (* base power)
568                 (exp (* power (log (coerce base '(complex double-double-float)))))))
569            #+double-double
570            (((foreach (complex double-double-float))
571              (foreach float (complex float)))
572           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
573               (* base power)               (* base power)
574               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))

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