/[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.2.1.2.1 by rtoy, Sat Jun 17 02:59:42 2006 UTC revision 1.45.2.1.2.1.2.2 by rtoy, Thu Jun 29 01:28:02 2006 UTC
# Line 49  Line 49 
49       (,function ,var))       (,function ,var))
50      #+double-double      #+double-double
51      ((double-double-float)      ((double-double-float)
52       ;; 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))))  
53    
54  ); eval-when (compile load eval)  ); eval-when (compile load eval)
55    
# Line 367  Line 364 
364           (real-expt base power 'double-float))           (real-expt base power 'double-float))
365          ((double-float single-float)          ((double-float single-float)
366           (real-expt base power 'double-float))           (real-expt base power 'double-float))
367            #+double-double
368            (((foreach fixnum (or bignum ratio) single-float double-float double-double-float)
369              double-double-float)
370             (dd-%pow (coerce base 'double-double-float) power))
371            #+double-double
372            ((double-double-float
373              (foreach fixnum (or bignum ratio) single-float double-float))
374             (dd-%pow base (coerce power 'double-double-float)))
375          (((foreach (complex rational) (complex float)) rational)          (((foreach (complex rational) (complex float)) rational)
376           (* (expt (abs base) power)           (* (expt (abs base) power)
377              (cis (* power (phase base)))))              (cis (* power (phase base)))))
# Line 491  Line 496 
496          ((double-double-float)          ((double-double-float)
497           ;; Hack!           ;; Hack!
498           (let ((hi (kernel:double-double-hi number)))           (let ((hi (kernel:double-double-hi number)))
499             (if (< (float-sign hi)             (if (< (float-sign hi) 0d0)
500                    (coerce 0 '(dispatch-type number)))                 (complex (dd-%log (- number)) dd-pi)
501                 (complex (coerce (log (- hi)) 'kernel:double-double-float)                 (dd-%log number))))
                         (coerce pi '(dispatch-type number)))  
                (coerce (%log hi) '(dispatch-type number)))))  
502          ((complex)          ((complex)
503           (complex-log number)))))           (complex-log number)))))
504    
# Line 513  Line 516 
516                   '(dispatch-type number))))                   '(dispatch-type number))))
517      #+double-double      #+double-double
518      ((double-double-float)      ((double-double-float)
519       (multiple-value-bind (hi lo)       (if (minusp number)
520           (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))           (dd-complex-sqrt number)
521         (kernel:make-double-double-float hi lo)))           (multiple-value-bind (hi lo)
522                 (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
523               (kernel:make-double-double-float hi lo))))
524      ((complex)      ((complex)
525       (complex-sqrt number))))       (complex-sqrt number))))
526    
# Line 581  Line 586 
586      #+double-double      #+double-double
587      (double-double-float      (double-double-float
588       (if (minusp (float-sign number))       (if (minusp (float-sign number))
589           (coerce pi 'double-double-float)           dd-pi
590           0w0))           0w0))
591      (complex      (complex
592       (atan (imagpart number) (realpart number)))))       (atan (imagpart number) (realpart number)))))
# Line 633  Line 638 
638                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
639           (coerce (%asin (coerce number 'double-float))           (coerce (%asin (coerce number 'double-float))
640                   '(dispatch-type number))                   '(dispatch-type number))
641                   (complex-asin number)))           (complex-asin number)))
642        #+double-double
643        ((double-double-float)
644         (if (or (float-nan-p number)
645                 (and (<= number 1w0)
646                      (>= number -1w0)))
647             (dd-%asin number)
648             (dd-complex-asin number)))
649      ((complex)      ((complex)
650       (complex-asin number))))       (complex-asin number))))
651    
# Line 650  Line 662 
662                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
663           (coerce (%acos (coerce number 'double-float))           (coerce (%acos (coerce number 'double-float))
664                   '(dispatch-type number))                   '(dispatch-type number))
665                   (complex-acos number)))           (complex-acos number)))
666        #+double-double
667        ((double-double-float)
668         (if (or (float-nan-p number)
669                 (and (<= number 1w0)
670                      (>= number -1w0)))
671             (dd-%acos number)
672             (complex-acos number)))
673      ((complex)      ((complex)
674       (complex-acos number))))       (complex-acos number))))
675    
# Line 679  Line 698 
698            (((foreach single-float fixnum bignum ratio)            (((foreach single-float fixnum bignum ratio)
699              (foreach single-float fixnum bignum ratio))              (foreach single-float fixnum bignum ratio))
700             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
701                     'single-float))))                     'single-float))
702              #+double-double
703              ((double-double-float
704                (foreach double-double-float double-float single-float fixnum bignum ratio))
705               (dd-%atan2 y (coerce x 'double-double-float)))
706              #+double-double
707              (((foreach double-float single-float fixnum bignum ratio)
708                double-double-float)
709               (dd-%atan2 (coerce y 'double-double-float) x))))
710        (number-dispatch ((y number))        (number-dispatch ((y number))
711          (handle-reals %atan y)          (handle-reals %atan y)
712          ((complex)          ((complex)
# Line 732  Line 759 
759           (complex-acosh number)           (complex-acosh number)
760           (coerce (%acosh (coerce number 'double-float))           (coerce (%acosh (coerce number 'double-float))
761                   '(dispatch-type number))))                   '(dispatch-type number))))
762        #+double-double
763        ((double-double-float)
764         (if (< number 1w0)
765             (complex-acosh number)
766             (dd-%acosh number)))
767      ((complex)      ((complex)
768       (complex-acosh number))))       (complex-acosh number))))
769    
# Line 749  Line 781 
781           (complex-atanh number)           (complex-atanh number)
782           (coerce (%atanh (coerce number 'double-float))           (coerce (%atanh (coerce number 'double-float))
783                   '(dispatch-type number))))                   '(dispatch-type number))))
784        #+double-double
785        ((double-double-float)
786         (if (or (> number 1w0)
787                 (< number -1w0))
788             (complex-atanh number)
789             (dd-%atanh (coerce number 'double-double-float))))
790      ((complex)      ((complex)
791       (complex-atanh number))))       (complex-atanh number))))
792    
# Line 818  Line 856 
856    
857  (declaim (inline square))  (declaim (inline square))
858  (defun square (x)  (defun square (x)
859    (declare (double-float x))    (declare (float x))
860    (* x x))    (* x x))
861    
862  ;; 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 829  Line 867 
867  (defun scalb (x n)  (defun scalb (x n)
868    "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
869  underlying floating-point format"  underlying floating-point format"
870    (declare (type double-float x)    (declare (type float x)
871             (type double-float-exponent n))             (type double-float-exponent n))
872    (scale-float x n))    (scale-float x n))
873    
# Line 837  underlying floating-point format" Line 875  underlying floating-point format"
875  (defun logb-finite (x)  (defun logb-finite (x)
876    "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
877  that we can always return an integer"  that we can always return an integer"
878    (declare (type double-float x))    (declare (type float x))
879    (multiple-value-bind (signif expon sign)    (multiple-value-bind (signif expon sign)
880        (decode-float x)        (decode-float x)
881      (declare (ignore signif sign))      (declare (ignore signif sign))
# Line 854  For the special cases, the following val Line 892  For the special cases, the following val
892     +/- infinity   +infinity     +/- infinity   +infinity
893     0              -infinity     0              -infinity
894  "  "
895    (declare (type double-float x))    (declare (type float x))
896    (cond ((float-nan-p x)    (cond ((float-nan-p x)
897           x)           x)
898          ((float-infinity-p x)          ((float-infinity-p x)
# Line 862  For the special cases, the following val Line 900  For the special cases, the following val
900          ((zerop x)          ((zerop x)
901           ;; The answer is negative infinity, but we are supposed to           ;; The answer is negative infinity, but we are supposed to
902           ;; signal divide-by-zero, so do the actual division           ;; signal divide-by-zero, so do the actual division
903           (/ -1.0d0 x)           (/ -1 x)
904           )           )
905          (t          (t
906           (logb-finite x))))           (logb-finite x))))
# Line 929  and Y are coerced to single-float." Line 967  and Y are coerced to single-float."
967    
968  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
969    (declare (number z))    (declare (number z))
970      #+double-double
971      (when (typep z '(or double-double-float (complex double-double-float)))
972        (return-from complex-sqrt (dd-complex-sqrt z)))
973    (multiple-value-bind (rho k)    (multiple-value-bind (rho k)
974        (cssqs z)        (cssqs z)
975      (declare (type (or (member 0d0) (double-float 0d0)) rho)      (declare (type (or (member 0d0) (double-float 0d0)) rho)
# Line 1005  This is for use with J /= 0 only when |z Line 1046  This is for use with J /= 0 only when |z
1046    
1047  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1048    (declare (number z))    (declare (number z))
1049      #+double-double
1050      (when (typep z '(or double-double-float (complex double-double-float)))
1051        (return-from complex-log (dd-complex-log-scaled z 0)))
1052    (complex-log-scaled z 0))    (complex-log-scaled z 0))
1053    
1054  ;; Let us note the following "strange" behavior.  atanh 1.0d0 is  ;; Let us note the following "strange" behavior.  atanh 1.0d0 is
# Line 1015  Z may be any number, but the result is a Line 1059  Z may be any number, but the result is a
1059  (defun complex-atanh (z)  (defun complex-atanh (z)
1060    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
1061    (declare (number z))    (declare (number z))
1062      #+double-double
1063      (when (typep z '(or double-double-float (complex double-double-float)))
1064        (return-from complex-atanh (dd-complex-atanh z)))
1065    
1066    (if (and (realp z) (< z -1))    (if (and (realp z) (< z -1))
1067        ;; atanh is continuous in quadrant III in this case.        ;; atanh is continuous in quadrant III in this case.
1068        (complex-atanh (complex z -0f0))        (complex-atanh (complex z -0f0))
# Line 1075  Z may be any number, but the result is a Line 1123  Z may be any number, but the result is a
1123  (defun complex-tanh (z)  (defun complex-tanh (z)
1124    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"
1125    (declare (number z))    (declare (number z))
1126      #+double-double
1127      (when (typep z '(or double-double-float (complex double-double-float)))
1128        (return-from complex-tanh (dd-complex-tanh z)))
1129    
1130    (let ((x (float (realpart z) 1.0d0))    (let ((x (float (realpart z) 1.0d0))
1131          (y (float (imagpart z) 1.0d0)))          (y (float (imagpart z) 1.0d0)))
1132      (locally      (locally
# Line 1158  Z may be any number, but the result is a Line 1210  Z may be any number, but the result is a
1210    
1211  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1212    (declare (number z))    (declare (number z))
1213      #+double-double
1214      (when (typep z '(or double-double-float (complex double-double-float)))
1215        (return-from complex-acos (dd-complex-acos z)))
1216    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1217        ;; acos is continuous in quadrant IV in this case.        ;; acos is continuous in quadrant IV in this case.
1218        (complex-acos (complex z -0f0))        (complex-acos (complex z -0f0))
# Line 1188  Z may be any number, but the result is a Line 1243  Z may be any number, but the result is a
1243    
1244  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1245    (declare (number z))    (declare (number z))
1246      #+double-double
1247      (when (typep z '(or double-double-float (complex double-double-float)))
1248        (return-from complex-asin (dd-complex-asin z)))
1249    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1250        ;; asin is continuous in quadrant IV in this case.        ;; asin is continuous in quadrant IV in this case.
1251        (complex-asin (complex z -0f0))        (complex-asin (complex z -0f0))
# Line 1205  Z may be any number, but the result is a Line 1263  Z may be any number, but the result is a
1263  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1264    (declare (number z))    (declare (number z))
1265    ;; asinh z = -i * asin (i*z)    ;; asinh z = -i * asin (i*z)
1266      #+double-double
1267      (when (typep z '(or double-double-float (complex double-double-float)))
1268        (return-from complex-asinh (dd-complex-asinh z)))
1269    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1270           (result (complex-asin iz)))           (result (complex-asin iz)))
1271      (complex (imagpart result)      (complex (imagpart result)
# Line 1216  Z may be any number, but the result is a Line 1277  Z may be any number, but the result is a
1277  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1278    (declare (number z))    (declare (number z))
1279    ;; atan z = -i * atanh (i*z)    ;; atan z = -i * atanh (i*z)
1280      #+double-double
1281      (when (typep z '(or double-double-float (complex double-double-float)))
1282        (return-from complex-atan (dd-complex-atan z)))
1283    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1284           (result (complex-atanh iz)))           (result (complex-atanh iz)))
1285      (complex (imagpart result)      (complex (imagpart result)
# Line 1227  Z may be any number, but the result is a Line 1291  Z may be any number, but the result is a
1291  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1292    (declare (number z))    (declare (number z))
1293    ;; tan z = -i * tanh(i*z)    ;; tan z = -i * tanh(i*z)
1294      #+double-double
1295      (when (typep z '(or double-double-float (complex double-double-float)))
1296        (return-from complex-tan (dd-complex-tan z)))
1297    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1298           (result (complex-tanh iz)))           (result (complex-tanh iz)))
1299      (complex (imagpart result)      (complex (imagpart result)

Legend:
Removed from v.1.45.2.1.2.1.2.1  
changed lines
  Added in v.1.45.2.1.2.1.2.2

  ViewVC Help
Powered by ViewVC 1.1.5