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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5