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

Diff of /src/code/irrat-dd.lisp

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

revision 1.1.2.5 by rtoy, Thu Jun 29 20:49:12 2006 UTC revision 1.20 by rtoy, Tue Apr 20 17:57:44 2010 UTC
# Line 19  Line 19 
19    
20  (in-package "KERNEL")  (in-package "KERNEL")
21    
22  ;;;; Random constants, utility functions.  (intl:textdomain "cmucl")
   
 #+nil  
 (defconstant max-log  
   1.1356523406294143949491931077970764891253w4  
   "log(most-positive-double-double-float)")  
23    
24  #+nil  ;;;; Random constants, utility functions.
 (defconstant min-log  
   -1.143276959615573793352782661133116431383730w4  
   "log(least-positive-double-double-float")  
25    
26  (defconstant max-log  (defconstant max-log
27    7.0978271289338399678773454114191w2    7.0978271289338399678773454114191w2
28    "log(most-positive-double-double-float)")    _N"log(most-positive-double-double-float)")
29    
30  (defconstant min-log  (defconstant min-log
31    -7.4444007192138126231410729844608w2    -7.4444007192138126231410729844608w2
32    "log(least-positive-double-double-float")    _N"log(least-positive-double-double-float")
33    
34    
35  (defconstant loge2  (defconstant loge2
36    0.6931471805599453094172321214581765680755001w0    0.6931471805599453094172321214581765680755001w0
37    "log(2)")    _N"log(2)")
38    
39  (defconstant log2e  (defconstant log2e
40    1.442695040888963407359924681001892137426646w0    1.442695040888963407359924681001892137426646w0
41    "Log base 2 of e")    _N"Log base 2 of e")
42    
43  (defconstant log2ea  (defconstant log2ea
44    4.4269504088896340735992468100189213742664595w-1    4.4269504088896340735992468100189213742664595w-1
45    "log2(e)-1")    _N"log2(e)-1")
46    
47  (defconstant dd-pi  (defconstant dd-pi
48    3.141592653589793238462643383279502884197169w0    3.141592653589793238462643383279502884197169w0
49    "Pi")    _N"Pi")
50    
51  (defconstant dd-pi/2  (defconstant dd-pi/2
52    1.570796326794896619231321691639751442098585w0    1.570796326794896619231321691639751442098585w0
53    "Pi/2")    _N"Pi/2")
54    
55  (defconstant dd-pi/4  (defconstant dd-pi/4
56    0.7853981633974483096156608458198757210492923w0    0.7853981633974483096156608458198757210492923w0
57    "Pi/4")    _N"Pi/4")
58    
59  ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +  ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +
60  ;; log2-c2 is log(2) to an accuracy greater than double-double-float.  ;; log2-c2 is log(2) to an accuracy greater than double-double-float.
# Line 74  Line 66 
66    
67  (defconstant sqrt-1/2  (defconstant sqrt-1/2
68    0.7071067811865475244008443621048490392848w0    0.7071067811865475244008443621048490392848w0
69    "Sqrt(2)")    _N"Sqrt(1/2)")
70    
71  ;; Evaluate polynomial  ;; Evaluate polynomial
72    (declaim (maybe-inline poly-eval poly-eval-1))
73  (defun poly-eval (x coef)  (defun poly-eval (x coef)
74    (declare (type double-double-float x)    (declare (type double-double-float x)
75             (type (simple-array double-double-float (*)) coef))             (type (simple-array double-double-float (*)) coef)
76               (optimize (speed 3) (space 0) (inhibit-warnings 3)))
77    ;; y = C0 + C1*x + C2*x^2 + ...    ;; y = C0 + C1*x + C2*x^2 + ...
78    ;;    ;;
79    ;; But coefficients are stored in reverse (descending powers) order:    ;; But coefficients are stored in reverse (descending powers) order:
# Line 93  Line 87 
87    
88  (defun poly-eval-1 (x coef)  (defun poly-eval-1 (x coef)
89    (declare (type double-double-float x)    (declare (type double-double-float x)
90             (type (simple-array double-double-float (*)) coef))             (type (simple-array double-double-float (*)) coef)
91               (optimize (speed 3) (space 0) (inhibit-warnings 3)))
92    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.
93    (let ((y 1w0))    (let ((y 1w0))
94      (declare (type double-double-float y))      (declare (type double-double-float y))
# Line 136  Line 131 
131        (q (make-array 8 :element-type 'double-double-float        (q (make-array 8 :element-type 'double-double-float
132                       :initial-contents '(                       :initial-contents '(
133                                           ;; 1.000000000000000000000000000000000000000w0                                           ;; 1.000000000000000000000000000000000000000w0
134                                           1.766112549341972444333352727998584753865w9                                           -8.802340681794263968892934703309274564037w1
                                          -7.848989743695296475743081255027098295771w8  
                                          1.615869009634292424463780387327037251069w8  
                                          -2.019684072836541751428967854947019415698w7  
                                          1.682912729190313538934190635536631941751w6  
                                          -9.615511549171441430850103489315371768998w4  
135                                           3.697714952261803935521187272204485251835w3                                           3.697714952261803935521187272204485251835w3
136                                           -8.802340681794263968892934703309274564037w1)))                                           -9.615511549171441430850103489315371768998w4
137                                             1.682912729190313538934190635536631941751w6
138                                             -2.019684072836541751428967854947019415698w7
139                                             1.615869009634292424463780387327037251069w8
140                                             -7.848989743695296475743081255027098295771w8
141                                             1.766112549341972444333352727998584753865w9)))
142        ;; ln 2^-114        ;; ln 2^-114
143        (minarg -7.9018778583833765273564461846232128760607w1))        (minarg -7.9018778583833765273564461846232128760607w1))
144    (declare (type double-double-float minarg))    (declare (type double-double-float minarg))
# Line 151  Line 146 
146    ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn    ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn
147    ;; above.  Theoretical peak relative error = 8.1e-36.    ;; above.  Theoretical peak relative error = 8.1e-36.
148    (defun dd-%expm1 (x)    (defun dd-%expm1 (x)
149      "exp(x) - 1"      _N"exp(x) - 1"
150      (declare (type double-double-float x)      (declare (type double-double-float x)
151               (optimize (speed 3)))               (optimize (speed 3) (space 0)
152                           (inhibit-warnings 3)))
153      ;; Range reduction is accomplished by separating the argument      ;; Range reduction is accomplished by separating the argument
154      ;; into an integer k and fraction f such that      ;; into an integer k and fraction f such that
155      ;;      ;;
# Line 171  Line 167 
167    
168      ;; Express x = ln(2)*(k+remainder), remainder not exceeding 1/2      ;; Express x = ln(2)*(k+remainder), remainder not exceeding 1/2
169      (let* ((xx (+ log2-c1 log2-c2))      (let* ((xx (+ log2-c1 log2-c2))
170             (k (floor (+ 1/2 (/ (the (double-double-float * 710w0) x) xx))))             (k (floor (+ 0.5w0 (/ (the (double-double-float * 710w0) x) xx))))
171             (px (coerce k 'double-double-float))             (px (coerce k 'double-double-float))
172             (qx 0w0))             (qx 0w0))
173        (declare (type double-double-float xx px qx))        (declare (type double-double-float xx px qx))
# Line 180  Line 176 
176        (decf x (* px log2-c2))        (decf x (* px log2-c2))
177    
178        ;; Approximate exp(remainder*ln(2))        ;; Approximate exp(remainder*ln(2))
179        #+nil        (setf px (* x (poly-eval x p)))
       (setf px (* x  
                   (+ p0  
                      (* x  
                         (+ p1  
                            (* x  
                               (+ p2  
                                  (* x  
                                     (+ p3  
                                        (* x  
                                           (+ p4  
                                              (* x  
                                                 (+ p5  
                                                    (* x  
                                                       (+ p6  
                                                          (* x p7))))))))))))))))  
       #+nil  
       (setf qx (+ q0  
                   (* x  
                      (+ q1  
                         (* x  
                            (+ q2  
                               (* x  
                                  (+ q3  
                                     (* x  
                                        (+ q4  
                                           (* x  
                                              (+ q5  
                                                 (* x  
                                                    (+ q6  
                                                       (* x  
                                                          (+ x q7))))))))))))))))  
       (setf px (poly-eval x p))  
180        (setf qx (poly-eval-1 x q))        (setf qx (poly-eval-1 x q))
181        (setf xx (* x x))        (setf xx (* x x))
182        (setf qx (+ x (+ (* 0.5w0 xx)        (setf qx (+ x (+ (* 0.5w0 xx)
# Line 282  Line 246 
246    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.
247    (defun dd-%exp (x)    (defun dd-%exp (x)
248      (declare (type double-double-float x)      (declare (type double-double-float x)
249               (optimize (speed 3)))               (optimize (speed 3) (space 0)
250                           (inhibit-warnings 3)))
251      (when (> x max-log)      (when (> x max-log)
252        (return-from dd-%exp        (return-from dd-%exp
253          (kernel:%make-double-double-float #.ext:double-float-positive-infinity          (kernel:%make-double-double-float #.ext:double-float-positive-infinity
# Line 394  Line 359 
359    ;;    ;;
360    (defun dd-%log (x)    (defun dd-%log (x)
361      (declare (type double-double-float x)      (declare (type double-double-float x)
362               (optimize (speed 3)))               (optimize (speed 3) (space 0)
363                           (inhibit-warnings 3)))
364      ;; Separate mantissa from exponent      ;; Separate mantissa from exponent
365      (multiple-value-bind (x e)      (multiple-value-bind (x e)
366          (decode-float x)          (decode-float x)
# Line 517  Line 483 
483    ;; Theoretical peak relative error = 5.3e-37,    ;; Theoretical peak relative error = 5.3e-37,
484    ;; relative peak error spread = 2.3e-14    ;; relative peak error spread = 2.3e-14
485    (defun dd-%log1p (xm1)    (defun dd-%log1p (xm1)
486      (declare (type double-double-float xm1))      (declare (type double-double-float xm1)
487                 (optimize (space 0)))
488      (let ((x (+ xm1 1))      (let ((x (+ xm1 1))
489            (z 0w0)            (z 0w0)
490            (y 0w0))            (y 0w0))
491        (declare (type double-double-float x y z)        (declare (type double-double-float x y z)
492                 (optimize (speed 3)))                 (optimize (speed 3)
493                             (inhibit-warnings 3)))
494        (multiple-value-bind (x e)        (multiple-value-bind (x e)
495            (decode-float x)            (decode-float x)
496          (declare (type double-double-float x)          (declare (type double-double-float x)
# Line 584  Line 552 
552                                           3.889701475261939961851358705515223019890w14))))                                           3.889701475261939961851358705515223019890w14))))
553    (defun dd-%sinh (x)    (defun dd-%sinh (x)
554      (declare (type double-double-float x)      (declare (type double-double-float x)
555               (optimize (speed 3)))               (optimize (speed 3) (space 0)
556                           (inhibit-warnings 3)))
557      (let ((a (abs x)))      (let ((a (abs x)))
558        (declare (type double-double-float a))        (declare (type double-double-float a))
559        (cond ((> a 1)        (cond ((> a 1)
# Line 618  Line 587 
587               (* 0.5w0 (+ d (/ d (+ 1 d)))))))))               (* 0.5w0 (+ d (/ d (+ 1 d)))))))))
588    
589  (defun dd-%cosh (x)  (defun dd-%cosh (x)
590    (declare (type double-double-float x))    (declare (type double-double-float x)
591               (optimize (speed 3) (space 0)
592                         (inhibit-warnings 3)))
593    (let ((y (dd-%exp x)))    (let ((y (dd-%exp x)))
594      (scale-float (+ y (/ y)) -1)))      (scale-float (+ y (/ y)) -1)))
595    
# Line 640  Line 611 
611                                           1.365413143842835040443257277862054198329w8))))                                           1.365413143842835040443257277862054198329w8))))
612    (defun dd-%tanh (x)    (defun dd-%tanh (x)
613      (declare (type double-double-float x)      (declare (type double-double-float x)
614               (optimize (speed 3)))               (optimize (speed 3) (space 0)
615                           (inhibit-warnings 3)))
616      (let ((z (abs x)))      (let ((z (abs x)))
617        (declare (type double-double-float z))        (declare (type double-double-float z))
618        (cond ((> z (* 0.5w0 max-log))        (cond ((> z (* 0.5w0 max-log))
# Line 690  Line 662 
662                                           ))))                                           ))))
663    (defun dd-%atanh (x)    (defun dd-%atanh (x)
664      (declare (type double-double-float x)      (declare (type double-double-float x)
665               (optimize (speed 3)))               (optimize (speed 3) (space 0)
666                           (inhibit-warnings 3)))
667      (cond ((minusp x)      (cond ((minusp x)
668             (- (the double-double-float (dd-%atanh (- x)))))             (- (the double-double-float (dd-%atanh (- x)))))
669            ((< x 1w-12)            ((< x 1w-12)
# Line 730  Line 703 
703                                           6.226145049170155830806967678679167550122w4))))                                           6.226145049170155830806967678679167550122w4))))
704    (defun dd-%asinh (x)    (defun dd-%asinh (x)
705      (declare (type double-double-float x)      (declare (type double-double-float x)
706               (optimize (speed 3)))               (optimize (speed 3) (space 0)
707                           (inhibit-warnings 3)))
708      (cond ((minusp x)      (cond ((minusp x)
709             (- (the double-double-float (dd-%asinh (- x)))))             (- (the double-double-float (dd-%asinh (- x)))))
710            #+nil            #+nil
# Line 771  Line 745 
745                                           3.190482951215438078279772140481195226621w7                                           3.190482951215438078279772140481195226621w7
746                                           ))))                                           ))))
747    (defun dd-%acosh (x)    (defun dd-%acosh (x)
748      (declare (type double-double-float x)      (declare (type (double-double-float 1w0) x)
749               (optimize (speed 3)))               (optimize (speed 3) (space 0)
750                           (inhibit-warnings 3)))
751      (cond ((> x 1w17)      (cond ((> x 1w17)
752             (+ loge2 (dd-%log x)))             (+ loge2 (dd-%log x)))
753            (t            (t
# Line 784  Line 759 
759                     (t                     (t
760                      (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))                      (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))
761    
762    
763  (let ((P (make-array 10 :element-type 'double-double-float  (let ((P (make-array 10 :element-type 'double-double-float
764                       :initial-contents '(                       :initial-contents '(
765                                           -8.067112765482705313585175280952515549833w-1                                           -8.067112765482705313585175280952515549833w-1
# Line 812  Line 788 
788                                           7.122795760168575261226395089432959614179w4                                           7.122795760168575261226395089432959614179w4
789                                           ))))                                           ))))
790    (defun dd-%asin (x)    (defun dd-%asin (x)
791      (declare (type double-double-float x)      (declare (type (double-double-float -1w0 1w0) x)
792               (optimize (speed 3)))               (optimize (speed 3) (space 0)
793                           (inhibit-warnings 3)))
794      (cond ((minusp x)      (cond ((minusp x)
795             (- (the double-double-float (dd-%asin (- x)))))             (- (the double-double-float (dd-%asin (- x)))))
796            #+nil            #+nil
797            ((< x 1w-8)            ((< x 1w-8)
798             x)             x)
799            (t            (t
800             (let ((flag nil)             (multiple-value-bind (flag z zz)
801                   (z 0w0)                 (cond ((> x 0.5w0)
802                   (zz 0w0))                        (let* ((zz1 (- 0.5w0 x))
803                                 (zz (scale-float (+ zz1 0.5w0) -1)))
804                            (values t (sqrt zz) zz)))
805                         (t
806                          (values nil x (* x x))))
807               (declare (type double-double-float z zz))               (declare (type double-double-float z zz))
808               (cond ((> x 0.5w0)  
                     (setf zz (- 0.5w0 x))  
                     (setf zz (scale-float (+ zz 0.5w0) -1))  
                     (setf z (sqrt zz))  
                     (setf flag t))  
                    (t  
                     (setf z x)  
                     (setf zz (* z z))  
                     (setf flag nil)))  
809               (let ((p (* zz (/ (poly-eval zz p)               (let ((p (* zz (/ (poly-eval zz p)
810                                 (poly-eval-1 zz q)))))                                 (poly-eval-1 zz q)))))
811                 (setf z (+ (* z p) z))                 (setf z (+ (* z p) z))
# Line 842  Line 815 
815                 z))))))                 z))))))
816    
817  (defun dd-%acos (x)  (defun dd-%acos (x)
818    (declare (type double-double-float x)    (declare (type (double-double-float -1w0 1w0) x)
819             (optimize (speed 3)))             (optimize (speed 3) (space 0)
820                         (inhibit-warnings 3)))
821    (cond ((< x -0.5w0)    (cond ((< x -0.5w0)
822           (- dd-pi           (- dd-pi
823              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))
# Line 884  Line 858 
858    (declare (type double-double-float t3p8 tp8))    (declare (type double-double-float t3p8 tp8))
859    (defun dd-%atan (x)    (defun dd-%atan (x)
860      (declare (type double-double-float x)      (declare (type double-double-float x)
861               (optimize (speed 3)))               (optimize (speed 3) (space 0)
862                           (inhibit-warnings 3)))
863      (when (minusp x)      (when (minusp x)
864        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))
865      ;; range reduction      ;; range reduction
# Line 910  Line 885 
885    
886  (defun dd-%atan2 (y x)  (defun dd-%atan2 (y x)
887    (declare (type double-double-float x y)    (declare (type double-double-float x y)
888             (optimize (speed 3)))             (optimize (speed 3) (space 0)
889                         (inhibit-warnings 3)))
890    (let ((code 0)    (let ((code 0)
891          (w 0w0))          (neg-x (minusp (float-sign x)))
892      (declare (type (integer 0 3) code)          (neg-y (minusp (float-sign y))))
893               (type double-double-float w))      (declare (type (integer 0 3) code))
894      (when (minusp x)      (when neg-x
895        (setf code 2))        (setf code 2))
896      (when (minusp y)      (when neg-y
897        (setf code (logior code 1)))        (setf code (logior code 1)))
898      (when (zerop x)      ;; Handle the case where x or y is zero, taking into account the
899        (unless (zerop (logand code 1))      ;; sign of the zero.
900          (return-from dd-%atan2 (- dd-pi/2)))      (cond ((zerop x)
901        (when (zerop y)             ;; x = 0.
902          (return-from dd-%atan2 0w0))             (cond ((zerop y)
903        (return-from dd-%atan2 dd-pi/2))                    ;; x = 0, y = 0
904      (when (zerop y)                    ;;
905        (return-from dd-%atan2                    ;; CLHS Figure 12-15 shows what the results should
906          (if (zerop (logand code 2))                    ;; be.
907              0w0                    (if neg-x
908              dd-pi)))                        (if neg-y
909      (setf w (ecase code                            (- dd-pi)
910                (0 0w0)                            dd-pi)
911                (1 0w0)                        (if neg-y
912                (2 dd-pi)                            -0w0
913                (3 (- dd-pi))))                            0w0)))
914                     (t
915                      ;; x = 0, y /= 0
916                      (if neg-y
917                          (- dd-pi/2)
918                          dd-pi/2))))
919              ((zerop y)
920               ;; y = 0, x /= 0
921               (cond (neg-x
922                      (if neg-y
923                          (- dd-pi)
924                          dd-pi))
925                     (t
926                      (if neg-y
927                          -0w0
928                          0w0))))
929              (t
930               (let ((w (ecase code
931                          (0
932                           ;; x > 0, y > 0
933                           0w0)
934                          (1
935                           ;; x > 0, y < 0
936                           -0w0)
937                          (2
938                           ;; x < 0, y > 0
939                           dd-pi)
940                          (3
941                           ;; x < 0, y < 0
942                           (- dd-pi)))))
943    
944      (+ w (dd-%atan (/ y x)))))               (+ w (dd-%atan (/ y x))))))))
945    
946    
947  ;;  ;;
# Line 1023  pi/4    11001001000011111101101010100010 Line 1028  pi/4    11001001000011111101101010100010
1028  (defconstant dp4  (defconstant dp4
1029    (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))    (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
1030    
1031  (defun dd-%sin (x)  (defun dd-%%sin (x)
1032    (declare (type double-double-float x))    (declare (type double-double-float x)
1033               (optimize (speed 2) (space 0)
1034                         (inhibit-warnings 3)))
1035    (when (minusp x)    (when (minusp x)
1036      (return-from dd-%sin (- (the double-double-float (dd-%sin (- x))))))      (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
1037    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
1038    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1039           (z (scale-float y -4)))           (z (scale-float y -4)))
# Line 1036  pi/4    11001001000011111101101010100010 Line 1043  pi/4    11001001000011111101101010100010
1043    
1044      (let ((j (truncate z))      (let ((j (truncate z))
1045            (sign 1))            (sign 1))
1046          (declare (type (integer -1 1) sign))
1047        (unless (zerop (logand j 1))        (unless (zerop (logand j 1))
1048          (incf j)          (incf j)
1049          (incf y))          (incf y))
# Line 1059  pi/4    11001001000011111101101010100010 Line 1067  pi/4    11001001000011111101101010100010
1067              (- y)              (- y)
1068              y)))))              y)))))
1069    
1070  (defun dd-%cos (x)  (defun dd-%%cos (x)
1071    (declare (type double-double-float x)    (declare (type double-double-float x)
1072             (optimize (speed 3)))             (optimize (speed 2) (space 0)
1073                         (inhibit-warnings 3)))
1074    (when (minusp x)    (when (minusp x)
1075      (return-from dd-%cos (dd-%cos (- x))))      (return-from dd-%%cos (dd-%%cos (- x))))
1076    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
1077    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1078           (z (scale-float y -4)))           (z (scale-float y -4)))
# Line 1125  pi/4    11001001000011111101101010100010 Line 1134  pi/4    11001001000011111101101010100010
1134                         ))))                         ))))
1135    (defun dd-tancot (xx cotflag)    (defun dd-tancot (xx cotflag)
1136      (declare (type double-double-float xx)      (declare (type double-double-float xx)
1137               (optimize (speed 3)))               (optimize (speed 2) (space 0)))
1138      (let ((x 0w0)      (let ((x 0w0)
1139            (sign 1))            (sign 1))
1140        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1168  pi/4    11001001000011111101101010100010 Line 1177  pi/4    11001001000011111101101010100010
1177                (- y)                (- y)
1178                y))))))                y))))))
1179    
1180  (defun dd-%tan (x)  (defun dd-%%tan (x)
1181    (declare (type double-double-float x))    (declare (type double-double-float x))
1182    (dd-tancot x nil))    (dd-tancot x nil))
1183    
1184    (declaim (inline %kernel-rem-pi/2))
1185    (alien:def-alien-routine ("__kernel_rem_pio2" %kernel-rem-pi/2) c-call:int
1186      (x (* double-float))
1187      (y (* double-float))
1188      (e0 c-call:int)
1189      (nx c-call:int)
1190      (prec c-call:int)
1191      (ipio2 (* c-call:int)))
1192    
1193    ;; This is taken from two_over_pi in fdlibm's e_rem_pio2.c.  We do
1194    ;; this here so that the Sparc version doesn't have to compile in
1195    ;; e_rem_pio2, which we don't need.  (But x86 and ppc do.)
1196    (defconstant two-over-pi
1197      (make-array 66 :element-type '(unsigned-byte 32)
1198                  :initial-contents
1199                  '(#xA2F983 #x6E4E44 #x1529FC #x2757D1 #xF534DD #xC0DB62
1200                    #x95993C #x439041 #xFE5163 #xABDEBB #xC561B7 #x246E3A
1201                    #x424DD2 #xE00649 #x2EEA09 #xD1921C #xFE1DEB #x1CB129
1202                    #xA73EE8 #x8235F5 #x2EBB44 #x84E99C #x7026B4 #x5F7E41
1203                    #x3991D6 #x398353 #x39F49C #x845F8B #xBDF928 #x3B1FF8
1204                    #x97FFDE #x05980F #xEF2F11 #x8B5A0A #x6D1F6D #x367ECF
1205                    #x27CB09 #xB74F46 #x3F669E #x5FEA2D #x7527BA #xC7EBE5
1206                    #xF17B3D #x0739F7 #x8A5292 #xEA6BFB #x5FB11F #x8D5D08
1207                    #x560330 #x46FC7B #x6BABF0 #xCFBC20 #x9AF436 #x1DA9E3
1208                    #x91615E #xE61B08 #x659985 #x5F14A0 #x68408D #xFFD880
1209                    #x4D7327 #x310606 #x1556CA #x73A8C9 #x60E27B #xC08C6B
1210                    ))
1211      _N"396 (hex) digits of 2/pi")
1212    
1213    
1214    (let ((y (make-array 3 :element-type 'double-float))
1215          (parts (make-array 5 :element-type 'double-float)))
1216      (declare (type (simple-array double-float (3)) y)
1217               (type (simple-array double-float (5)) parts))
1218      ;; Take the double-double-float number and break it into 24-bit
1219      ;; chunks.  Each chunk is an integer, which is coerced to a
1220      ;; double-float and stored in PARTS.
1221      (defun dd-expand (x)
1222        (declare (double-double-float x)
1223                 (optimize (speed 3) (space 0)
1224                           (inhibit-warnings 3)))
1225        (multiple-value-bind (frac exp)
1226            (decode-float x)
1227          (declare (double-double-float frac)
1228                   (type (signed-byte 16) exp))
1229          (setf frac (scale-float frac 24))
1230          (decf exp 24)
1231          (dotimes (k 5)
1232            (setf (aref parts k) (coerce (ffloor frac) 'double-float))
1233            (setf frac (scale-float (- frac (aref parts k)) 24)))
1234          exp))
1235      (defun reduce-arg (x)
1236        (declare (double-double-float x)
1237                 (optimize (speed 3)
1238                           (inhibit-warnings 3)))
1239        (let* ((e0 (dd-expand x))
1240               (n (sys:without-gcing
1241                   (%kernel-rem-pi/2 (vector-sap parts)
1242                                     (vector-sap y)
1243                                     e0
1244                                     (length parts)
1245                                     3
1246                                     (vector-sap two-over-pi))))
1247               (sum (+ (coerce (aref y 2) 'double-double-float)
1248                       (coerce (aref y 1) 'double-double-float)
1249                       (coerce (aref y 0) 'double-double-float))))
1250          (values n sum))))
1251    
1252    
1253    (declaim (ftype (function (double-double-float) double-double-float)
1254                    dd-%sin))
1255    (defun dd-%sin (x)
1256      (declare (double-double-float x))
1257      (cond ((minusp (float-sign x))
1258             (- (dd-%sin (- x))))
1259            ((< (abs x) (/ pi 4))
1260             (dd-%%sin x))
1261            (t
1262             ;; Argument reduction needed
1263             (multiple-value-bind (n reduced)
1264                 (reduce-arg x)
1265               (ecase (logand n 3)
1266                 (0 (dd-%%sin reduced))
1267                 (1 (dd-%%cos reduced))
1268                 (2 (- (dd-%%sin reduced)))
1269                 (3 (- (dd-%%cos reduced))))))))
1270    
1271    (declaim (ftype (function (double-double-float) double-double-float)
1272                    dd-%cos))
1273    (defun dd-%cos (x)
1274      (declare (double-double-float x))
1275      (cond ((minusp x)
1276             (dd-%cos (- x)))
1277            ((< (abs x) (/ pi 4))
1278             (dd-%%cos x))
1279            (t
1280             ;; Argument reduction needed
1281             (multiple-value-bind (n reduced)
1282                 (reduce-arg x)
1283               (ecase (logand n 3)
1284                 (0 (dd-%%cos reduced))
1285                 (1 (- (dd-%%sin reduced)))
1286                 (2 (- (dd-%%cos reduced)))
1287                 (3 (dd-%%sin reduced)))))))
1288    
1289    (declaim (ftype (function (double-double-float) double-double-float)
1290                    dd-%tan))
1291    (defun dd-%tan (x)
1292      (declare (double-double-float x))
1293      (cond ((minusp (float-sign x))
1294             (- (dd-%tan (- x))))
1295            ((< (abs x) (/ pi 4))
1296             (dd-%%tan x))
1297            (t
1298             ;; Argument reduction needed
1299             (multiple-value-bind (n reduced)
1300                 (reduce-arg x)
1301               (if (evenp n)
1302                   (dd-%%tan reduced)
1303                   (- (/ (dd-%%tan reduced))))))))
1304    
1305  ;;; dd-%log2  ;;; dd-%log2
1306  ;;; Base 2 logarithm.  ;;; Base 2 logarithm.
1307    
# Line 1237  pi/4    11001001000011111101101010100010 Line 1367  pi/4    11001001000011111101101010100010
1367                         ))))                         ))))
1368    (defun dd-%log2 (x)    (defun dd-%log2 (x)
1369      (declare (type double-double-float x)      (declare (type double-double-float x)
1370               (optimize (speed 3)))               (optimize (speed 3) (space 0)
1371                           (inhibit-warnings 3)))
1372      (multiple-value-bind (x e)      (multiple-value-bind (x e)
1373          (decode-float x)          (decode-float x)
1374        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1305  pi/4    11001001000011111101101010100010 Line 1436  pi/4    11001001000011111101101010100010
1436                         ))))                         ))))
1437    (defun dd-%exp2 (x)    (defun dd-%exp2 (x)
1438      (declare (type double-double-float x)      (declare (type double-double-float x)
1439               (optimize (speed 3)))               (optimize (speed 3) (space 0)
1440                           (inhibit-warnings 3)))
1441      (when (>= x 1024w0)      (when (>= x 1024w0)
1442        (return-from dd-%exp2        (return-from dd-%exp2
1443          (%make-double-double-float ext:double-float-positive-infinity          (%make-double-double-float ext:double-float-positive-infinity
# Line 1327  pi/4    11001001000011111101101010100010 Line 1459  pi/4    11001001000011111101101010100010
1459  (defun dd-%powil (x nn)  (defun dd-%powil (x nn)
1460    (declare (type double-double-float x)    (declare (type double-double-float x)
1461             (fixnum nn)             (fixnum nn)
1462             (optimize (speed 3)))             (optimize (speed 3) (space 0)
1463                         (inhibit-warnings 3)))
1464    (when (zerop x)    (when (zerop x)
1465      (return-from dd-%powil      (return-from dd-%powil
1466        (cond ((zerop nn)        (cond ((zerop nn)
# Line 1375  pi/4    11001001000011111101101010100010 Line 1508  pi/4    11001001000011111101101010100010
1508                 (setf s (* loge2 e))))                 (setf s (* loge2 e))))
1509          (when (> s max-log)          (when (> s max-log)
1510            ;; Overflow.  What to do?            ;; Overflow.  What to do?
1511            (error "Overflow"))            (error (intl:gettext "Overflow")))
1512          (when (< s min-log)          (when (< s min-log)
1513            (return-from dd-%powil 0w0))            (return-from dd-%powil 0w0))
1514    
# Line 1415  pi/4    11001001000011111101101010100010 Line 1548  pi/4    11001001000011111101101010100010
1548    
1549  (defun dd-real-pow (x y)  (defun dd-real-pow (x y)
1550    (declare (type double-double-float x y)    (declare (type double-double-float x y)
1551             (optimize (speed 3)))             (optimize (speed 3) (space 0)
1552                         (inhibit-warnings 3)))
1553    (let ((nflg 0)    (let ((nflg 0)
1554          (w (floor y)))          (w (floor y)))
1555      ;; nflg = 1 if x < 0 raised to integer power      ;; nflg = 1 if x < 0 raised to integer power
# Line 1471  pi/4    11001001000011111101101010100010 Line 1605  pi/4    11001001000011111101101010100010
1605      ;; catch the overflow or underflow signal?  For now, we turn all      ;; catch the overflow or underflow signal?  For now, we turn all
1606      ;; traps off and look at the accrued exceptions to see if any      ;; traps off and look at the accrued exceptions to see if any
1607      ;; signal would have been raised.      ;; signal would have been raised.
1608      (with-float-traps-masked (:underflow :overflow)      ;;
1609        ;; Actually, for double-double-floats, we should probably
1610        ;; explicitly check for overflow instead of disabling the traps.
1611        ;; Why?  Because instead of overflow, double-double signals
1612        ;; invalid operation.
1613        (with-float-traps-masked (:underflow :overflow :invalid)
1614        (let ((rho (+ (square x) (square y))))        (let ((rho (+ (square x) (square y))))
1615          (declare (optimize (speed 3) (space 0)))          (declare (optimize (speed 3) (space 0)
1616                               (inhibit-warnings 3)))
1617          (cond ((and (or (float-nan-p rho)          (cond ((and (or (float-nan-p rho)
1618                          (float-infinity-p rho))                          (float-infinity-p rho))
1619                      (or (float-infinity-p (abs x))                      (or (float-infinity-p (abs x))
1620                          (float-infinity-p (abs y))))                          (float-infinity-p (abs y))))
1621                 (values (%make-double-double-float ext:double-float-positive-infinity 0w0)                 (values (%make-double-double-float ext:double-float-positive-infinity 0d0)
1622                         0))                         0))
1623                ((let ((threshold #.(/ least-positive-double-float                ((let ((threshold #.(/ least-positive-double-float
1624                                       double-float-epsilon))                                       double-float-epsilon))
# Line 1500  pi/4    11001001000011111101101010100010 Line 1640  pi/4    11001001000011111101101010100010
1640                 (values rho 0)))))))                 (values rho 0)))))))
1641    
1642  (defun dd-complex-sqrt (z)  (defun dd-complex-sqrt (z)
1643    "Principle square root of Z    _N"Principle square root of Z
1644    
1645  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1646    (declare (number z))    (declare (number z))
# Line 1516  Z may be any number, but the result is a Line 1656  Z may be any number, but the result is a
1656    
1657        (locally        (locally
1658            ;; space 0 to get maybe-inline functions inlined.            ;; space 0 to get maybe-inline functions inlined.
1659            (declare (optimize (speed 3) (space 0)))            (declare (optimize (speed 3) (space 0)
1660                                 (inhibit-warnings 3)))
1661    
1662          (if (not (float-nan-p x))          (if (not (float-nan-p x))
1663              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
# Line 1541  Z may be any number, but the result is a Line 1682  Z may be any number, but the result is a
1682          (complex eta nu)))))          (complex eta nu)))))
1683    
1684  (defun dd-complex-log-scaled (z j)  (defun dd-complex-log-scaled (z j)
1685    "Compute log(2^j*z).    _N"Compute log(2^j*z).
1686    
1687  This is for use with J /= 0 only when |z| is huge."  This is for use with J /= 0 only when |z| is huge."
1688    (declare (number z)    (declare (number z)
# Line 1559  This is for use with J /= 0 only when |z Line 1700  This is for use with J /= 0 only when |z
1700          (y (float (imagpart z) 1.0w0)))          (y (float (imagpart z) 1.0w0)))
1701      (multiple-value-bind (rho k)      (multiple-value-bind (rho k)
1702          (dd-cssqs z)          (dd-cssqs z)
1703        (declare (optimize (speed 3)))        (declare (optimize (speed 3)
1704                             (inhibit-warnings 3)))
1705        (let ((beta (max (abs x) (abs y)))        (let ((beta (max (abs x) (abs y)))
1706              (theta (min (abs x) (abs y))))              (theta (min (abs x) (abs y))))
1707          (complex (if (and (zerop k)          (complex (if (and (zerop k)
# Line 1575  This is for use with J /= 0 only when |z Line 1717  This is for use with J /= 0 only when |z
1717                   (atan y x))))))                   (atan y x))))))
1718    
1719  (defun dd-complex-log (z)  (defun dd-complex-log (z)
1720    "Log of Z = log |Z| + i * arg Z    _N"Log of Z = log |Z| + i * arg Z
1721    
1722  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1723    (declare (number z))    (declare (number z))
# Line 1587  Z may be any number, but the result is a Line 1729  Z may be any number, but the result is a
1729  ;; never 0 since we have positive and negative zeroes.  ;; never 0 since we have positive and negative zeroes.
1730    
1731  (defun dd-complex-atanh (z)  (defun dd-complex-atanh (z)
1732    "Compute atanh z = (log(1+z) - log(1-z))/2"    _N"Compute atanh z = (log(1+z) - log(1-z))/2"
1733    (declare (number z))    (declare (number z))
1734    (if (and (realp z) (< z -1))    (cond ((and (realp z) (< z -1))
1735        ;; atanh is continuous in quadrant III in this case.           ;; ATANH is continuous with quadrant III in this case.
1736        (dd-complex-atanh (complex z -0f0))           (dd-complex-atanh (complex z -0d0)))
1737        (let* ( ;; Constants          (t
1738               (theta (/ (sqrt most-positive-double-float) 4.0w0))           (let* ( ;; Constants
1739               (rho (/ 4.0w0 (sqrt most-positive-double-float)))                  (theta (/ (sqrt most-positive-double-float) 4.0w0))
1740               (half-pi dd-pi/2)                  (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1741               (rp (float (realpart z) 1.0w0))                  (half-pi dd-pi/2)
1742               (beta (float-sign rp 1.0w0))                  (rp (float (realpart z) 1.0w0))
1743               (x (* beta rp))                  (beta (float-sign rp 1.0w0))
1744               (y (* beta (- (float (imagpart z) 1.0w0))))                  (x (* beta rp))
1745               (eta 0.0w0)                  (y (* beta (- (float (imagpart z) 1.0w0))))
1746               (nu 0.0w0))                  (eta 0.0w0)
1747          ;; Shouldn't need this declare.                  (nu 0.0w0))
1748          (declare (double-double-float x y))             ;; Shouldn't need this declare.
1749          (locally             (declare (double-double-float x y))
1750              (declare (optimize (speed 3)))             (locally
1751            (cond ((or (> x theta)                 (declare (optimize (speed 3)
1752                       (> (abs y) theta))                                    (inhibit-warnings 3)))
1753                   ;; To avoid overflow...               (cond ((or (> x theta)
1754                   (setf nu (float-sign y half-pi))                          (> (abs y) theta))
1755                   ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),                      ;; To avoid overflow...
1756                   ;; which can cause overflow.  Arrange this computation so                      (setf nu (float-sign y half-pi))
1757                   ;; that it won't overflow.                      ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
1758                   (setf eta (let* ((x-bigger (> x (abs y)))                      ;; which can cause overflow.  Arrange this computation so
1759                                    (r (if x-bigger (/ y x) (/ x y)))                      ;; that it won't overflow.
1760                                    (d (+ 1.0d0 (* r r))))                      (setf eta (let* ((x-bigger (> x (abs y)))
1761                               (if x-bigger                                       (r (if x-bigger (/ y x) (/ x y)))
1762                                   (/ (/ x) d)                                       (d (+ 1.0d0 (* r r))))
1763                                   (/ (/ r y) d)))))                                  (if x-bigger
1764                  ((= x 1.0w0)                                      (/ (/ x) d)
1765                   ;; Should this be changed so that if y is zero, eta is set                                      (/ (/ r y) d)))))
1766                   ;; to +infinity instead of approx 176?  In any case                     ((= x 1.0w0)
1767                   ;; tanh(176) is 1.0d0 within working precision.                      ;; Should this be changed so that if y is zero, eta is set
1768                   (let ((t1 (+ 4w0 (square y)))                      ;; to +infinity instead of approx 176?  In any case
1769                         (t2 (+ (abs y) rho)))                      ;; tanh(176) is 1.0d0 within working precision.
1770                     (setf eta (dd-%log (/ (sqrt (sqrt t1))                      (let ((t1 (+ 4w0 (square y)))
1771                                           (sqrt t2))))                            (t2 (+ (abs y) rho)))
1772                     (setf nu (* 0.5d0                        (setf eta (dd-%log (/ (sqrt (sqrt t1))
1773                                 (float-sign y                                              (sqrt t2))))
1774                                             (+ half-pi (dd-%atan (* 0.5d0 t2))))))))                        (setf nu (* 0.5d0
1775                  (t                                    (float-sign y
1776                   (let ((t1 (+ (abs y) rho)))                                                (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1777                     ;; Normal case using log1p(x) = log(1 + x)                     (t
1778                     (setf eta (* 0.25d0                      (let ((t1 (+ (abs y) rho)))
1779                                  (dd-%log1p (/ (* 4.0d0 x)                        ;; Normal case using log1p(x) = log(1 + x)
1780                                                (+ (square (- 1.0d0 x))                        (setf eta (* 0.25d0
1781                                                   (square t1))))))                                     (dd-%log1p (/ (* 4.0d0 x)
1782                     (setf nu (* 0.5d0                                                   (+ (square (- 1.0d0 x))
1783                                 (dd-%atan2 (* 2.0d0 y)                                                      (square t1))))))
1784                                            (- (* (- 1.0d0 x)                        (setf nu (* 0.5d0
1785                                                  (+ 1.0d0 x))                                    (dd-%atan2 (* 2.0d0 y)
1786                                               (square t1))))))))                                               (- (* (- 1.0d0 x)
1787            (complex (* beta eta)                                                     (+ 1.0d0 x))
1788                     (- (* beta nu)))))))                                                  (square t1))))))))
1789                 (complex (* beta eta)
1790                          (- (* beta nu))))))))
1791    
1792  (defun dd-complex-tanh (z)  (defun dd-complex-tanh (z)
1793    "Compute tanh z = sinh z / cosh z"    _N"Compute tanh z = sinh z / cosh z"
1794    (declare (number z))    (declare (number z))
1795    (let ((x (float (realpart z) 1.0w0))    (let ((x (float (realpart z) 1.0w0))
1796          (y (float (imagpart z) 1.0w0)))          (y (float (imagpart z) 1.0w0)))
1797      (locally      (locally
1798          ;; space 0 to get maybe-inline functions inlined          ;; space 0 to get maybe-inline functions inlined
1799          (declare (optimize (speed 3) (space 0)))          (declare (optimize (speed 3) (space 0)
1800                               (inhibit-warnings 3)))
1801        (cond ((> (abs x)        (cond ((> (abs x)
1802                  #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)                  #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1803                  ;; This is more accurate under linux.                  ;; This is more accurate under linux.
# Line 1725  Z may be any number, but the result is a Line 1870  Z may be any number, but the result is a
1870    (complex (+ (realpart z) 1) (imagpart z)))    (complex (+ (realpart z) 1) (imagpart z)))
1871    
1872  (defun dd-complex-acos (z)  (defun dd-complex-acos (z)
1873    "Compute acos z = pi/2 - asin z    _N"Compute acos z = pi/2 - asin z
1874    
1875  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1876    (declare (number z))    (declare (number z))
# Line 1751  Z may be any number, but the result is a Line 1896  Z may be any number, but the result is a
1896                                              sqrt-1-z)))))))))                                              sqrt-1-z)))))))))
1897    
1898  (defun dd-complex-acosh (z)  (defun dd-complex-acosh (z)
1899    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))    _N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1900    
1901  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1902    (declare (number z))    (declare (number z))
# Line 1777  Z may be any number, but the result is a Line 1922  Z may be any number, but the result is a
1922    
1923    
1924  (defun dd-complex-asin (z)  (defun dd-complex-asin (z)
1925    "Compute asin z = asinh(i*z)/i    _N"Compute asin z = asinh(i*z)/i
1926    
1927  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1928    (declare (number z))    (declare (number z))
1929    (if (and (realp z) (> z 1))    (cond ((and (realp z) (> z 1))
1930        ;; asin is continuous in quadrant IV in this case.           (dd-complex-asin (complex z -0w0)))
1931        (dd-complex-asin (complex z -0f0))          (t
1932        (let* ((sqrt-1-z (complex-sqrt (1-z z)))           (let* ((sqrt-1-z (complex-sqrt (1-z z)))
1933               (sqrt-1+z (complex-sqrt (1+z z)))                  (sqrt-1+z (complex-sqrt (1+z z)))
1934               (den (realpart (* sqrt-1-z sqrt-1+z))))                  (den (realpart (* sqrt-1-z sqrt-1+z))))
1935          (cond ((zerop den)             (cond ((zerop den)
1936                 ;; Like below but we handle atan part ourselves.                    ;; Like below but we handle atan part ourselves.
1937                 (complex (if (minusp (float-sign den))                    ;; Must be sure to take into account the sign of
1938                              (- dd-pi/2)                    ;; (realpart z) and den!
1939                              dd-pi/2)                    (complex (if (minusp (* (float-sign (realpart z))
1940                     (asinh (imagpart (* (conjugate sqrt-1-z)                                            (float-sign den)))
1941                                         sqrt-1+z)))))                                 (- dd-pi/2)
1942                (t                                 dd-pi/2)
1943                 (with-float-traps-masked (:divide-by-zero)                             (asinh (imagpart (* (conjugate sqrt-1-z)
1944                   ;; We get a invalid operation here when z is real and |z| > 1.                                                 sqrt-1+z)))))
1945                   (complex (atan (/ (realpart z)                   (t
1946                                     (realpart (* sqrt-1-z sqrt-1+z))))                    (with-float-traps-masked (:divide-by-zero)
1947                            (asinh (imagpart (* (conjugate sqrt-1-z)                      ;; We get a invalid operation here when z is real and |z| > 1.
1948                                                sqrt-1+z))))))))))                      (complex (atan (/ (realpart z)
1949                                          den))
1950                                 (asinh (imagpart (* (conjugate sqrt-1-z)
1951                                                     sqrt-1+z)))))))))))
1952    
1953  (defun dd-complex-asinh (z)  (defun dd-complex-asinh (z)
1954    "Compute asinh z = log(z + sqrt(1 + z*z))    _N"Compute asinh z = log(z + sqrt(1 + z*z))
1955    
1956  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1957    (declare (number z))    (declare (number z))
# Line 1814  Z may be any number, but the result is a Line 1962  Z may be any number, but the result is a
1962               (- (realpart result)))))               (- (realpart result)))))
1963    
1964  (defun dd-complex-atan (z)  (defun dd-complex-atan (z)
1965    "Compute atan z = atanh (i*z) / i    _N"Compute atan z = atanh (i*z) / i
1966    
1967  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1968    (declare (number z))    (declare (number z))
# Line 1825  Z may be any number, but the result is a Line 1973  Z may be any number, but the result is a
1973               (- (realpart result)))))               (- (realpart result)))))
1974    
1975  (defun dd-complex-tan (z)  (defun dd-complex-tan (z)
1976    "Compute tan z = -i * tanh(i * z)    _N"Compute tan z = -i * tanh(i * z)
1977    
1978  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1979    (declare (number z))    (declare (number z))

Legend:
Removed from v.1.1.2.5  
changed lines
  Added in v.1.20

  ViewVC Help
Powered by ViewVC 1.1.5