/[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 by rtoy, Sun Jun 11 20:11:45 2006 UTC revision 1.45.2.1.2.1.2.3 by rtoy, Thu Jun 29 14:50:27 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)))))
378          (((foreach fixnum (or bignum ratio) single-float double-float)          (((foreach fixnum (or bignum ratio) single-float double-float
379                       #+double-double double-double-float)
380            complex)            complex)
381           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
382               (* base power)               (* base power)
383               (exp (* power (log base)))))               (exp (* power (log base)))))
384          (((foreach (complex float) (complex rational))          (((foreach (complex float) (complex rational))
385            (foreach complex double-float single-float))            (foreach complex double-float single-float #+double-double double-double-float))
386           (if (and (zerop base) (plusp (realpart power)))           (if (and (zerop base) (plusp (realpart power)))
387               (* base power)               (* base power)
388               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))
# Line 428  Line 434 
434                 (((foreach single-float fixnum bignum ratio)                 (((foreach single-float fixnum bignum ratio)
435                   double-float)                   double-float)
436                  (/ (log (coerce number 'double-float)) (log base)))                  (/ (log (coerce number 'double-float)) (log base)))
437                   #+double-double
438                   ((double-double-float
439                     (foreach double-double-float double-float single-float fixnum bignum ratio))
440                    (/ (log number) (log (coerce base 'double-double-float))))
441                   #+double-double
442                   (((foreach double-float single-float fixnum bignum ratio)
443                     double-double-float)
444                    (/ (log (coerce number 'double-double-float)) (log base)))
445                 (((foreach single-float fixnum bignum ratio)                 (((foreach single-float fixnum bignum ratio)
446                   (foreach single-float fixnum bignum ratio))                   (foreach single-float fixnum bignum ratio))
447                  ;; Converting everything to double-float helps the                  ;; Converting everything to double-float helps the
# Line 489  Line 503 
503                       '(dispatch-type number))))                       '(dispatch-type number))))
504          #+double-double          #+double-double
505          ((double-double-float)          ((double-double-float)
          ;; Hack!  
506           (let ((hi (kernel:double-double-hi number)))           (let ((hi (kernel:double-double-hi number)))
507             (if (< (float-sign hi)             (if (< (float-sign hi) 0d0)
508                    (coerce 0 '(dispatch-type number)))                 (complex (dd-%log (- number)) dd-pi)
509                 (complex (coerce (log (- hi)) 'kernel:double-double-float)                 (dd-%log number))))
                         (coerce pi '(dispatch-type number)))  
                (coerce (%log hi) '(dispatch-type number)))))  
510          ((complex)          ((complex)
511           (complex-log number)))))           (complex-log number)))))
512    
# Line 513  Line 524 
524                   '(dispatch-type number))))                   '(dispatch-type number))))
525      #+double-double      #+double-double
526      ((double-double-float)      ((double-double-float)
527       (multiple-value-bind (hi lo)       (if (minusp number)
528           (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))           (dd-complex-sqrt number)
529         (kernel:make-double-double-float hi lo)))           (multiple-value-bind (hi lo)
530                 (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
531               (kernel:make-double-double-float hi lo))))
532      ((complex)      ((complex)
533       (complex-sqrt number))))       (complex-sqrt number))))
534    
# Line 581  Line 594 
594      #+double-double      #+double-double
595      (double-double-float      (double-double-float
596       (if (minusp (float-sign number))       (if (minusp (float-sign number))
597           (coerce pi 'double-double-float)           dd-pi
598           (kernel:make-double-double-float 0d0 0d0)))           0w0))
599      (complex      (complex
600       (atan (imagpart number) (realpart number)))))       (atan (imagpart number) (realpart number)))))
601    
# Line 633  Line 646 
646                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
647           (coerce (%asin (coerce number 'double-float))           (coerce (%asin (coerce number 'double-float))
648                   '(dispatch-type number))                   '(dispatch-type number))
649                   (complex-asin number)))           (complex-asin number)))
650        #+double-double
651        ((double-double-float)
652         (if (or (float-nan-p number)
653                 (and (<= number 1w0)
654                      (>= number -1w0)))
655             (dd-%asin number)
656             (dd-complex-asin number)))
657      ((complex)      ((complex)
658       (complex-asin number))))       (complex-asin number))))
659    
# Line 650  Line 670 
670                    (>= number (coerce -1 '(dispatch-type number)))))                    (>= number (coerce -1 '(dispatch-type number)))))
671           (coerce (%acos (coerce number 'double-float))           (coerce (%acos (coerce number 'double-float))
672                   '(dispatch-type number))                   '(dispatch-type number))
673                   (complex-acos number)))           (complex-acos number)))
674        #+double-double
675        ((double-double-float)
676         (if (or (float-nan-p number)
677                 (and (<= number 1w0)
678                      (>= number -1w0)))
679             (dd-%acos number)
680             (complex-acos number)))
681      ((complex)      ((complex)
682       (complex-acos number))))       (complex-acos number))))
683    
# Line 679  Line 706 
706            (((foreach single-float fixnum bignum ratio)            (((foreach single-float fixnum bignum ratio)
707              (foreach single-float fixnum bignum ratio))              (foreach single-float fixnum bignum ratio))
708             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))             (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
709                     'single-float))))                     'single-float))
710              #+double-double
711              ((double-double-float
712                (foreach double-double-float double-float single-float fixnum bignum ratio))
713               (dd-%atan2 y (coerce x 'double-double-float)))
714              #+double-double
715              (((foreach double-float single-float fixnum bignum ratio)
716                double-double-float)
717               (dd-%atan2 (coerce y 'double-double-float) x))))
718        (number-dispatch ((y number))        (number-dispatch ((y number))
719          (handle-reals %atan y)          (handle-reals %atan y)
720          ((complex)          ((complex)
# Line 732  Line 767 
767           (complex-acosh number)           (complex-acosh number)
768           (coerce (%acosh (coerce number 'double-float))           (coerce (%acosh (coerce number 'double-float))
769                   '(dispatch-type number))))                   '(dispatch-type number))))
770        #+double-double
771        ((double-double-float)
772         (if (< number 1w0)
773             (complex-acosh number)
774             (dd-%acosh number)))
775      ((complex)      ((complex)
776       (complex-acosh number))))       (complex-acosh number))))
777    
# Line 749  Line 789 
789           (complex-atanh number)           (complex-atanh number)
790           (coerce (%atanh (coerce number 'double-float))           (coerce (%atanh (coerce number 'double-float))
791                   '(dispatch-type number))))                   '(dispatch-type number))))
792        #+double-double
793        ((double-double-float)
794         (if (or (> number 1w0)
795                 (< number -1w0))
796             (complex-atanh number)
797             (dd-%atanh (coerce number 'double-double-float))))
798      ((complex)      ((complex)
799       (complex-atanh number))))       (complex-atanh number))))
800    
# Line 818  Line 864 
864    
865  (declaim (inline square))  (declaim (inline square))
866  (defun square (x)  (defun square (x)
867    (declare (double-float x))    (declare (float x))
868    (* x x))    (* x x))
869    
870  ;; 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 875 
875  (defun scalb (x n)  (defun scalb (x n)
876    "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
877  underlying floating-point format"  underlying floating-point format"
878    (declare (type double-float x)    (declare (type float x)
879             (type double-float-exponent n))             (type double-float-exponent n))
880    (scale-float x n))    (scale-float x n))
881    
# Line 837  underlying floating-point format" Line 883  underlying floating-point format"
883  (defun logb-finite (x)  (defun logb-finite (x)
884    "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
885  that we can always return an integer"  that we can always return an integer"
886    (declare (type double-float x))    (declare (type float x))
887    (multiple-value-bind (signif expon sign)    (multiple-value-bind (signif expon sign)
888        (decode-float x)        (decode-float x)
889      (declare (ignore signif sign))      (declare (ignore signif sign))
# Line 854  For the special cases, the following val Line 900  For the special cases, the following val
900     +/- infinity   +infinity     +/- infinity   +infinity
901     0              -infinity     0              -infinity
902  "  "
903    (declare (type double-float x))    (declare (type float x))
904    (cond ((float-nan-p x)    (cond ((float-nan-p x)
905           x)           x)
906          ((float-infinity-p x)          ((float-infinity-p x)
# Line 862  For the special cases, the following val Line 908  For the special cases, the following val
908          ((zerop x)          ((zerop x)
909           ;; The answer is negative infinity, but we are supposed to           ;; The answer is negative infinity, but we are supposed to
910           ;; signal divide-by-zero, so do the actual division           ;; signal divide-by-zero, so do the actual division
911           (/ -1.0d0 x)           (/ -1 x)
912           )           )
913          (t          (t
914           (logb-finite x))))           (logb-finite x))))
# Line 929  and Y are coerced to single-float." Line 975  and Y are coerced to single-float."
975    
976  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
977    (declare (number z))    (declare (number z))
978      #+double-double
979      (when (typep z '(or double-double-float (complex double-double-float)))
980        (return-from complex-sqrt (dd-complex-sqrt z)))
981    (multiple-value-bind (rho k)    (multiple-value-bind (rho k)
982        (cssqs z)        (cssqs z)
983      (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 1054  This is for use with J /= 0 only when |z
1054    
1055  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1056    (declare (number z))    (declare (number z))
1057      #+double-double
1058      (when (typep z '(or double-double-float (complex double-double-float)))
1059        (return-from complex-log (dd-complex-log-scaled z 0)))
1060    (complex-log-scaled z 0))    (complex-log-scaled z 0))
1061    
1062  ;; 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 1067  Z may be any number, but the result is a
1067  (defun complex-atanh (z)  (defun complex-atanh (z)
1068    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
1069    (declare (number z))    (declare (number z))
1070      #+double-double
1071      (when (typep z '(or double-double-float (complex double-double-float)))
1072        (return-from complex-atanh (dd-complex-atanh z)))
1073    
1074    (if (and (realp z) (< z -1))    (if (and (realp z) (< z -1))
1075        ;; atanh is continuous in quadrant III in this case.        ;; atanh is continuous in quadrant III in this case.
1076        (complex-atanh (complex z -0f0))        (complex-atanh (complex z -0f0))
# Line 1075  Z may be any number, but the result is a Line 1131  Z may be any number, but the result is a
1131  (defun complex-tanh (z)  (defun complex-tanh (z)
1132    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"
1133    (declare (number z))    (declare (number z))
1134      #+double-double
1135      (when (typep z '(or double-double-float (complex double-double-float)))
1136        (return-from complex-tanh (dd-complex-tanh z)))
1137    
1138    (let ((x (float (realpart z) 1.0d0))    (let ((x (float (realpart z) 1.0d0))
1139          (y (float (imagpart z) 1.0d0)))          (y (float (imagpart z) 1.0d0)))
1140      (locally      (locally
# Line 1158  Z may be any number, but the result is a Line 1218  Z may be any number, but the result is a
1218    
1219  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1220    (declare (number z))    (declare (number z))
1221      #+double-double
1222      (when (typep z '(or double-double-float (complex double-double-float)))
1223        (return-from complex-acos (dd-complex-acos z)))
1224    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1225        ;; acos is continuous in quadrant IV in this case.        ;; acos is continuous in quadrant IV in this case.
1226        (complex-acos (complex z -0f0))        (complex-acos (complex z -0f0))
# Line 1188  Z may be any number, but the result is a Line 1251  Z may be any number, but the result is a
1251    
1252  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1253    (declare (number z))    (declare (number z))
1254      #+double-double
1255      (when (typep z '(or double-double-float (complex double-double-float)))
1256        (return-from complex-asin (dd-complex-asin z)))
1257    (if (and (realp z) (> z 1))    (if (and (realp z) (> z 1))
1258        ;; asin is continuous in quadrant IV in this case.        ;; asin is continuous in quadrant IV in this case.
1259        (complex-asin (complex z -0f0))        (complex-asin (complex z -0f0))
# Line 1205  Z may be any number, but the result is a Line 1271  Z may be any number, but the result is a
1271  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1272    (declare (number z))    (declare (number z))
1273    ;; asinh z = -i * asin (i*z)    ;; asinh z = -i * asin (i*z)
1274      #+double-double
1275      (when (typep z '(or double-double-float (complex double-double-float)))
1276        (return-from complex-asinh (dd-complex-asinh z)))
1277    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1278           (result (complex-asin iz)))           (result (complex-asin iz)))
1279      (complex (imagpart result)      (complex (imagpart result)
# Line 1216  Z may be any number, but the result is a Line 1285  Z may be any number, but the result is a
1285  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1286    (declare (number z))    (declare (number z))
1287    ;; atan z = -i * atanh (i*z)    ;; atan z = -i * atanh (i*z)
1288      #+double-double
1289      (when (typep z '(or double-double-float (complex double-double-float)))
1290        (return-from complex-atan (dd-complex-atan z)))
1291    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1292           (result (complex-atanh iz)))           (result (complex-atanh iz)))
1293      (complex (imagpart result)      (complex (imagpart result)
# Line 1227  Z may be any number, but the result is a Line 1299  Z may be any number, but the result is a
1299  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1300    (declare (number z))    (declare (number z))
1301    ;; tan z = -i * tanh(i*z)    ;; tan z = -i * tanh(i*z)
1302      #+double-double
1303      (when (typep z '(or double-double-float (complex double-double-float)))
1304        (return-from complex-tan (dd-complex-tan z)))
1305    (let* ((iz (complex (- (imagpart z)) (realpart z)))    (let* ((iz (complex (- (imagpart z)) (realpart z)))
1306           (result (complex-tanh iz)))           (result (complex-tanh iz)))
1307      (complex (imagpart result)      (complex (imagpart result)

Legend:
Removed from v.1.45.2.1.2.1  
changed lines
  Added in v.1.45.2.1.2.1.2.3

  ViewVC Help
Powered by ViewVC 1.1.5