/[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 by rtoy, Thu Jun 29 01:28:02 2006 UTC revision 1.1.2.5 by rtoy, Thu Jun 29 20:49:12 2006 UTC
# Line 0  Line 1 
1    ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2    ;;;
3    ;;; **********************************************************************
4    ;;; This code was written as part of the CMU Common Lisp project at
5    ;;; Carnegie Mellon University, and has been placed in the public domain.
6    ;;;
7    (ext:file-comment
8      "$Header$")
9    ;;;
10    ;;; **********************************************************************
11    ;;;
12    ;;; This file contains all the irrational functions for double-double
13    ;;; float.  The algorithms and coefficients come from the Cephes math
14    ;;; library.  See http://www.moshier.net/#Cephes and
15    ;;; http://www.netlib.org/cephes/.
16    ;;;
17    ;;; Author: Raymond Toy.
18    ;;;
19    
20    (in-package "KERNEL")
21    
22    ;;;; Random constants, utility functions.
23    
24    #+nil
25    (defconstant max-log
26      1.1356523406294143949491931077970764891253w4
27      "log(most-positive-double-double-float)")
28    
29    #+nil
30    (defconstant min-log
31      -1.143276959615573793352782661133116431383730w4
32      "log(least-positive-double-double-float")
33    
34    (defconstant max-log
35      7.0978271289338399678773454114191w2
36      "log(most-positive-double-double-float)")
37    
38    (defconstant min-log
39      -7.4444007192138126231410729844608w2
40      "log(least-positive-double-double-float")
41    
42    
43    (defconstant loge2
44      0.6931471805599453094172321214581765680755001w0
45      "log(2)")
46    
47    (defconstant log2e
48      1.442695040888963407359924681001892137426646w0
49      "Log base 2 of e")
50    
51    (defconstant log2ea
52      4.4269504088896340735992468100189213742664595w-1
53      "log2(e)-1")
54    
55    (defconstant dd-pi
56      3.141592653589793238462643383279502884197169w0
57      "Pi")
58    
59    (defconstant dd-pi/2
60      1.570796326794896619231321691639751442098585w0
61      "Pi/2")
62    
63    (defconstant dd-pi/4
64      0.7853981633974483096156608458198757210492923w0
65      "Pi/4")
66    
67    ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +
68    ;; log2-c2 is log(2) to an accuracy greater than double-double-float.
69    (defconstant log2-c1
70      6.93145751953125w-1)
71    
72    (defconstant log2-c2
73      1.428606820309417232121458176568075500134w-6)
74    
75    (defconstant sqrt-1/2
76      0.7071067811865475244008443621048490392848w0
77      "Sqrt(2)")
78    
79    ;; Evaluate polynomial
80    (defun poly-eval (x coef)
81      (declare (type double-double-float x)
82               (type (simple-array double-double-float (*)) coef))
83      ;; y = C0 + C1*x + C2*x^2 + ...
84      ;;
85      ;; But coefficients are stored in reverse (descending powers) order:
86      ;; coef[0] = CN, ..., coef[N] = C0.
87      (let ((y 0w0))
88        (declare (type double-double-float y))
89        (loop for c across coef do
90             (setf y (+ (* y x)
91                        c)))
92        y))
93    
94    (defun poly-eval-1 (x coef)
95      (declare (type double-double-float x)
96               (type (simple-array double-double-float (*)) coef))
97      ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.
98      (let ((y 1w0))
99        (declare (type double-double-float y))
100        (loop for c across coef do
101             (setf y (+ (* y x)
102                        c)))
103        y))
104    
105    
106    ;;; dd-%expm1
107    ;;;
108    ;;; exp(x)-1
109    ;;; Range reduction is accomplished by separating the argument
110    ;;; into an integer k and fraction f such that
111    ;;;
112    ;;;     x    k  f
113    ;;;    e  = 2  e.
114    ;;;
115    ;;; An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
116    ;;; in the basic range [-0.5 ln 2, 0.5 ln 2].
117    ;;;
118    ;;;
119    ;;; ACCURACY:
120    ;;;
121    ;;;                      Relative error:
122    ;;; arithmetic   domain     # trials      peak         rms
123    ;;;    IEEE    -79,+MAXLOG    100,000     1.7e-34     4.5e-35
124    ;;;
125    
126    (let ((p (make-array 8 :element-type 'double-double-float
127                         :initial-contents '(
128                                             -4.888737542888633647784737721812546636240w-1
129                                             4.401308817383362136048032038528753151144w1
130                                             -1.716772506388927649032068540558788106762w3
131                                             4.578962475841642634225390068461943438441w4
132                                             -7.212432713558031519943281748462837065308w5
133                                             8.944630806357575461578107295909719817253w6
134                                             -5.722847283900608941516165725053359168840w7
135                                             2.943520915569954073888921213330863757240w8)))
136          (q (make-array 8 :element-type 'double-double-float
137                         :initial-contents '(
138                                             ;; 1.000000000000000000000000000000000000000w0
139                                             1.766112549341972444333352727998584753865w9
140                                             -7.848989743695296475743081255027098295771w8
141                                             1.615869009634292424463780387327037251069w8
142                                             -2.019684072836541751428967854947019415698w7
143                                             1.682912729190313538934190635536631941751w6
144                                             -9.615511549171441430850103489315371768998w4
145                                             3.697714952261803935521187272204485251835w3
146                                             -8.802340681794263968892934703309274564037w1)))
147          ;; ln 2^-114
148          (minarg -7.9018778583833765273564461846232128760607w1))
149      (declare (type double-double-float minarg))
150      ;; exp(x) - 1 = x + 0.5*x^2+x^3*P(x)/Q(x), for -log(2)/2 < x <
151      ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn
152      ;; above.  Theoretical peak relative error = 8.1e-36.
153      (defun dd-%expm1 (x)
154        "exp(x) - 1"
155        (declare (type double-double-float x)
156                 (optimize (speed 3)))
157        ;; Range reduction is accomplished by separating the argument
158        ;; into an integer k and fraction f such that
159        ;;
160        ;;     x    k  f
161        ;;    e  = 2  e.
162        ;;
163        ;; An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
164        ;; in the basic range [-0.5 ln 2, 0.5 ln 2].
165        (when (> x max-log)
166          (return-from dd-%expm1
167            (kernel:%make-double-double-float #.ext:double-float-positive-infinity
168                                              #.ext:double-float-positive-infinity)))
169        (when (< x minarg)
170          (return-from dd-%expm1 -1w0))
171    
172        ;; Express x = ln(2)*(k+remainder), remainder not exceeding 1/2
173        (let* ((xx (+ log2-c1 log2-c2))
174               (k (floor (+ 1/2 (/ (the (double-double-float * 710w0) x) xx))))
175               (px (coerce k 'double-double-float))
176               (qx 0w0))
177          (declare (type double-double-float xx px qx))
178          ;; remainder times ln 2
179          (decf x (* px log2-c1))
180          (decf x (* px log2-c2))
181    
182          ;; Approximate exp(remainder*ln(2))
183          #+nil
184          (setf px (* x
185                      (+ p0
186                         (* x
187                            (+ p1
188                               (* x
189                                  (+ p2
190                                     (* x
191                                        (+ p3
192                                           (* x
193                                              (+ p4
194                                                 (* x
195                                                    (+ p5
196                                                       (* x
197                                                          (+ p6
198                                                             (* x p7))))))))))))))))
199          #+nil
200          (setf qx (+ q0
201                      (* x
202                         (+ q1
203                            (* x
204                               (+ q2
205                                  (* x
206                                     (+ q3
207                                        (* x
208                                           (+ q4
209                                              (* x
210                                                 (+ q5
211                                                    (* x
212                                                       (+ q6
213                                                          (* x
214                                                             (+ x q7))))))))))))))))
215          (setf px (poly-eval x p))
216          (setf qx (poly-eval-1 x q))
217          (setf xx (* x x))
218          (setf qx (+ x (+ (* 0.5w0 xx)
219                           (/  (* xx px)
220                               qx))))
221          ;; exp(x) = exp(k*ln(2))*exp(remainder*ln(2)) = 2^k*exp(remainder*ln(2))
222          ;;
223          ;; We have qx = exp(remainder*ln(2))-1, so exp(x) - 1 =
224          ;; 2^k*(qx+1)-1 = 2^k*qx + 2^k - 1
225          (setf px (scale-float 1w0 k))
226          (+ (* px qx)
227             (- px 1)))))
228    
229    ;;; dd-%exp
230    ;;; exp(x)
231    ;;; Range reduction is accomplished by separating the argument
232    ;;; into an integer k and fraction f such that
233    ;;;
234    ;;;     x    k  f
235    ;;;    e  = 2  e.
236    ;;;
237    ;;; A Pade' form of degree 2/3 is used to approximate exp(f) - 1
238    ;;; in the basic range [-0.5 ln 2, 0.5 ln 2].
239    ;;;
240    ;;;
241    ;;; ACCURACY:
242    ;;;
243    ;;;                      Relative error:
244    ;;; arithmetic   domain     # trials      peak         rms
245    ;;;    IEEE      +-MAXLOG    100,000     2.6e-34     8.6e-35
246    ;;;
247    ;;;
248    ;;; Error amplification in the exponential function can be
249    ;;; a serious matter.  The error propagation involves
250    ;;; exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
251    ;;; which shows that a 1 lsb error in representing X produces
252    ;;; a relative error of X times 1 lsb in the function.
253    ;;; While the routine gives an accurate result for arguments
254    ;;; that are exactly represented by a long double precision
255    ;;; computer number, the result contains amplified roundoff
256    ;;; error for large arguments not exactly represented.
257    
258    (let ((p (make-array 5 :element-type 'double-double-float
259                         :initial-contents
260                         '(
261                           3.279723985560247033712687707263393506266w-10
262                           6.141506007208645008909088812338454698548w-7
263                           2.708775201978218837374512615596512792224w-4
264                           3.508710990737834361215404761139478627390w-2
265                           9.999999999999999999999999999999999998502w-1
266                           )))
267          (q (make-array 6 :element-type 'double-double-float
268                         :initial-contents
269                         '(
270                           2.980756652081995192255342779918052538681w-12
271                           1.771372078166251484503904874657985291164w-8
272                           1.504792651814944826817779302637284053660w-5
273                           3.611828913847589925056132680618007270344w-3
274                           2.368408864814233538909747618894558968880w-1
275                           2.000000000000000000000000000000000000150w0
276                           )))
277          ;; C1 + C2 = ln 2
278          (c1 (- log2-c1))
279          (c2 (- log2-c2)))
280      (declare (type double-double-float c1 c2))
281      ;; p and q are Pade coefficients for exp(x) - 1.  Theoretical peak
282      ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.
283      (defun dd-%exp (x)
284        (declare (type double-double-float x)
285                 (optimize (speed 3)))
286        (when (> x max-log)
287          (return-from dd-%exp
288            (kernel:%make-double-double-float #.ext:double-float-positive-infinity
289                                              #.ext:double-float-positive-infinity)))
290        (when (< x min-log)
291          (return-from dd-%exp 0w0))
292        ;; Express
293        ;;   exp(x) = exp(g)*2^n
294        ;;          = exp(g)*exp(n*log(2))
295        ;;          = exp(g + n * log2)
296        (let* ((n (floor (+ 0.5w0 (* x log2e))))
297               (px (coerce n 'double-double-float)))
298          (declare (type double-double-float px))
299          (incf x (* px c1))
300          (incf x (* px c2))
301          ;; Rational approx for exponential of fractional part:
302          ;; exp(x) = 1+2*x*p(x^2)/(q(x^2)-p(x^2))
303          ;;
304          ;; [The above comment seems not to match the code below.  But
305          ;; perhaps q(x) isn't what I think it is.]
306          (let ((xx (* x x)))
307            (setf px (* x (poly-eval xx p)))
308            (setf xx (poly-eval xx q))
309            (setf x (/ px (- xx px)))
310            (setf x (+ 1 x x))
311            (scale-float x n)))))
312    
313    
314    ;;; dd-%log
315    ;;; log(x), natural logarithm.
316    ;;; The argument is separated into its exponent and fractional
317    ;;; parts.  If the exponent is between -1 and +1, the logarithm
318    ;;; of the fraction is approximated by
319    ;;;
320    ;;;     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
321    ;;;
322    ;;; Otherwise, setting  z = 2(x-1)/x+1),
323    ;;;
324    ;;;     log(x) = z + z**3 P(z)/Q(z).
325    ;;;
326    ;;;
327    ;;;
328    ;;; ACCURACY:
329    ;;;
330    ;;;                      Relative error:
331    ;;; arithmetic   domain     # trials      peak         rms
332    ;;;    IEEE   exp(+-MAXLOGL) 36,000      9.5e-35     4.1e-35
333    ;;;
334    
335    (let ((p (make-array 13 :element-type 'double-double-float
336                         :initial-contents
337                         '(1.538612243596254322971797716843006400388w-6
338                           4.998469661968096229986658302195402690910w-1
339                           2.321125933898420063925789532045674660756w1
340                           4.114517881637811823002128927449878962058w2
341                           3.824952356185897735160588078446136783779w3
342                           2.128857716871515081352991964243375186031w4
343                           7.594356839258970405033155585486712125861w4
344                           1.797628303815655343403735250238293741397w5
345                           2.854829159639697837788887080758954924001w5
346                           3.007007295140399532324943111654767187848w5
347                           2.014652742082537582487669938141683759923w5
348                           7.771154681358524243729929227226708890930w4
349                           1.313572404063446165910279910527789794488w4)))
350          (q (make-array 12 :element-type 'double-double-float
351                         :initial-contents
352                         '(
353                           ;; 1.000000000000000000000000000000000000000w0
354                           4.839208193348159620282142911143429644326w1
355                           9.104928120962988414618126155557301584078w2
356                           9.147150349299596453976674231612674085381w3
357                           5.605842085972455027590989944010492125825w4
358                           2.248234257620569139969141618556349415120w5
359                           6.132189329546557743179177159925690841200w5
360                           1.158019977462989115839826904108208787040w6
361                           1.514882452993549494932585972882995548426w6
362                           1.347518538384329112529391120390701166528w6
363                           7.777690340007566932935753241556479363645w5
364                           2.626900195321832660448791748036714883242w5
365                           3.940717212190338497730839731583397586124w4)))
366          (r (make-array 6 :element-type 'double-double-float
367                         :initial-contents
368                         '(-8.828896441624934385266096344596648080902w-1
369                           8.057002716646055371965756206836056074715w1
370                           -2.024301798136027039250415126250455056397w3
371                           2.048819892795278657810231591630928516206w4
372                           -8.977257995689735303686582344659576526998w4
373                           1.418134209872192732479751274970992665513w5
374                           )))
375          (s (make-array 6 :element-type 'double-double-float
376                         :initial-contents
377                         '( ;; 1.000000000000000000000000000000000000000w0
378                           -1.186359407982897997337150403816839480438w2
379                           3.998526750980007367835804959888064681098w3
380                           -5.748542087379434595104154610899551484314w4
381                           4.001557694070773974936904547424676279307w5
382                           -1.332535117259762928288745111081235577029w6
383                           1.701761051846631278975701529965589676574w6))))
384      ;; p and q are coefficients in the expansion for log(1+x) = x -
385      ;; x^2/2 + x^3*p(x)/q(x), where 1/sqrt(2) <= 1+x < sqrt(2).
386      ;; Theoretical peak relative error = 5.3e-37, relative peak error
387      ;; spread = 2.3e-14.
388      ;;
389      ;; R and S are coefficients for log(x) = z + z^3*r(z^2)/s(z^2),
390      ;; where z = 2*(x-1)/(x+1), where 1/sqrt(2) <= x < sqrt(2).
391      ;; Theoretical peak relative error = 1.1e-35, relative peak error
392      ;; spread 1.1e-9.
393      ;;
394      ;;
395      (defun dd-%log (x)
396        (declare (type double-double-float x)
397                 (optimize (speed 3)))
398        ;; Separate mantissa from exponent
399        (multiple-value-bind (x e)
400            (decode-float x)
401          (declare (type double-double-float x)
402                   (type double-float-exponent e))
403          (let ((z 0w0)
404                (y 0w0))
405            (declare (type double-double-float z y))
406            ;; Logarithm using log(x) = z + z^3*P(z^2)/Q(z^2); where z = 2*(x-1)/x + 1;
407            (cond ((or (> e 2)
408                       (< e -2))
409                   (cond ((< x sqrt-1/2)
410                          ;; 2*(2*x-1)/(2*x+1)
411                          (decf e)
412                          (setf z (- x 0.5w0))
413                          (setf y (+ (* 0.5w0 z) 0.5w0)))
414                         (t
415                          ;; 2*(x-1)/(x+1)
416                          (setf z (- x 0.5w0))
417                          (decf z 0.5w0)
418                          (setf y (+ (* 0.5w0 x) 0.5w0))))
419                   (setf x (/ z y))
420                   (setf z (* x x))
421                   (setf z (* x (/ (* z (poly-eval z r))
422                                     (poly-eval-1 z s))))
423                   (incf z (* e log2-c2))
424                   (incf z x)
425                   (incf z (* e log2-c1)))
426                  (t
427                   ;; Log using log(1+x) = x - 0.5*x^2+x^3*p(x)/q(x)
428                   (cond ((< x sqrt-1/2)
429                          (decf e)
430                          (setf x (- (scale-float x 1) 1)))
431                         (t
432                          (decf x)))
433                   (setf z (* x x))
434                   (setf y (* x (* z (/ (poly-eval x p)
435                                        (poly-eval-1 x q)))))
436                   (incf y (* e log2-c2))
437                   ;; z = y - 0.5*z
438                   (setf z (- y (scale-float z -1)))
439                   (incf z x)
440                   (incf z (* e log2-c1))
441                   ))))))
442    
443    ;;; dd-%log1p(x)
444    ;;; log(1+x)
445    ;;;
446    ;;; The argument 1+x is separated into its exponent and fractional
447    ;;; parts.  If the exponent is between -1 and +1, the logarithm
448    ;;; of the fraction is approximated by
449    ;;;
450    ;;;     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
451    ;;;
452    ;;; Otherwise, setting  z = 2(x-1)/x+1),
453    ;;;
454    ;;;     log(x) = z + z^3 P(z)/Q(z).
455    ;;;
456    ;;;
457    ;;;
458    ;;; ACCURACY:
459    ;;;
460    ;;;                      Relative error:
461    ;;; arithmetic   domain     # trials      peak         rms
462    ;;;    IEEE      -1, 8       100000      1.9e-34     4.3e-35
463    
464    (let ((p (make-array 13 :element-type 'double-double-float
465                         :initial-contents '(
466                                             1.538612243596254322971797716843006400388w-6
467                                             4.998469661968096229986658302195402690910w-1
468                                             2.321125933898420063925789532045674660756w1
469                                             4.114517881637811823002128927449878962058w2
470                                             3.824952356185897735160588078446136783779w3
471                                             2.128857716871515081352991964243375186031w4
472                                             7.594356839258970405033155585486712125861w4
473                                             1.797628303815655343403735250238293741397w5
474                                             2.854829159639697837788887080758954924001w5
475                                             3.007007295140399532324943111654767187848w5
476                                             2.014652742082537582487669938141683759923w5
477                                             7.771154681358524243729929227226708890930w4
478                                             1.313572404063446165910279910527789794488w4
479                                             )))
480          (q (make-array 12 :element-type 'double-double-float
481                         :initial-contents '(
482                                             ;; 1.000000000000000000000000000000000000000w0
483                                             4.839208193348159620282142911143429644326w1
484                                             9.104928120962988414618126155557301584078w2
485                                             9.147150349299596453976674231612674085381w3
486                                             5.605842085972455027590989944010492125825w4
487                                             2.248234257620569139969141618556349415120w5
488                                             6.132189329546557743179177159925690841200w5
489                                             1.158019977462989115839826904108208787040w6
490                                             1.514882452993549494932585972882995548426w6
491                                             1.347518538384329112529391120390701166528w6
492                                             7.777690340007566932935753241556479363645w5
493                                             2.626900195321832660448791748036714883242w5
494                                             3.940717212190338497730839731583397586124w4
495                                             )))
496          (r (make-array 6 :element-type 'double-double-float
497                         :initial-contents '(-8.828896441624934385266096344596648080902w-1
498                                             8.057002716646055371965756206836056074715w1
499                                             -2.024301798136027039250415126250455056397w3
500                                             2.048819892795278657810231591630928516206w4
501                                             -8.977257995689735303686582344659576526998w4
502                                             1.418134209872192732479751274970992665513w5
503                                             )))
504          (s (make-array 6 :element-type 'double-double-float
505                         :initial-contents '(
506                                             ;; 1.000000000000000000000000000000000000000w0
507                                             -1.186359407982897997337150403816839480438w2
508                                             3.998526750980007367835804959888064681098w3
509                                             -5.748542087379434595104154610899551484314w4
510                                             4.001557694070773974936904547424676279307w5
511                                             -1.332535117259762928288745111081235577029w6
512                                             1.701761051846631278975701529965589676574w6
513                                             )))
514    )
515      ;; Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
516      ;; 1/sqrt(2) <= 1+x < sqrt(2)
517      ;; Theoretical peak relative error = 5.3e-37,
518      ;; relative peak error spread = 2.3e-14
519      (defun dd-%log1p (xm1)
520        (declare (type double-double-float xm1))
521        (let ((x (+ xm1 1))
522              (z 0w0)
523              (y 0w0))
524          (declare (type double-double-float x y z)
525                   (optimize (speed 3)))
526          (multiple-value-bind (x e)
527              (decode-float x)
528            (declare (type double-double-float x)
529                     (type double-float-exponent e))
530            (cond ((or (> e 2)
531                       (< e -2))
532                   ;; Log using log(x) = z + z^3*P(z^2)/Q(z^2)
533                   (cond ((< x sqrt-1/2)
534                          (decf e)
535                          (setf z (- x 0.5w0))
536                          (setf y (+ 0.5w0 (* 0.5w0 z))))
537                         (t
538                          (setf z (- x 0.5w0))
539                          (decf z 0.5w0)
540                          (setf y (+ 0.5w0 (* 0.5w0 x)))))
541                   (setf x (/ z y))
542                   (setf z (* x x))
543                   (setf z (* x (/ (* z (poly-eval z r))
544                                   (poly-eval-1 z s))))
545                   (incf z (* e log2-c2))
546                   (incf z x)
547                   (incf z (* e log2-c1)))
548                  (t
549                   ;; Log using log(1+x) = x - 0.5*x^2 + x^3*p(x)/q(x)
550                   (cond ((< x sqrt-1/2)
551                          (decf e)
552                          (if (/= e 0)
553                              (setf x (- (* 2 x) 1))
554                              (setf x xm1)))
555                         (t
556                          (if (/= e 0)
557                              (decf x 1)
558                              (setf x xm1))))
559                   (setf z (* x x))
560                   (setf y (* x (/ (* z (poly-eval x p))
561                                   (poly-eval-1 x q))))
562                   (incf y (* e log2-c2))
563                   (setf z (- y (* 0.5w0 z)))
564                   (incf z x)
565                   (incf z (* e log2-c1))))))))
566    
567    (let ((P (make-array 6 :element-type 'double-double-float
568                         :initial-contents '(
569                                             1.622194395724068297909052717437740288268w3
570                                             1.124862584587770079742188354390171794549w6
571                                             3.047548980769660162696832999871894196102w8
572                                             3.966215348072348368191433063260384329745w10
573                                             2.375869584584371194838551715348965605295w12
574                                             6.482835792103233269752264509192030816323w13
575                                             )))
576          (Q (make-array 6 :element-type 'double-double-float
577                         :initial-contents '(
578                                             ;; 1.000000000000000000000000000000000000000w0 */
579                                             -9.101683853129357776079049616394849086007w2
580                                             4.486400519836461218634448973793765123186w5
581                                             -1.492531313030440305095318968983514314656w8
582                                             3.457771488856930054902696708717192082887w10
583                                             -5.193289868803472640225483235513427062460w12
584                                             3.889701475261939961851358705515223019890w14))))
585      (defun dd-%sinh (x)
586        (declare (type double-double-float x)
587                 (optimize (speed 3)))
588        (let ((a (abs x)))
589          (declare (type double-double-float a))
590          (cond ((> a 1)
591                 (setf a (dd-%exp a))
592                 (setf a (- (* 0.5w0 a) (/ 0.5w0 a)))
593                 (if (< x 0)
594                     (- a)
595                     a))
596                (t
597                 (setf a (* a a))
598                 (let ((pp (poly-eval a p))
599                       (qq (poly-eval-1 a q)))
600                   (+ x (* x a
601                           (/ pp qq)))))))))
602    
603    #+nil
604    (defun dd-%sinh (x)
605      (declare (type double-double-float x))
606      (let ((a (abs x)))
607        (declare (type double-double-float a))
608        (cond ((> a 1)
609               (setf a (dd-%exp a))
610               (setf a (- (* 0.5w0 a) (/ 0.5w0 a)))
611               (if (< x 0)
612                   (- a)
613                   a))
614              (t
615               ;; Use sinh(x) = 1/2*(D(x)+D(x)/(1+D(x))), here D(x) =
616               ;; exp(x)-1.
617               (let ((d (dd-%expm1 x)))
618                 (* 0.5w0 (+ d (/ d (+ 1 d)))))))))
619    
620    (defun dd-%cosh (x)
621      (declare (type double-double-float x))
622      (let ((y (dd-%exp x)))
623        (scale-float (+ y (/ y)) -1)))
624    
625    (let ((P (make-array 6 :element-type 'double-double-float
626                         :initial-contents '(-6.505693197948351084912624750702492767503w-6
627                                             -9.804083860188429726356968570322356183383w-1
628                                             -5.055287638900473250703725789725376004355w2
629                                             -7.307477148073823966594990496301416814519w4
630                                             -3.531606586182691280701462523692471322688w6
631                                             -4.551377146142783468144190926206842300707w7
632                                             )))
633          (Q (make-array 5 :element-type 'double-double-float
634                         :initial-contents '(
635                                             ;; 1.000000000000000000000000000000000000000w0 */
636                                             5.334865598460027935735737253027154828002w2
637                                             8.058475607422391042912151298751537172870w4
638                                             4.197073523796142343374222405869721575491w6
639                                             6.521134551226147545983467868553677881771w7
640                                             1.365413143842835040443257277862054198329w8))))
641      (defun dd-%tanh (x)
642        (declare (type double-double-float x)
643                 (optimize (speed 3)))
644        (let ((z (abs x)))
645          (declare (type double-double-float z))
646          (cond ((> z (* 0.5w0 max-log))
647                 (if (> x 0)
648                     1w0
649                     -1w0))
650                ((> z 0.625w0)
651                 (let ((s (dd-%exp (* 2 z))))
652                   (setf z (- 1 (/ 2 (+ s 1))))
653                   (if (minusp x)
654                       (- z)
655                       z)))
656                (t
657                 (let* ((s (* x x)))
658                   (declare (optimize speed))
659                   (setf z (/ (poly-eval s p)
660                              (poly-eval-1 s q)))
661                   (setf z (* x s z))
662                   (+ x z)))))))
663    
664    (let ((P (make-array 10 :element-type 'double-double-float
665                         :initial-contents '(
666                                             -9.217569843805850417698565442251656375681w-1
667                                             5.321929116410615470118183794063211260728w1
668                                             -9.139522976807685333981548145417830690552w2
669                                             7.204314536952949779101646454146682033772w3
670                                             -3.097809640165146436529075324081668598891w4
671                                             7.865376554210973897486215630898496100534w4
672                                             -1.211716814094785128366087489224821937203w5
673                                             1.112669508789123834670923967462068457013w5
674                                             -5.600242872292477863751728708249167956542w4
675                                             1.188901082233997739779618679364295772810w4
676                                             )))
677          (Q (make-array 10 :element-type 'double-double-float
678                         :initial-contents '(
679                                             ;; 1.000000000000000000000000000000000000000w0 */
680                                             -6.807348436010016270202879229504392062418w1
681                                             1.386763299649315831625106608182196351693w3
682                                             -1.310805752656879543134785263832907269320w4
683                                             6.872174720355764193772953852564737816928w4
684                                             -2.181008360536226513009076189881617939380w5
685                                             4.362736119602298592874941767284979857248w5
686                                             -5.535251007539393347687001489396152923502w5
687                                             4.321594849688346708841188057241308805551w5
688                                             -1.894075056489862952285849974761239845873w5
689                                             3.566703246701993219338856038092901974725w4
690                                             ))))
691      (defun dd-%atanh (x)
692        (declare (type double-double-float x)
693                 (optimize (speed 3)))
694        (cond ((minusp x)
695               (- (the double-double-float (dd-%atanh (- x)))))
696              ((< x 1w-12)
697               x)
698              ((< x  0.5w0)
699               (let ((z (* x x)))
700                 (+ x (* x z (/ (poly-eval z p)
701                                (poly-eval-1 z q))))))
702              (t
703               (* 0.5w0 (dd-%log (/ (+ 1 x)
704                                   (- 1 x))))))))
705    
706    
707    (let ((P (make-array 9 :element-type 'double-double-float
708                         :initial-contents '(
709                                             -8.104404283317298189545629468767571317688w-1
710                                             -4.954206127425209147110732546633675599008w1
711                                             -8.438175619831548439550086251740438689853w2
712                                             -6.269710069245210459536983820505214648057w3
713                                             -2.418935474493501382372711518024193326434w4
714                                             -5.208121780431312783866941311277024486498w4
715                                             -6.302755086521614763280617114866439227971w4
716                                             -4.003566436224198252093684987323233921339w4
717                                             -1.037690841528359305134494613113086980551w4
718                                             )))
719          (Q (make-array 9 :element-type 'double-double-float
720                         :initial-contents '(
721                                             ;; 1.000000000000000000000000000000000000000w0 */
722                                             8.175806439951395194771977809279448392548w1
723                                             1.822215299975696008284027212745010251320w3
724                                             1.772040003462901790853111853838978236828w4
725                                             9.077625379864046240143413577745818879353w4
726                                             2.675554475070211205153169988669677418808w5
727                                             4.689758557916492969463473819426544383586w5
728                                             4.821923684550711724710891114802924039911w5
729                                             2.682316388947175963642524537892687560973w5
730                                             6.226145049170155830806967678679167550122w4))))
731      (defun dd-%asinh (x)
732        (declare (type double-double-float x)
733                 (optimize (speed 3)))
734        (cond ((minusp x)
735               (- (the double-double-float (dd-%asinh (- x)))))
736              #+nil
737              ((> x 1w10)
738               (+ loge2 (dd-%log x)))
739              ((< x 0.5w0)
740               (let* ((z (* x x))
741                      (a (* z (/ (poly-eval z p)
742                                 (poly-eval-1 z q)))))
743                 (+ (* a x) x)))
744              (t
745               (dd-%log (+ x (sqrt (1+ (* x x)))))))))
746    
747    (let ((P (make-array 10 :element-type 'double-double-float
748                         :initial-contents '(
749                                             1.895467874386341763387398084072833727168w-1
750                                             6.443902084393244878979969557171256604767w1
751                                             3.914593556594721458616408528941154205393w3
752                                             9.164040999602964494412169748897754668733w4
753                                             1.065909694792026382660307834723001543839w6
754                                             6.899169896709615182428217047370629406305w6
755                                             2.599781868717579447900896150777162652518w7
756                                             5.663733059389964024656501196827345337766w7
757                                             6.606302846870644033621560858582696134512w7
758                                             3.190482951215438078279772140481195200593w7
759                                             )))
760          (Q (make-array 9 :element-type 'double-double-float
761                         :initial-contents '(
762                                             ;; 1.000000000000000000000000000000000000000w0 */
763                                             1.635418024331924674147953764918262009321w2
764                                             7.290983678312632723073455563799692165828w3
765                                             1.418207894088607063257675159183397062114w5
766                                             1.453154285419072886840913424715826321357w6
767                                             8.566841438576725234955968880501739464425w6
768                                             3.003448667795089562511136059766833630017w7
769                                             6.176592872899557661256383958395266919654w7
770                                             6.872176426138597206811541870289420510034w7
771                                             3.190482951215438078279772140481195226621w7
772                                             ))))
773      (defun dd-%acosh (x)
774        (declare (type double-double-float x)
775                 (optimize (speed 3)))
776        (cond ((> x 1w17)
777               (+ loge2 (dd-%log x)))
778              (t
779               (let ((z (- x 1)))
780                 (cond ((< z 0.5w0)
781                        (* (sqrt (* 2 z))
782                           (/ (poly-eval z p)
783                              (poly-eval-1 z q))))
784                       (t
785                        (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))
786    
787    (let ((P (make-array 10 :element-type 'double-double-float
788                         :initial-contents '(
789                                             -8.067112765482705313585175280952515549833w-1
790                                             4.845649797786849136525020822000172350977w1
791                                             -8.510195404865297879959793548843395926847w2
792                                             6.815196841370292688574521445731895826485w3
793                                             -2.967135182120339728996157454994675519735w4
794                                             7.612250656518818109652985996692466409670w4
795                                             -1.183360579752620455689557157684221905030w5
796                                             1.095432262510413338755837156377401348063w5
797                                             -5.554124580991113991999636773382495788705w4
798                                             1.187132626694762543537732514905488896985w4
799                                             )))
800          (Q (make-array 10 :element-type 'double-double-float
801                         :initial-contents '(
802                                             ;;  1.000000000000000000000000000000000000000w0 */
803                                             -8.005471061732009595694099899234272342478w1
804                                             1.817324228942812880965069608562483918025w3
805                                             -1.867017317425756524289537002141956583706w4
806                                             1.048196619402464497478959760337779705622w5
807                                             -3.527040897897253459022458866536165564103w5
808                                             7.426302422018858001691440351763370029242w5
809                                             -9.863068411558756277454631976667880674474w5
810                                             8.025654653926121907774766642393757364326w5
811                                             -3.653000557802254281954969843055623398839w5
812                                             7.122795760168575261226395089432959614179w4
813                                             ))))
814      (defun dd-%asin (x)
815        (declare (type double-double-float x)
816                 (optimize (speed 3)))
817        (cond ((minusp x)
818               (- (the double-double-float (dd-%asin (- x)))))
819              #+nil
820              ((< x 1w-8)
821               x)
822              (t
823               (let ((flag nil)
824                     (z 0w0)
825                     (zz 0w0))
826                 (declare (type double-double-float z zz))
827                 (cond ((> x 0.5w0)
828                        (setf zz (- 0.5w0 x))
829                        (setf zz (scale-float (+ zz 0.5w0) -1))
830                        (setf z (sqrt zz))
831                        (setf flag t))
832                       (t
833                        (setf z x)
834                        (setf zz (* z z))
835                        (setf flag nil)))
836                 (let ((p (* zz (/ (poly-eval zz p)
837                                   (poly-eval-1 zz q)))))
838                   (setf z (+ (* z p) z))
839                   (when flag
840                     (setf z (+ z z))
841                     (setf z (- dd-pi/2 z)))
842                   z))))))
843    
844    (defun dd-%acos (x)
845      (declare (type double-double-float x)
846               (optimize (speed 3)))
847      (cond ((< x -0.5w0)
848             (- dd-pi
849                (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))
850            ((> x 0.5w0)
851             (* 2 (dd-%asin (sqrt (* 0.5w0 (- 1w0 x))))))
852            (t
853             (- dd-pi/2 (dd-%asin x)))))
854    
855    (let ((P (make-array 9 :element-type 'double-double-float
856                         :initial-contents '(
857                                             -6.635810778635296712545011270011752799963w-4
858                                             -8.768423468036849091777415076702113400070w-1
859                                             -2.548067867495502632615671450650071218995w1
860                                             -2.497759878476618348858065206895055957104w2
861                                             -1.148164399808514330375280133523543970854w3
862                                             -2.792272753241044941703278827346430350236w3
863                                             -3.696264445691821235400930243493001671932w3
864                                             -2.514829758941713674909996882101723647996w3
865                                             -6.880597774405940432145577545328795037141w2
866                                             )))
867          (Q (make-array 8 :element-type 'double-double-float
868                         :initial-contents '(
869                                             ;; 1.000000000000000000000000000000000000000w0 */
870                                             3.566239794444800849656497338030115886153w1
871                                             4.308348370818927353321556740027020068897w2
872                                             2.494680540950601626662048893678584497900w3
873                                             7.928572347062145288093560392463784743935w3
874                                             1.458510242529987155225086911411015961174w4
875                                             1.547394317752562611786521896296215170819w4
876                                             8.782996876218210302516194604424986107121w3
877                                             2.064179332321782129643673263598686441900w3
878                                             )))
879    
880          ;; tan( 3*pi/8 )
881          (T3P8  2.414213562373095048801688724209698078569672w0)
882          ;; tan( pi/8 )
883          (TP8 0.414213562373095048801688724209698078569672w0))
884      (declare (type double-double-float t3p8 tp8))
885      (defun dd-%atan (x)
886        (declare (type double-double-float x)
887                 (optimize (speed 3)))
888        (when (minusp x)
889          (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))
890        ;; range reduction
891        (let ((y 0w0))
892          (declare (type double-double-float y))
893          (cond ((> x t3p8)
894                 (setf y dd-pi/2)
895                 (setf x (/ -1 x)))
896                ((> x tp8)
897                 (setf y dd-pi/4)
898                 (setf x (/ (- x 1)
899                            (+ x 1))))
900                (t
901                 (setf y 0w0)))
902          ;; Rational form in x^2
903          (let ((z (* x x)))
904            (setf y (+ y
905                       (* (/ (poly-eval z p)
906                             (poly-eval-1 z q))
907                          z x)
908                       x))
909            y))))
910    
911    (defun dd-%atan2 (y x)
912      (declare (type double-double-float x y)
913               (optimize (speed 3)))
914      (let ((code 0)
915            (w 0w0))
916        (declare (type (integer 0 3) code)
917                 (type double-double-float w))
918        (when (minusp x)
919          (setf code 2))
920        (when (minusp y)
921          (setf code (logior code 1)))
922        (when (zerop x)
923          (unless (zerop (logand code 1))
924            (return-from dd-%atan2 (- dd-pi/2)))
925          (when (zerop y)
926            (return-from dd-%atan2 0w0))
927          (return-from dd-%atan2 dd-pi/2))
928        (when (zerop y)
929          (return-from dd-%atan2
930            (if (zerop (logand code 2))
931                0w0
932                dd-pi)))
933        (setf w (ecase code
934                  (0 0w0)
935                  (1 0w0)
936                  (2 dd-pi)
937                  (3 (- dd-pi))))
938    
939        (+ w (dd-%atan (/ y x)))))
940    
941    
942    ;;
943    ;; Here are the fractional digits of pi, in hex.
944    ;;
945    ;;  243f6a8885 a308d31319 8a2e037073 44a4093822 299f31d008 2efa98ec4e 6c89452821 e638d01377 be5466cf34 e90c6cc0ac
946    ;; 29b7c97c50 dd3f84d5b5 b547091792 16d5d98979 fb1bd1310b a698dfb5ac 2ffd72dbd0 1adfb7b8e1 afed6a267e 96ba7c9045
947    ;;
948    ;; We want to express pi/4 in 3 parts such that the sum of the parts
949    ;; is exactly the value of pi to the higher precision.  For
950    ;; double-double, we want 53 bits in each part.
951    ;;
952    ;; pi/4 in hex is
953    ;;
954    ;; C90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411
955    ;;
956    ;;
957    #||
958    (let* ((frac #xC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411)
959           (len (integer-length frac)))
960      (format t "len = ~D bits~%" len)
961      (format t "pi/4    ~,,' ,106:b~%" (ldb (byte (* 4 106) (- len (* 4 106))) frac)))
962    ->
963    pi/4    11001001000011111101101010100010001000010110100011000 01000110100110001001100011001100010100010111000000011 01110000011100110100010010100100000010010011100000100 01000101001100111110011000111010000000010000010111011
964    
965    ||#
966    ;; d1 = (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106)
967    ;; d2 = (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106))
968    ;; d3 = (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106))
969    ;; d4 = (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106))
970    
971    ;; sin(x) = x + x^3 P(x^2)
972    ;; Theoretical peak relative error = 5.6e-39
973    ;; relative peak error spread = 1.7e-9
974    
975    (defconstant sincof
976      (make-array 12 :element-type 'double-double-float
977                  :initial-contents
978                  '(
979                    6.410290407010279602425714995528976754871w-26
980                    -3.868105354403065333804959405965295962871w-23
981                    1.957294039628045847156851410307133941611w-20
982                    -8.220635246181818130416407184286068307901w-18
983                    2.811457254345322887443598804951004537784w-15
984                    -7.647163731819815869711749952353081768709w-13
985                    1.605904383682161459812515654720205050216w-10
986                    -2.505210838544171877505034150892770940116w-8
987                    2.755731922398589065255731765498970284004w-6
988                    -1.984126984126984126984126984045294307281w-4
989                    8.333333333333333333333333333333119885283w-3
990                    -1.666666666666666666666666666666666647199w-1
991                    )))
992    
993    ;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2))
994    ;; Theoretical peak relative error = 2.1e-37,
995    ;; relative peak error spread = 1.4e-8
996    
997    (defconstant coscof
998      (make-array 11 :element-type 'double-double-float
999                  :initial-contents
1000                  '(
1001                    1.601961934248327059668321782499768648351w-24
1002                    -8.896621117922334603659240022184527001401w-22
1003                    4.110317451243694098169570731967589555498w-19
1004                    -1.561920696747074515985647487260202922160w-16
1005                    4.779477332386900932514186378501779328195w-14
1006                    -1.147074559772972328629102981460088437917w-11
1007                    2.087675698786809897637922200570559726116w-9
1008                    -2.755731922398589065255365968070684102298w-7
1009                    2.480158730158730158730158440896461945271w-5
1010                    -1.388888888888888888888888888765724370132w-3
1011                    4.166666666666666666666666666666459301466w-2
1012                    )))
1013    
1014    (defconstant dp1
1015      (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106))
1016    
1017    (defconstant dp2
1018      (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106)))
1019    
1020    (defconstant dp3
1021      (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106)))
1022    
1023    (defconstant dp4
1024      (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
1025    
1026    (defun dd-%sin (x)
1027      (declare (type double-double-float x))
1028      (when (minusp x)
1029        (return-from dd-%sin (- (the double-double-float (dd-%sin (- x))))))
1030      ;; y = integer part of x/(pi/4).
1031      (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1032             (z (scale-float y -4)))
1033        (declare (type double-double-float y z))
1034        (setf z (float (floor z) 1w0))      ; integer part of y/8
1035        (setf z (- y (scale-float z 4)))    ; y - 16*(y/16)
1036    
1037        (let ((j (truncate z))
1038              (sign 1))
1039          (unless (zerop (logand j 1))
1040            (incf j)
1041            (incf y))
1042          (setf j (logand j 7))
1043    
1044          (when (> j 3)
1045            (setf sign (- sign))
1046            (decf j 4))
1047    
1048          ;; Extended precision modular arithmetic
1049          (setf z (- (- (- x (* y dp1))
1050                        (* y dp2))
1051                     (* y dp3)))
1052          (let ((zz (* z z)))
1053            (if (or (= j 1)
1054                    (= j 2))
1055                (setf y (+ (- 1 (scale-float zz -1))
1056                           (* zz zz (poly-eval zz coscof))))
1057                (setf y (+ z (* z (* zz (poly-eval zz sincof))))))
1058            (if (< sign 0)
1059                (- y)
1060                y)))))
1061    
1062    (defun dd-%cos (x)
1063      (declare (type double-double-float x)
1064               (optimize (speed 3)))
1065      (when (minusp x)
1066        (return-from dd-%cos (dd-%cos (- x))))
1067      ;; y = integer part of x/(pi/4).
1068      (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1069             (z (scale-float y -4)))
1070        (declare (type double-double-float y z))
1071        (setf z (float (floor z) 1w0))      ; integer part of y/8
1072        (setf z (- y (scale-float z 4)))    ; y - 16*(y/16)
1073    
1074        (let ((i (truncate z))
1075              (j 0)
1076              (sign 1))
1077          (declare (type (integer 0 7) j)
1078                   (type (integer -1 1) sign))
1079          (unless (zerop (logand i 1))
1080            (incf i)
1081            (incf y))
1082          (setf j (logand i 7))
1083    
1084          (when (> j 3)
1085            (setf sign (- sign))
1086            (decf j 4))
1087          (when (> j 1)
1088            (setf sign (- sign)))
1089    
1090          ;; Extended precision modular arithmetic.  This is basically
1091          ;; computing x - y*(pi/4) accurately so that |z| < pi/4.
1092          (setf z (- (- (- x (* y dp1))
1093                        (* y dp2))
1094                     (* y dp3)))
1095          (let ((zz (* z z)))
1096            (if (or (= j 1)
1097                    (= j 2))
1098                (setf y (+ z (* z (* zz (poly-eval zz sincof)))))
1099                (setf y (+ (- 1 (scale-float zz -1))
1100                           (* zz (poly-eval zz coscof) zz))))
1101            (if (< sign 0)
1102                (- y)
1103                y)))))
1104    
1105    (let ((P (make-array 6 :element-type 'double-double-float
1106                         :initial-contents
1107                         '(
1108                           -9.889929415807650724957118893791829849557w-1
1109                           1.272297782199996882828849455156962260810w3
1110                           -4.249691853501233575668486667664718192660w5
1111                           5.160188250214037865511600561074819366815w7
1112                           -2.307030822693734879744223131873392503321w9
1113                           2.883414728874239697964612246732416606301w10
1114                           )))
1115          (Q (make-array 6 :element-type 'double-double-float
1116                         :initial-contents
1117                         '(
1118                           ;; 1.000000000000000000000000000000000000000w0 */
1119                           -1.317243702830553658702531997959756728291w3
1120                           4.529422062441341616231663543669583527923w5
1121                           -5.733709132766856723608447733926138506824w7
1122                           2.758476078803232151774723646710890525496w9
1123                           -4.152206921457208101480801635640958361612w10
1124                           8.650244186622719093893836740197250197602w10
1125                           ))))
1126      (defun dd-tancot (xx cotflag)
1127        (declare (type double-double-float xx)
1128                 (optimize (speed 3)))
1129        (let ((x 0w0)
1130              (sign 1))
1131          (declare (type double-double-float x)
1132                   (type (integer -1 1) sign))
1133          (cond ((minusp xx)
1134                 (setf x (- xx))
1135                 (setf sign -1))
1136                (t
1137                 (setf x xx)))
1138          (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1139                 (z (scale-float y -4))
1140                 (j 0))
1141            (declare (type double-double-float y z)
1142                     (type fixnum j))
1143            (setf z (float (floor z) 1w0))
1144            (setf z (- y (scale-float z 4)))
1145    
1146            (setf j (truncate z))
1147    
1148            (unless (zerop (logand j 1))
1149              (incf j)
1150              (incf y))
1151    
1152            (setf z (- (- (- x (* y dp1))
1153                          (* y dp2))
1154                       (* y dp3)))
1155            (let ((zz (* z z)))
1156              (if (> zz 1w-40)
1157                  (setf y (+ z
1158                             (* z (* zz (/ (poly-eval zz p)
1159                                           (poly-eval-1 zz q))))))
1160                  (setf y z))
1161              (if (not (zerop (logand j 2)))
1162                  (if cotflag
1163                      (setf y (- y))
1164                      (setf y (/ -1 y)))
1165                  (if cotflag
1166                      (setf y (/ y))))
1167              (if (< sign 0)
1168                  (- y)
1169                  y))))))
1170    
1171    (defun dd-%tan (x)
1172      (declare (type double-double-float x))
1173      (dd-tancot x nil))
1174    
1175    ;;; dd-%log2
1176    ;;; Base 2 logarithm.
1177    
1178    (let ((P (make-array 13 :element-type 'double-double-float
1179                         :initial-contents
1180                         '(
1181                           1.538612243596254322971797716843006400388w-6
1182                           4.998469661968096229986658302195402690910w-1
1183                           2.321125933898420063925789532045674660756w1
1184                           4.114517881637811823002128927449878962058w2
1185                           3.824952356185897735160588078446136783779w3
1186                           2.128857716871515081352991964243375186031w4
1187                           7.594356839258970405033155585486712125861w4
1188                           1.797628303815655343403735250238293741397w5
1189                           2.854829159639697837788887080758954924001w5
1190                           3.007007295140399532324943111654767187848w5
1191                           2.014652742082537582487669938141683759923w5
1192                           7.771154681358524243729929227226708890930w4
1193                           1.313572404063446165910279910527789794488w4
1194                           )))
1195          (Q (make-array 12 :element-type 'double-double-float
1196                         :initial-contents
1197                         '(
1198                           ;; 1.000000000000000000000000000000000000000w0
1199                           4.839208193348159620282142911143429644326w1
1200                           9.104928120962988414618126155557301584078w2
1201                           9.147150349299596453976674231612674085381w3
1202                           5.605842085972455027590989944010492125825w4
1203                           2.248234257620569139969141618556349415120w5
1204                           6.132189329546557743179177159925690841200w5
1205                           1.158019977462989115839826904108208787040w6
1206                           1.514882452993549494932585972882995548426w6
1207                           1.347518538384329112529391120390701166528w6
1208                           7.777690340007566932935753241556479363645w5
1209                           2.626900195321832660448791748036714883242w5
1210                           3.940717212190338497730839731583397586124w4
1211                           )))
1212          ;; Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
1213          ;; where z = 2(x-1)/(x+1)
1214          ;; 1/sqrt(2) <= x < sqrt(2)
1215          ;; Theoretical peak relative error = 1.1e-35,
1216          ;; relative peak error spread 1.1e-9
1217          (R (make-array 6 :element-type 'double-double-float
1218                         :initial-contents
1219                         '(
1220                           -8.828896441624934385266096344596648080902w-1
1221                           8.057002716646055371965756206836056074715w1
1222                           -2.024301798136027039250415126250455056397w3
1223                           2.048819892795278657810231591630928516206w4
1224                           -8.977257995689735303686582344659576526998w4
1225                           1.418134209872192732479751274970992665513w5
1226                           )))
1227          (S (make-array 6 :element-type 'double-double-float
1228                         :initial-contents
1229                         '(
1230                           ;; 1.000000000000000000000000000000000000000w0 */
1231                           -1.186359407982897997337150403816839480438w2
1232                           3.998526750980007367835804959888064681098w3
1233                           -5.748542087379434595104154610899551484314w4
1234                           4.001557694070773974936904547424676279307w5
1235                           -1.332535117259762928288745111081235577029w6
1236                           1.701761051846631278975701529965589676574w6
1237                           ))))
1238      (defun dd-%log2 (x)
1239        (declare (type double-double-float x)
1240                 (optimize (speed 3)))
1241        (multiple-value-bind (x e)
1242            (decode-float x)
1243          (declare (type double-double-float x)
1244                   (type double-float-exponent e))
1245          (let ((z 0w0)
1246                (y 0w0))
1247            (declare (type double-double-float z y))
1248            (cond ((or (> e 2)
1249                       (< e -2))
1250                   (cond ((< x sqrt-1/2)
1251                          ;; 2*(2*x-1)/(2*x+1)
1252                          (decf e)
1253                          (setf z (- x 0.5w0))
1254                          (setf y (+ (* 0.5w0 z) 0.5w0)))
1255                         (t
1256                          ;; 2*(x-1)/(x+1)
1257                          (setf z (- x 0.5w0))
1258                          (decf z 0.5w0)
1259                          (setf y (+ (* 0.5w0 z) 0.5w0))))
1260                   (setf x (/ z y))
1261                   (setf z (* x x))
1262                   (setf y (* x (/ (* z (poly-eval z r))
1263                                   (poly-eval-1 z s)))))
1264                  (t
1265                   (cond ((< x sqrt-1/2)
1266                          (decf e)
1267                          (setf x (- (scale-float x 1) 1)))
1268                         (t
1269                          (decf x)))
1270                   (setf z (* x x))
1271                   (setf y (* x (/ (* z (poly-eval x p))
1272                                   (poly-eval-1 x q))))
1273                   (decf y (scale-float z -1))))
1274            ;; Multiply log of fraction by log2(e) and base 2 exponent by 1
1275            ;;
1276            ;; This sequence of operations is critical
1277            (setf z (* y log2ea))
1278            (setf z (+ z (* x log2ea)))
1279            (setf z (+ z y))
1280            (setf z (+ z x))
1281            (setf z (+ z e))
1282            z))))
1283    
1284    ;;; dd-%exp2
1285    ;;; 2^x
1286    
1287    (let ((P (make-array 5 :element-type 'double-double-float
1288                         :initial-contents
1289                         '(
1290                           1.587171580015525194694938306936721666031w2
1291                           6.185032670011643762127954396427045467506w5
1292                           5.677513871931844661829755443994214173883w8
1293                           1.530625323728429161131811299626419117557w11
1294                           9.079594442980146270952372234833529694788w12
1295                           )))
1296          (Q (make-array 5 :element-type 'double-double-float
1297                         :initial-contents
1298                         '(
1299                           ;; 1.000000000000000000000000000000000000000w0 */
1300                           1.236602014442099053716561665053645270207w4
1301                           2.186249607051644894762167991800811827835w7
1302                           1.092141473886177435056423606755843616331w10
1303                           1.490560994263653042761789432690793026977w12
1304                           2.619817175234089411411070339065679229869w13
1305                           ))))
1306      (defun dd-%exp2 (x)
1307        (declare (type double-double-float x)
1308                 (optimize (speed 3)))
1309        (when (>= x 1024w0)
1310          (return-from dd-%exp2
1311            (%make-double-double-float ext:double-float-positive-infinity
1312                                       ext:double-float-positive-infinity)))
1313        (when (<= x -1024w0)
1314          (return-from dd-%exp2 0w0))
1315        (multiple-value-bind (n x)
1316            (floor (the (double-double-float -1024w0 1024w0) x))
1317          (declare (type double-double-float x))
1318          (let* ((xx (* x x))
1319                 (px (* x (poly-eval xx p))))
1320            (setf x (/ px (- (poly-eval-1 xx q) px)))
1321            (setf x (+ 1 (scale-float x 1)))
1322            (scale-float x n)))))
1323    
1324    ;;; dd-%powil
1325    ;;; x^n where n is an integer
1326    
1327    (defun dd-%powil (x nn)
1328      (declare (type double-double-float x)
1329               (fixnum nn)
1330               (optimize (speed 3)))
1331      (when (zerop x)
1332        (return-from dd-%powil
1333          (cond ((zerop nn)
1334                 1w0)
1335                ((minusp nn)
1336                 (%make-double-double-float ext:double-float-positive-infinity
1337                                            ext:double-float-positive-infinity))
1338                (t
1339                 0w0))))
1340      (when (zerop nn)
1341        (return-from dd-%powil 1w0))
1342    
1343      (let ((asign 0)
1344            (sign 0)
1345            (n nn))
1346        (declare (type (integer -1 0) asign sign)
1347                 (fixnum n))
1348        (when (minusp x)
1349          (setf asign -1)
1350          (setf x (- x)))
1351        (cond ((minusp nn)
1352               (setf sign -1)
1353               (setf n (- nn)))
1354              (t
1355               (setf sign 0)
1356               (setf n nn)))
1357        ;; Overflow detection
1358    
1359        ;; Calculate approximate log of answer
1360        (multiple-value-bind (s lx)
1361            (decode-float x)
1362          (declare (type double-double-float s)
1363                   (type double-float-exponent lx))
1364          (let ((e (* n (1- lx))))
1365            (cond ((or (zerop e)
1366                       (> e 64)
1367                       (< e -64))
1368                   (setf s (/ (- s 7.0710678118654752w-1)
1369                              (+ s 7.0710678118654752w-1)))
1370                   (setf s (* (+ (- (* 2.9142135623730950w0 s)
1371                                    0.5w0)
1372                                 lx)
1373                              nn log2e)))
1374                  (t
1375                   (setf s (* loge2 e))))
1376            (when (> s max-log)
1377              ;; Overflow.  What to do?
1378              (error "Overflow"))
1379            (when (< s min-log)
1380              (return-from dd-%powil 0w0))
1381    
1382            ;; Handle denormal answer
1383            (when (< s (- 2 max-log))
1384              (setf x (/ x))
1385              (setf sign 0))
1386    
1387            ;; First bit of the power
1388            (let ((y 0w0))
1389              (declare (type double-double-float y))
1390              (cond ((zerop (logand n 1))
1391                     (setf y 1w0)
1392                     (setf asign 0))
1393                    (t
1394                     (setf y x)))
1395              (let ((w x))
1396                (declare (type double-double-float w))
1397                (setf n (ash n -1))
1398                (loop while (not (zerop n))
1399                   do
1400                   (setf w (* w w))
1401                   (unless (zerop (logand n 1))
1402                     (setf y (* y w)))
1403                   (setf n (ash n -1))))
1404              ;; Odd power of a negative number
1405              (unless (zerop asign)
1406                (setf y (- y)))
1407    
1408              (unless (zerop sign)
1409                (setf y (/ y)))
1410              y)))))
1411    
1412    
1413    ;;; dd-real-pow
1414    ;;; x^y, for x and y real, and real result.
1415    
1416    (defun dd-real-pow (x y)
1417      (declare (type double-double-float x y)
1418               (optimize (speed 3)))
1419      (let ((nflg 0)
1420            (w (floor y)))
1421        ;; nflg = 1 if x < 0 raised to integer power
1422        (when (and (= w y)
1423                   (< (abs w) 32768))
1424          (return-from dd-real-pow (dd-%powil x w)))
1425    
1426        (when (<= x 0)
1427          (cond ((zerop x)
1428                 (if (zerop y)
1429                     ;; 0^0
1430                     (return-from dd-real-pow 1w0)
1431                     ;; 0^y
1432                     (return-from dd-real-pow 0w0)))
1433                (t
1434                 (when (/= w y)
1435                   ;; noninteger power of negative number
1436                   (let ((p (the double-double-float (dd-real-pow (abs x) y)))
1437                         (y*pi (* y dd-pi)))
1438                     (return-from dd-real-pow (complex (* p (dd-%cos y*pi))
1439                                                       (* p (dd-%sin y*pi))))))
1440    
1441                 ;; For negative x, find out if the integer exponent is odd or even.
1442                 (let ((w (scale-float y -1)))
1443                   (declare (type double-double-float w))
1444                   (setf w (ffloor w))
1445                   (setf w (scale-float w 1))
1446                   (when (/= w y)
1447                     (setf nflg 1))
1448                   (setf x (abs x))))))
1449        (let ((z (dd-%log2 x)))
1450          (declare (type double-double-float z))
1451          (setf z (dd-%exp2 (* y z)))
1452          (unless (zerop nflg)
1453            (setf z (- z)))
1454          z)))
1455    
1456    (defun dd-%pow (x y)
1457      (declare (type double-double-float x y))
1458      (dd-real-pow x y))
1459    
1460    
1461    ;; These are essentially the same as in irrat.lisp, but very slightly
1462    ;; modified to operate on double-double-floats.
1463    (defun dd-cssqs (z)
1464      ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1465      ;; result is r + i*k, where k is an integer.
1466    
1467      ;; Save all FP flags
1468      (let ((x (float (realpart z) 1w0))
1469            (y (float (imagpart z) 1w0)))
1470        ;; Would this be better handled using an exception handler to
1471        ;; catch the overflow or underflow signal?  For now, we turn all
1472        ;; traps off and look at the accrued exceptions to see if any
1473        ;; signal would have been raised.
1474        (with-float-traps-masked (:underflow :overflow)
1475          (let ((rho (+ (square x) (square y))))
1476            (declare (optimize (speed 3) (space 0)))
1477            (cond ((and (or (float-nan-p rho)
1478                            (float-infinity-p rho))
1479                        (or (float-infinity-p (abs x))
1480                            (float-infinity-p (abs y))))
1481                   (values (%make-double-double-float ext:double-float-positive-infinity 0w0)
1482                           0))
1483                  ((let ((threshold #.(/ least-positive-double-float
1484                                         double-float-epsilon))
1485                         (traps (ldb vm::float-sticky-bits
1486                                     (vm:floating-point-modes))))
1487                     ;; Overflow raised or (underflow raised and rho <
1488                     ;; lambda/eps)
1489                     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1490                         (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1491                              (< rho threshold))))
1492                   ;; If we're here, neither x nor y are infinity and at
1493                   ;; least one is non-zero.. Thus logb returns a nice
1494                   ;; integer.
1495                   (let ((k (- (logb-finite (max (abs x) (abs y))))))
1496                     (values (+ (square (scalb x k))
1497                                (square (scalb y k)))
1498                             (- k))))
1499                  (t
1500                   (values rho 0)))))))
1501    
1502    (defun dd-complex-sqrt (z)
1503      "Principle square root of Z
1504    
1505    Z may be any number, but the result is always a complex."
1506      (declare (number z))
1507      (multiple-value-bind (rho k)
1508          (dd-cssqs z)
1509        (declare (type (or (member 0w0) (double-double-float 0w0)) rho)
1510                 (type fixnum k))
1511        (let ((x (float (realpart z) 1.0w0))
1512              (y (float (imagpart z) 1.0w0))
1513              (eta 0w0)
1514              (nu 0w0))
1515          (declare (type double-double-float x y eta nu))
1516    
1517          (locally
1518              ;; space 0 to get maybe-inline functions inlined.
1519              (declare (optimize (speed 3) (space 0)))
1520    
1521            (if (not (float-nan-p x))
1522                (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1523    
1524            (cond ((oddp k)
1525                   (setf k (ash k -1)))
1526                  (t
1527                   (setf k (1- (ash k -1)))
1528                   (setf rho (+ rho rho))))
1529    
1530            (setf rho (scalb (sqrt rho) k))
1531    
1532            (setf eta rho)
1533            (setf nu y)
1534    
1535            (when (/= rho 0d0)
1536              (when (not (float-infinity-p (abs nu)))
1537                (setf nu (/ (/ nu rho) 2d0)))
1538              (when (< x 0d0)
1539                (setf eta (abs nu))
1540                (setf nu (float-sign y rho))))
1541            (complex eta nu)))))
1542    
1543    (defun dd-complex-log-scaled (z j)
1544      "Compute log(2^j*z).
1545    
1546    This is for use with J /= 0 only when |z| is huge."
1547      (declare (number z)
1548               (fixnum j))
1549      ;; The constants t0, t1, t2 should be evaluated to machine
1550      ;; precision.  In addition, Kahan says the accuracy of log1p
1551      ;; influences the choices of these constants but doesn't say how to
1552      ;; choose them.  We'll just assume his choices matches our
1553      ;; implementation of log1p.
1554      (let ((t0 #.(/ 1 (sqrt 2.0w0)))
1555            (t1 1.2w0)
1556            (t2 3w0)
1557            (ln2 #.(log 2w0))
1558            (x (float (realpart z) 1.0w0))
1559            (y (float (imagpart z) 1.0w0)))
1560        (multiple-value-bind (rho k)
1561            (dd-cssqs z)
1562          (declare (optimize (speed 3)))
1563          (let ((beta (max (abs x) (abs y)))
1564                (theta (min (abs x) (abs y))))
1565            (complex (if (and (zerop k)
1566                              (< t0 beta)
1567                              (or (<= beta t1)
1568                                  (< rho t2)))
1569                         (/ (dd-%log1p (+ (* (- beta 1.0d0)
1570                                             (+ beta 1.0d0))
1571                                          (* theta theta)))
1572                            2d0)
1573                         (+ (/ (log rho) 2d0)
1574                            (* (+ k j) ln2)))
1575                     (atan y x))))))
1576    
1577    (defun dd-complex-log (z)
1578      "Log of Z = log |Z| + i * arg Z
1579    
1580    Z may be any number, but the result is always a complex."
1581      (declare (number z))
1582      (dd-complex-log-scaled z 0))
1583    
1584    ;; Let us note the following "strange" behavior.  atanh 1.0d0 is
1585    ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1586    ;; reason for the imaginary part is caused by the fact that arg i*y is
1587    ;; never 0 since we have positive and negative zeroes.
1588    
1589    (defun dd-complex-atanh (z)
1590      "Compute atanh z = (log(1+z) - log(1-z))/2"
1591      (declare (number z))
1592      (if (and (realp z) (< z -1))
1593          ;; atanh is continuous in quadrant III in this case.
1594          (dd-complex-atanh (complex z -0f0))
1595          (let* ( ;; Constants
1596                 (theta (/ (sqrt most-positive-double-float) 4.0w0))
1597                 (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1598                 (half-pi dd-pi/2)
1599                 (rp (float (realpart z) 1.0w0))
1600                 (beta (float-sign rp 1.0w0))
1601                 (x (* beta rp))
1602                 (y (* beta (- (float (imagpart z) 1.0w0))))
1603                 (eta 0.0w0)
1604                 (nu 0.0w0))
1605            ;; Shouldn't need this declare.
1606            (declare (double-double-float x y))
1607            (locally
1608                (declare (optimize (speed 3)))
1609              (cond ((or (> x theta)
1610                         (> (abs y) theta))
1611                     ;; To avoid overflow...
1612                     (setf nu (float-sign y half-pi))
1613                     ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
1614                     ;; which can cause overflow.  Arrange this computation so
1615                     ;; that it won't overflow.
1616                     (setf eta (let* ((x-bigger (> x (abs y)))
1617                                      (r (if x-bigger (/ y x) (/ x y)))
1618                                      (d (+ 1.0d0 (* r r))))
1619                                 (if x-bigger
1620                                     (/ (/ x) d)
1621                                     (/ (/ r y) d)))))
1622                    ((= x 1.0w0)
1623                     ;; Should this be changed so that if y is zero, eta is set
1624                     ;; to +infinity instead of approx 176?  In any case
1625                     ;; tanh(176) is 1.0d0 within working precision.
1626                     (let ((t1 (+ 4w0 (square y)))
1627                           (t2 (+ (abs y) rho)))
1628                       (setf eta (dd-%log (/ (sqrt (sqrt t1))
1629                                             (sqrt t2))))
1630                       (setf nu (* 0.5d0
1631                                   (float-sign y
1632                                               (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1633                    (t
1634                     (let ((t1 (+ (abs y) rho)))
1635                       ;; Normal case using log1p(x) = log(1 + x)
1636                       (setf eta (* 0.25d0
1637                                    (dd-%log1p (/ (* 4.0d0 x)
1638                                                  (+ (square (- 1.0d0 x))
1639                                                     (square t1))))))
1640                       (setf nu (* 0.5d0
1641                                   (dd-%atan2 (* 2.0d0 y)
1642                                              (- (* (- 1.0d0 x)
1643                                                    (+ 1.0d0 x))
1644                                                 (square t1))))))))
1645              (complex (* beta eta)
1646                       (- (* beta nu)))))))
1647    
1648    (defun dd-complex-tanh (z)
1649      "Compute tanh z = sinh z / cosh z"
1650      (declare (number z))
1651      (let ((x (float (realpart z) 1.0w0))
1652            (y (float (imagpart z) 1.0w0)))
1653        (locally
1654            ;; space 0 to get maybe-inline functions inlined
1655            (declare (optimize (speed 3) (space 0)))
1656          (cond ((> (abs x)
1657                    #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1658                    ;; This is more accurate under linux.
1659                    #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1660                                              (%log most-positive-double-float)) 4d0))
1661                 (complex (float-sign x)
1662                          (float-sign y)))
1663                (t
1664                 (let* ((tv (dd-%tan y))
1665                        (beta (+ 1.0d0 (* tv tv)))
1666                        (s (sinh x))
1667                        (rho (sqrt (+ 1.0d0 (* s s)))))
1668                   (if (float-infinity-p (abs tv))
1669                       (complex (/ rho s)
1670                                (/ tv))
1671                       (let ((den (+ 1.0d0 (* beta s s))))
1672                         (complex (/ (* beta rho s)
1673                                     den)
1674                                  (/ tv den))))))))))
1675    
1676    ;; Kahan says we should only compute the parts needed.  Thus, the
1677    ;; realpart's below should only compute the real part, not the whole
1678    ;; complex expression.  Doing this can be important because we may get
1679    ;; spurious signals that occur in the part that we are not using.
1680    ;;
1681    ;; However, we take a pragmatic approach and just use the whole
1682    ;; expression.
1683    
1684    ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1685    ;; it's the conjugate of the square root or the square root of the
1686    ;; conjugate.  This needs to be checked.
1687    
1688    ;; I checked.  It doesn't matter because (conjugate (sqrt z)) is the
1689    ;; same as (sqrt (conjugate z)) for all z.  This follows because
1690    ;;
1691    ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1692    ;;
1693    ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1694    ;;
1695    ;; and these two expressions are equal if and only if arg conj z =
1696    ;; -arg z, which is clearly true for all z.
1697    
1698    ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1699    ;; complex, the real is converted to a complex before performing the
1700    ;; operation.  However, Kahan says in this paper (pg 176):
1701    ;;
1702    ;; (iii) Careless handling can turn infinity or the sign of zero into
1703    ;;       misinformation that subsequently disappears leaving behind
1704    ;;       only a plausible but incorrect result.  That is why compilers
1705    ;;       must not transform z-1 into z-(1+i*0), as we have seen above,
1706    ;;       nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1707    ;;       subsequent logarithm or square root produce a non-zero
1708    ;;       imaginary part whose sign is opposite to what was intended.
1709    ;;
1710    ;; The interesting examples are too long and complicated to reproduce
1711    ;; here.  We refer the reader to his paper.
1712    ;;
1713    ;; The functions below are intended to handle the cases where a real
1714    ;; is mixed with a complex and we don't want CL complex contagion to
1715    ;; occur..
1716    
1717    (declaim (inline 1+z 1-z z-1 z+1))
1718    (defun dd-1+z (z)
1719      (complex (+ 1 (realpart z)) (imagpart z)))
1720    (defun dd-1-z (z)
1721      (complex (- 1 (realpart z)) (- (imagpart z))))
1722    (defun dd-z-1 (z)
1723      (complex (- (realpart z) 1) (imagpart z)))
1724    (defun dd-z+1 (z)
1725      (complex (+ (realpart z) 1) (imagpart z)))
1726    
1727    (defun dd-complex-acos (z)
1728      "Compute acos z = pi/2 - asin z
1729    
1730    Z may be any number, but the result is always a complex."
1731      (declare (number z))
1732      (if (and (realp z) (> z 1))
1733          ;; acos is continuous in quadrant IV in this case.
1734          (complex-acos (complex z -0f0))
1735          (let ((sqrt-1+z (complex-sqrt (1+z z)))
1736                (sqrt-1-z (complex-sqrt (1-z z))))
1737            (cond ((zerop (realpart sqrt-1+z))
1738                   ;; Same as below, but we compute atan ourselves (because we
1739                   ;; have atan +/- infinity).
1740                   (complex
1741                            (if (minusp (float-sign (* (realpart sqrt-1-z)
1742                                                       (realpart sqrt-1+z))))
1743                                (- dd-pi)
1744                                dd-pi)
1745                            (asinh (imagpart (* (conjugate sqrt-1+z)
1746                                                sqrt-1-z)))))
1747                  (t
1748                   (complex (* 2 (atan (/ (realpart sqrt-1-z)
1749                                          (realpart sqrt-1+z))))
1750                            (asinh (imagpart (* (conjugate sqrt-1+z)
1751                                                sqrt-1-z)))))))))
1752    
1753    (defun dd-complex-acosh (z)
1754      "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1755    
1756    Z may be any number, but the result is always a complex."
1757      (declare (number z))
1758      (let* ((sqrt-z-1 (complex-sqrt (z-1 z)))
1759             (sqrt-z+1 (complex-sqrt (z+1 z))))
1760        ;; We need to handle the case where real part of sqrt-z+1 is zero,
1761        ;; because division by zero with double-double-floats doesn't
1762        ;; produce infinity.
1763        (cond ((zerop (realpart sqrt-z+1))
1764               ;; Same as below, but we compute atan ourselves (because we
1765               ;; have atan +/- infinity).
1766               (complex (asinh (realpart (* (conjugate sqrt-z-1)
1767                                            sqrt-z+1)))
1768                        (if (minusp (float-sign (* (imagpart sqrt-z-1)
1769                                                   (realpart sqrt-z+1))))
1770                            (- dd-pi)
1771                            dd-pi)))
1772              (t
1773               (complex (asinh (realpart (* (conjugate sqrt-z-1)
1774                                            sqrt-z+1)))
1775                        (* 2 (atan (/ (imagpart sqrt-z-1)
1776                                      (realpart sqrt-z+1)))))))))
1777    
1778    
1779    (defun dd-complex-asin (z)
1780      "Compute asin z = asinh(i*z)/i
1781    
1782    Z may be any number, but the result is always a complex."
1783      (declare (number z))
1784      (if (and (realp z) (> z 1))
1785          ;; asin is continuous in quadrant IV in this case.
1786          (dd-complex-asin (complex z -0f0))
1787          (let* ((sqrt-1-z (complex-sqrt (1-z z)))
1788                 (sqrt-1+z (complex-sqrt (1+z z)))
1789                 (den (realpart (* sqrt-1-z sqrt-1+z))))
1790            (cond ((zerop den)
1791                   ;; Like below but we handle atan part ourselves.
1792                   (complex (if (minusp (float-sign den))
1793                                (- dd-pi/2)
1794                                dd-pi/2)
1795                       (asinh (imagpart (* (conjugate sqrt-1-z)
1796                                           sqrt-1+z)))))
1797                  (t
1798                   (with-float-traps-masked (:divide-by-zero)
1799                     ;; We get a invalid operation here when z is real and |z| > 1.
1800                     (complex (atan (/ (realpart z)
1801                                       (realpart (* sqrt-1-z sqrt-1+z))))
1802                              (asinh (imagpart (* (conjugate sqrt-1-z)
1803                                                  sqrt-1+z))))))))))
1804    
1805    (defun dd-complex-asinh (z)
1806      "Compute asinh z = log(z + sqrt(1 + z*z))
1807    
1808    Z may be any number, but the result is always a complex."
1809      (declare (number z))
1810      ;; asinh z = -i * asin (i*z)
1811      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1812             (result (complex-asin iz)))
1813        (complex (imagpart result)
1814                 (- (realpart result)))))
1815    
1816    (defun dd-complex-atan (z)
1817      "Compute atan z = atanh (i*z) / i
1818    
1819    Z may be any number, but the result is always a complex."
1820      (declare (number z))
1821      ;; atan z = -i * atanh (i*z)
1822      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1823             (result (complex-atanh iz)))
1824        (complex (imagpart result)
1825                 (- (realpart result)))))
1826    
1827    (defun dd-complex-tan (z)
1828      "Compute tan z = -i * tanh(i * z)
1829    
1830    Z may be any number, but the result is always a complex."
1831      (declare (number z))
1832      ;; tan z = -i * tanh(i*z)
1833      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1834             (result (complex-tanh iz)))
1835        (complex (imagpart result)
1836                 (- (realpart result)))))

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

  ViewVC Help
Powered by ViewVC 1.1.5