/[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.17 - (hide annotations)
Thu Aug 2 18:18:18 2007 UTC (6 years, 8 months ago) by rtoy
Branch: MAIN
CVS Tags: snapshot-2007-09, snapshot-2008-08, snapshot-2008-09, sse2-packed-2008-11-12, snapshot-2008-05, snapshot-2008-06, snapshot-2008-07, snapshot-2008-01, snapshot-2008-02, snapshot-2008-03, sse2-base, sse2-packed-base, release-19f-pre1, snapshot-2008-12, snapshot-2008-11, release-19e, unicode-utf16-sync-2008-12, label-2009-03-16, release-19f-base, merge-sse2-packed, merge-with-19f, unicode-utf16-sync-2008-07, unicode-utf16-sync-2008-09, unicode-utf16-extfmts-sync-2008-12, snapshot-2008-04, unicode-utf16-sync-label-2009-03-16, RELEASE_19f, unicode-utf16-extfmts-pre-sync-2008-11, snapshot-2008-10, unicode-utf16-sync-2008-11, release-19e-pre1, release-19e-pre2, sse2-checkpoint-2008-10-01, sse2-merge-with-2008-11, sse2-merge-with-2008-10, unicode-utf16-string-support, release-19e-base, unicode-utf16-base, snapshot-2007-12, snapshot-2007-10, snapshot-2007-11, snapshot-2009-02, snapshot-2009-01, pre-telent-clx
Branch point for: RELEASE-19F-BRANCH, sse2-packed-branch, unicode-utf16-branch, release-19e-branch, sse2-branch, unicode-utf16-extfmt-branch
Changes since 1.16: +55 -28 lines
Add (INHIBIT-WARNINGS 3) to various bits of double-double float code.
I checked the notes and they all have to do with either boxing or
generic operations (usually from FLOOR).

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

  ViewVC Help
Powered by ViewVC 1.1.5